lds.hpp

Go to the documentation of this file.

#pragma once
 
#include <array>
#include <atomic>
#include <cmath>
#include <concepts>
#include <cstddef>
#include <iterator>
#include <limits>
 
#ifndef M_PI
#    define M_PI 3.14159265358979323846264338327950288
#endif
 
namespace ldsgen {
 
    constexpr const auto TWO_PI = 2.0 * M_PI;
 
    // Constants for magic numbers
    constexpr unsigned int MAX_REVERSE_BITS = 64;
    constexpr double MAPPING_FACTOR = 2.0;
 
    namespace detail {
 
        template <typename T, typename Table>
        constexpr auto vdc_digit_sum(unsigned long n, unsigned long base, const Table& weights)
            -> T {
            T reslt{};
            std::size_t idx = 0;
            while (n != 0) {
                const auto remainder = n % base;
                n /= base;
                reslt += static_cast<T>(remainder) * weights[idx];
                ++idx;
            }
            return reslt;
        }
 
    }  // namespace detail
 
    template <typename Generator, typename Value> class GeneratorIterator {
        Generator* gen;
        unsigned long index;
 
      public:
        using iterator_category = std::input_iterator_tag;
        using value_type = Value;
        using difference_type = std::ptrdiff_t;
        using pointer = const value_type*;
        using reference = value_type;
 
        explicit GeneratorIterator(Generator* g = nullptr, unsigned long idx = 0)
            : gen{g}, index{idx} {}
 
        auto operator*() const -> Value {
            if (gen) {
                auto temp_idx = gen->get_index();
                gen->reseed(index);
                auto value = gen->pop();
                gen->reseed(temp_idx);
                return value;
            }
            return Value{};
        }
 
        auto operator++() -> GeneratorIterator& {
            ++index;
            return *this;
        }
 
        auto operator++(int) -> GeneratorIterator {
            auto temp = *this;
            ++index;
            return temp;
        }
 
        auto operator==(const GeneratorIterator& other) const -> bool {
            return index == other.index;
        }
 
        auto operator!=(const GeneratorIterator& other) const -> bool {
            return index != other.index;
        }
 
        [[nodiscard]] auto get_index() const -> unsigned long { return index; }
    };
 
    template <typename G, typename V>
    concept SequenceGenerator = requires(G g, unsigned long n) {
        { g.pop() } -> std::convertible_to<V>;
        { g.peek() } -> std::convertible_to<V>;
        g.skip(static_cast<unsigned int>(n));
        g.reseed(n);
        { g.get_index() } -> std::convertible_to<unsigned long>;
    };
 
    template <typename Derived, typename Value> class GeneratorBase {
      public:
        auto pop() -> Value {
            auto count_value
                = this->count_.fetch_add(1, std::memory_order_relaxed) + 1;  // ignore 0
            return derived().value_at(count_value);
        }
 
        [[nodiscard]] auto peek() -> Value {
            auto count_value = this->count_.load(std::memory_order_relaxed) + 1;
            return derived().value_at(count_value);
        }
 
        auto skip(unsigned int n) -> void { this->count_.fetch_add(n, std::memory_order_relaxed); }
 
        auto reseed(const unsigned long& seed) -> void {
            this->count_.store(seed, std::memory_order_relaxed);
        }
 
        [[nodiscard]] auto get_index() const -> unsigned long {
            return this->count_.load(std::memory_order_relaxed);
        }
 
      protected:
        std::atomic<unsigned long> count_{0};  
 
      private:
        friend Derived;  // for classes deriving directly from GeneratorBase (ilds.hpp, lds_n.hpp)
        template <typename, typename> friend class GeneratorIterable;
        GeneratorBase() = default;
 
        auto derived() -> Derived& { return static_cast<Derived&>(*this); }
    };
 
    template <typename Derived, typename Value> class GeneratorIterable
        : public GeneratorBase<Derived, Value> {
 
      private:
        friend Derived;
        GeneratorIterable() = default;
 
      public:
        auto begin() -> GeneratorIterator<Derived, Value> {
            return GeneratorIterator<Derived, Value>(static_cast<Derived*>(this));
        }
 
        [[nodiscard]] auto end() const -> GeneratorIterator<Derived, Value> {
            return GeneratorIterator<Derived, Value>(nullptr,
                                                     std::numeric_limits<unsigned long>::max());
        }
    };
 
    constexpr auto vdc(unsigned long count, unsigned long base) -> double {
        auto reslt = 0.0;
        auto denom = 1.0;
        while (count != 0) {
            const auto remainder = count % base;
            count /= base;
            denom *= static_cast<double>(base);
            reslt += static_cast<double>(remainder) / denom;
        }
        return reslt;
    }
 
