diff --git a/s2ga.hpp b/s2ga.hpp index f8bed49..3adad1b 100644 --- a/s2ga.hpp +++ b/s2ga.hpp @@ -1,4 +1,74 @@ +#include +#include +#include +#include + +namespace s2ga +{ + template + 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 + T random(T min, T max) + { + if constexpr(std::is_integral_v) + { + return std::uniform_int_distribution(min, max)(*this); + } + else if constexpr(std::is_floating_point_v) + { + return std::uniform_real_distribution(min, max)(*this); + } + else + { + static_assert(always_false::value, + "unsupported distribution type"); + } + } + + private: + __uint128_t s; + }; + + static_assert(std::uniform_random_bit_generator); +}; // 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 #include #include #include #include #include #include +#include +#include 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 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 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(num_bins); + const double critical_value = 16.92; // χ²(0.05, 9) + + std::vector 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 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);