From 72da4beff0656b1762cc000280dc03d5d4219e4e Mon Sep 17 00:00:00 2001
From: NukeBird <nukebird.dev@gmail.com>
Date: Sun, 9 Mar 2025 23:28:43 +0300
Subject: [PATCH] add lehmer64 rng
---
s2ga.hpp | 72 ++++++++++++++++++++++-
tests/test.cpp | 157 +++++++++++++++++++++++++++++++++++++++++++++++++
2 files changed, 228 insertions(+), 1 deletion(-)
diff --git a/s2ga.hpp b/s2ga.hpp
index f8bed49..3adad1b 100644
--- a/s2ga.hpp
+++ b/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()
{
-
}
diff --git a/tests/test.cpp b/tests/test.cpp
index 8810825..7bd637c 100644
--- a/tests/test.cpp
+++ b/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);