add lehmer64 rng

s2ga.hpp

@@ -1,4 +1,74 @@
+#include <array>
+#include <bit>
+#include <random>
+#include <regex>
+
+namespace s2ga
+{
+    template<class... T>
+    struct always_false: std::false_type
+    {
+    };
+
+    class lehmer64
+    {
+    public:
+        lehmer64(__uint128_t seed = 0)
+        {
+            this->seed(seed);
+        }
+
+        void seed(__uint128_t z)
+        {
+            if(z == 0)
+            {
+                z = 0xC0FFEE;
+            }
+
+            s = z;
+        }
+
+        uint64_t operator()() noexcept
+        {
+            s *= 0xda942042e4dd58b5;
+            return s >> 64;
+        }
+
+        static constexpr uint64_t min()
+        {
+            return 0;
+        }
+
+        static constexpr uint64_t max()
+        {
+            return UINT64_MAX;
+        }
+
+        template<typename T>
+        T random(T min, T max)
+        {
+            if constexpr(std::is_integral_v<T>)
+            {
+                return std::uniform_int_distribution<T>(min, max)(*this);
+            }
+            else if constexpr(std::is_floating_point_v<T>)
+            {
+                return std::uniform_real_distribution<T>(min, max)(*this);
+            }
+            else
+            {
+                static_assert(always_false<T>::value,
+                              "unsupported distribution type");
+            }
+        }
+
+    private:
+        __uint128_t s;
+    };
+
+    static_assert(std::uniform_random_bit_generator<lehmer64>);
+}; // namespace s2ga
+
 inline void foo()
 {
-
 }

tests/test.cpp

@@ -1,15 +1,172 @@
 #include "../s2ga.hpp"
+#include <catch2/benchmark/catch_benchmark.hpp>
 #include <catch2/catch_all.hpp>
 #include <catch2/catch_approx.hpp>
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/matchers/catch_matchers_range_equals.hpp>
 #include <cmath>
 #include <numbers>
+#include <random>
+#include <vector>
 
 using namespace Catch;
 
 constexpr double TOLERANCE = 1e-4;
 
+TEST_CASE("lehmer64 rng")
+{
+    s2ga::lehmer64 rng(0);
+
+    REQUIRE(rng() != 0);
+
+    const int N = 100;
+
+    {
+        std::vector<int> ivalues;
+        rng.seed(4);
+
+        for(int i = 0; i < N; ++i)
+        {
+            ivalues.emplace_back(rng.random(0, 256));
+        }
+
+        rng.seed(4);
+
+        for(int i = 0; i < N; ++i)
+        {
+            REQUIRE(ivalues[i] == rng.random(0, 256));
+        }
+    }
+
+    {
+        std::vector<float> fvalues;
+        rng.seed(5);
+
+        for(int i = 0; i < N; ++i)
+        {
+            fvalues.emplace_back(rng.random(-256.0f, 256.0f));
+        }
+
+        rng.seed(5);
+
+        for(int i = 0; i < N; ++i)
+        {
+            REQUIRE(fvalues[i] == rng.random(-256.0f, 256.0f));
+        }
+    }
+}
+
+TEST_CASE("Uniform Distribution Tests")
+{
+    s2ga::lehmer64 rng(42); // Fixed seed for reproducibility
+
+    SECTION("Integer Uniformity - Chi-Squared Test")
+    {
+        const int min = 0;
+        const int max = 9;
+        const int num_samples = 1'000'000;
+        const int num_bins = max - min + 1;
+        const double expected = num_samples / static_cast<double>(num_bins);
+        const double critical_value = 16.92; // χ²(0.05, 9)
+
+        std::vector<int> counts(num_bins, 0);
+        for(int i = 0; i < num_samples; ++i)
+        {
+            int val = rng.random(min, max);
+            ++counts[val - min];
+        }
+
+        double chi_sq = 0.0;
+        for(int count: counts)
+        {
+            double diff = count - expected;
+            chi_sq += (diff * diff) / expected;
+        }
+
+        CHECK(chi_sq < critical_value);
+    }
+
+    SECTION("Floating Point Uniformity - Kolmogorov-Smirnov Test")
+    {
+        const double min = 0.0;
+        const double max = 1.0;
+        const int num_samples = 100'000;
+        const double ks_critical = 1.36 / std::sqrt(num_samples); // α=0.05
+
+        std::vector<double> samples(num_samples);
+        std::generate(samples.begin(), samples.end(),
+                      [&] { return rng.random(min, max); });
+
+        std::sort(samples.begin(), samples.end());
+
+        double d_plus = 0.0;
+        double d_minus = 0.0;
+        for(size_t i = 0; i < samples.size(); ++i)
+        {
+            double fn = (i + 1.0) / num_samples;
+            double f = samples[i];
+            d_plus = std::max(d_plus, fn - f);
+            d_minus = std::max(d_minus, f - (i / (num_samples - 1.0)));
+        }
+
+        double d_stat = std::max(d_plus, d_minus);
+        CHECK(d_stat < ks_critical);
+    }
+
+    SECTION("Edge Case Coverage")
+    {
+        const int iterations = 10'000;
+
+        // Test minimum and maximum inclusion
+        bool hit_min = false;
+        bool hit_max = false;
+        for(int i = 0; i < iterations; ++i)
+        {
+            int val = rng.random(1, 10);
+            if(val == 1)
+                hit_min = true;
+            if(val == 10)
+                hit_max = true;
+        }
+        CHECK(hit_min);
+        CHECK(hit_max);
+
+        // Floating point bounds check
+        for(int i = 0; i < iterations; ++i)
+        {
+            double val = rng.random(0.0, 1.0);
+            CHECK(val >= 0.0);
+            CHECK(val < 1.0);
+        }
+    }
+}
+
+TEST_CASE("random benchmarking", "[benchmark]")
+{
+    s2ga::lehmer64 rng(0);
+
+    const int N = 100'000'000;
+
+    std::mt19937 std_rng(4);
+    std::uniform_int_distribution dist(-512, 512);
+
+    BENCHMARK("mt19937 100'000'000")
+    {
+        for(int i = 0; i < N; ++i)
+        {
+            (void)dist(std_rng);
+        }
+    };
+
+    BENCHMARK("lehmer64 100'000'000")
+    {
+        for(int i = 0; i < N; ++i)
+        {
+            (void)dist(rng);
+        }
+    };
+}
+
 inline double sphere(double x, double y)
 {
     return std::pow(x, 2.0) + std::pow(y, 2.0);