    class VdCorput : public GeneratorIterable<VdCorput, double> {
        unsigned long base;
        std::array<double, MAX_REVERSE_BITS> rev_lst{};
        static_assert(MAX_REVERSE_BITS >= sizeof(unsigned long) * 8,
                      "MAX_REVERSE_BITS must be at least the number of bits in unsigned long");
 
      public:
        explicit VdCorput(const unsigned long base) : base{base} {
            double reverse = 1.0;
            for (unsigned int i = 0; i < MAX_REVERSE_BITS; ++i) {
                reverse /= static_cast<double>(base);
                this->rev_lst[i] = reverse;
            }
        }
 
        auto value_at(unsigned long n) const -> double {
            return detail::vdc_digit_sum<double>(n, this->base, this->rev_lst);
        }
    };
 
    class Halton : public GeneratorIterable<Halton, std::array<double, 2>> {
        VdCorput vdc0;
        VdCorput vdc1;
 
      public:
        Halton(const unsigned long base0, const unsigned long base1) : vdc0(base0), vdc1(base1) {}
 
        auto value_at(unsigned long n) const -> std::array<double, 2> {
            return {this->vdc0.value_at(n), this->vdc1.value_at(n)};
        }
    };
 
    class Circle : public GeneratorIterable<Circle, std::array<double, 2>> {
        VdCorput vdc;
 
      public:
        explicit Circle(const unsigned long base) : vdc(base) {}
 
        auto value_at(unsigned long n) const -> std::array<double, 2> {
            auto theta = this->vdc.value_at(n) * TWO_PI;  // map to [0, 2*pi];
            return {std::cos(theta), std::sin(theta)};
        }
    };
 
    class Disk : public GeneratorIterable<Disk, std::array<double, 2>> {
        VdCorput vdc0;
        VdCorput vdc1;
 
      public:
        Disk(const unsigned long base0, const unsigned long base1) : vdc0(base0), vdc1(base1) {}
 
        auto value_at(unsigned long n) const -> std::array<double, 2> {
            auto theta = this->vdc0.value_at(n) * TWO_PI;  // map to [0, 2*pi];
            auto radius = std::sqrt(this->vdc1.value_at(n));
            return {radius * std::cos(theta), radius * std::sin(theta)};
        }
    };
 
    class Sphere : public GeneratorIterable<Sphere, std::array<double, 3>> {
        VdCorput vdcgen;
        Circle cirgen;
 
      public:
        Sphere(const unsigned long base0, const unsigned long base1)
            : vdcgen(base0), cirgen(base1) {}
 
        auto value_at(unsigned long n) const -> std::array<double, 3> {
            auto cosphi = (MAPPING_FACTOR * this->vdcgen.value_at(n)) - 1.0;  // map to [-1, 1];
            auto sinphi = std::sqrt(1.0 - (cosphi * cosphi));
            auto arr = this->cirgen.value_at(n);
            return {sinphi * arr[0], sinphi * arr[1], cosphi};
        }
    };
 
    class Sphere3Hopf : public GeneratorIterable<Sphere3Hopf, std::array<double, 4>> {
        VdCorput vdc0;
        VdCorput vdc1;
        VdCorput vdc2;
 
      public:
        Sphere3Hopf(const unsigned long base0, const unsigned long base1, const unsigned long base2)
            : vdc0(base0), vdc1(base1), vdc2(base2) {}
 
        auto value_at(unsigned long n) const -> std::array<double, 4> {
            auto phi = this->vdc0.value_at(n) * TWO_PI;  // map to [0, 2*pi];
            auto psy = this->vdc1.value_at(n) * TWO_PI;  // map to [0, 2*pi];
            auto vdc = this->vdc2.value_at(n);
            auto cos_eta = std::sqrt(vdc);
            auto sin_eta = std::sqrt(1.0 - vdc);
            return {
                cos_eta * std::cos(psy),
                cos_eta * std::sin(psy),
                sin_eta * std::cos(phi + psy),
                sin_eta * std::sin(phi + psy),
            };
        }
    };
 
    // Compile-time contract checks: every generator satisfies the protocol concept.
    static_assert(SequenceGenerator<VdCorput, double>);
    static_assert(SequenceGenerator<Halton, std::array<double, 2>>);
    static_assert(SequenceGenerator<Circle, std::array<double, 2>>);
    static_assert(SequenceGenerator<Disk, std::array<double, 2>>);
    static_assert(SequenceGenerator<Sphere, std::array<double, 3>>);
    static_assert(SequenceGenerator<Sphere3Hopf, std::array<double, 4>>);
 
    extern unsigned long dummy(unsigned int index);
}  // namespace ldsgen