Compare commits

..

2 commits

Author SHA1 Message Date
092e7c2102 add [test] tag 2025-03-09 23:31:09 +03:00
72da4beff0 add lehmer64 rng 2025-03-09 23:28:43 +03:00
2 changed files with 244 additions and 17 deletions

View file

@ -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() inline void foo()
{ {
} }

View file

@ -1,21 +1,178 @@
#include "../s2ga.hpp" #include "../s2ga.hpp"
#include <catch2/benchmark/catch_benchmark.hpp>
#include <catch2/catch_all.hpp> #include <catch2/catch_all.hpp>
#include <catch2/catch_approx.hpp> #include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp> #include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_range_equals.hpp> #include <catch2/matchers/catch_matchers_range_equals.hpp>
#include <cmath> #include <cmath>
#include <numbers> #include <numbers>
#include <random>
#include <vector>
using namespace Catch; using namespace Catch;
constexpr double TOLERANCE = 1e-4; constexpr double TOLERANCE = 1e-4;
TEST_CASE("lehmer64 rng", "[test]")
{
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", "[test]")
{
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) inline double sphere(double x, double y)
{ {
return std::pow(x, 2.0) + std::pow(y, 2.0); return std::pow(x, 2.0) + std::pow(y, 2.0);
} }
TEST_CASE("sphere(0, 0) = 0") TEST_CASE("sphere(0, 0) = 0", "[test]")
{ {
REQUIRE(sphere(0.0, 0.0) == Approx(0.0).margin(TOLERANCE)); REQUIRE(sphere(0.0, 0.0) == Approx(0.0).margin(TOLERANCE));
} }
@ -28,7 +185,7 @@ inline double ackley(double x, double y)
std::numbers::e + 20.0; std::numbers::e + 20.0;
} }
TEST_CASE("ackley(0, 0) = 0") TEST_CASE("ackley(0, 0) = 0", "[test]")
{ {
REQUIRE(ackley(0.0, 0.0) == Approx(0.0).margin(TOLERANCE)); REQUIRE(ackley(0.0, 0.0) == Approx(0.0).margin(TOLERANCE));
} }
@ -41,7 +198,7 @@ inline double rastrigin(double x, double y)
(std::pow(y, 2.0) - A * std::cos(2.0 * std::numbers::pi * y)); (std::pow(y, 2.0) - A * std::cos(2.0 * std::numbers::pi * y));
} }
TEST_CASE("rastrigin(0, 0) = 0") TEST_CASE("rastrigin(0, 0) = 0", "[test]")
{ {
REQUIRE(rastrigin(0.0, 0.0) == Approx(0.0).margin(TOLERANCE)); REQUIRE(rastrigin(0.0, 0.0) == Approx(0.0).margin(TOLERANCE));
} }
@ -51,7 +208,7 @@ inline double rosenbrock(double x, double y)
return 100.0 * std::pow(y - std::pow(x, 2.0), 2.0) + std::pow(1.0 - x, 2.0); return 100.0 * std::pow(y - std::pow(x, 2.0), 2.0) + std::pow(1.0 - x, 2.0);
} }
TEST_CASE("rosenbrock(1, 1) = 0") TEST_CASE("rosenbrock(1, 1) = 0", "[test]")
{ {
REQUIRE(rosenbrock(1.0, 1.0) == Approx(0.0).margin(TOLERANCE)); REQUIRE(rosenbrock(1.0, 1.0) == Approx(0.0).margin(TOLERANCE));
} }
@ -63,7 +220,7 @@ inline double bill(double x, double y)
std::pow(2.625 - x + x * std::pow(y, 3.0), 2.0); std::pow(2.625 - x + x * std::pow(y, 3.0), 2.0);
} }
TEST_CASE("bill(3, 0.5) = 0") TEST_CASE("bill(3, 0.5) = 0", "[test]")
{ {
REQUIRE(bill(3.0, 0.5) == Approx(0.0).margin(TOLERANCE)); REQUIRE(bill(3.0, 0.5) == Approx(0.0).margin(TOLERANCE));
} }
@ -78,7 +235,7 @@ inline double goldstein_price(double x, double y)
36.0 * x * y + 27.0 * std::pow(y, 2.0))); 36.0 * x * y + 27.0 * std::pow(y, 2.0)));
} }
TEST_CASE("goldstein_price(0, -1) = 3") TEST_CASE("goldstein_price(0, -1) = 3", "[test]")
{ {
REQUIRE(goldstein_price(0.0, -1.0) == Approx(3.0).margin(TOLERANCE)); REQUIRE(goldstein_price(0.0, -1.0) == Approx(3.0).margin(TOLERANCE));
} }
@ -88,7 +245,7 @@ inline double booth(double x, double y)
return std::pow(x + 2.0 * y - 7.0, 2.0) + std::pow(2.0 * x + y - 5.0, 2.0); return std::pow(x + 2.0 * y - 7.0, 2.0) + std::pow(2.0 * x + y - 5.0, 2.0);
} }
TEST_CASE("booth(1, 3) = 0") TEST_CASE("booth(1, 3) = 0", "[test]")
{ {
REQUIRE(booth(1.0, 3.0) == Approx(0.0).margin(TOLERANCE)); REQUIRE(booth(1.0, 3.0) == Approx(0.0).margin(TOLERANCE));
} }
@ -99,7 +256,7 @@ inline double bukin_n6(double x, double y)
0.01 * std::abs(x + 10.0); 0.01 * std::abs(x + 10.0);
} }
TEST_CASE("bukin_n6(-10, 1) = 0") TEST_CASE("bukin_n6(-10, 1) = 0", "[test]")
{ {
REQUIRE(bukin_n6(-10.0, 1.0) == Approx(0.0).margin(TOLERANCE)); REQUIRE(bukin_n6(-10.0, 1.0) == Approx(0.0).margin(TOLERANCE));
} }
@ -109,7 +266,7 @@ inline double matyas(double x, double y)
return 0.26 * sphere(x, y) - 0.48 * x * y; return 0.26 * sphere(x, y) - 0.48 * x * y;
} }
TEST_CASE("matyas(0, 0) = 0") TEST_CASE("matyas(0, 0) = 0", "[test]")
{ {
REQUIRE(matyas(0.0, 0.0) == Approx(0.0).margin(TOLERANCE)); REQUIRE(matyas(0.0, 0.0) == Approx(0.0).margin(TOLERANCE));
} }
@ -126,7 +283,7 @@ inline double levi_n13(double x, double y)
std::pow(y - 1.0, 2.0) * (1.0 + sin2(2.0 * std::numbers::pi * y)); std::pow(y - 1.0, 2.0) * (1.0 + sin2(2.0 * std::numbers::pi * y));
} }
TEST_CASE("levi_n13(1, 1) = 0") TEST_CASE("levi_n13(1, 1) = 0", "[test]")
{ {
REQUIRE(levi_n13(1.0, 1.0) == Approx(0.0).margin(TOLERANCE)); REQUIRE(levi_n13(1.0, 1.0) == Approx(0.0).margin(TOLERANCE));
} }
@ -137,7 +294,7 @@ inline double three_hump_camel(double x, double y)
std::pow(x, 6.0) / 6.0 + x * y + std::pow(y, 2.0); std::pow(x, 6.0) / 6.0 + x * y + std::pow(y, 2.0);
} }
TEST_CASE("three_hump_camel(0, 0) = 0") TEST_CASE("three_hump_camel(0, 0) = 0", "[test]")
{ {
REQUIRE(three_hump_camel(0.0, 0.0) == Approx(0.0).margin(TOLERANCE)); REQUIRE(three_hump_camel(0.0, 0.0) == Approx(0.0).margin(TOLERANCE));
} }
@ -148,7 +305,7 @@ inline double eggholder(double x, double y)
x * std::sin(std::sqrt(std::abs(x - (y + 47.0)))); x * std::sin(std::sqrt(std::abs(x - (y + 47.0))));
} }
TEST_CASE("eggholder(512, 404.2319) = -959.6407") TEST_CASE("eggholder(512, 404.2319) = -959.6407", "[test]")
{ {
REQUIRE(eggholder(512.0, 404.2319) == Approx(-959.6407).margin(TOLERANCE)); REQUIRE(eggholder(512.0, 404.2319) == Approx(-959.6407).margin(TOLERANCE));
} }
@ -158,7 +315,7 @@ inline double mccormick(double x, double y)
return std::sin(x + y) + std::pow(x - y, 2.0) - 1.5 * x + 2.5 * y + 1.0; return std::sin(x + y) + std::pow(x - y, 2.0) - 1.5 * x + 2.5 * y + 1.0;
} }
TEST_CASE("mccormick(-0.54719, -1.54719) = -1.9133") TEST_CASE("mccormick(-0.54719, -1.54719) = -1.9133", "[test]")
{ {
REQUIRE(mccormick(-0.54719, -1.54719) == Approx(-1.9133).margin(TOLERANCE)); REQUIRE(mccormick(-0.54719, -1.54719) == Approx(-1.9133).margin(TOLERANCE));
} }
@ -169,7 +326,7 @@ inline double schaffer_n2(double x, double y)
std::pow(1.0 + 0.001 * sphere(x, y), 2.0); std::pow(1.0 + 0.001 * sphere(x, y), 2.0);
} }
TEST_CASE("schaffer_n2(0, 0) = 0") TEST_CASE("schaffer_n2(0, 0) = 0", "[test]")
{ {
REQUIRE(schaffer_n2(0.0, 0.0) == Approx(0.0).margin(TOLERANCE)); REQUIRE(schaffer_n2(0.0, 0.0) == Approx(0.0).margin(TOLERANCE));
} }
@ -187,12 +344,12 @@ inline double schaffer_n4(double x, double y)
std::pow(1.0 + 0.001 * sphere(x, y), 2.0); std::pow(1.0 + 0.001 * sphere(x, y), 2.0);
} }
TEST_CASE("schaffer_n4(0, 1.25313) = 0.292579") TEST_CASE("schaffer_n4(0, 1.25313) = 0.292579", "[test]")
{ {
REQUIRE(schaffer_n4(0.0, 1.25313) == Approx(0.292579).margin(TOLERANCE)); REQUIRE(schaffer_n4(0.0, 1.25313) == Approx(0.292579).margin(TOLERANCE));
} }
TEST_CASE("Should pass") TEST_CASE("Should pass", "[test]")
{ {
REQUIRE_NOTHROW(foo()); REQUIRE_NOTHROW(foo());
} }