|
1 | 1 | #pragma once |
2 | 2 |
|
3 | | -#include <cmath> |
4 | 3 | #include <cstdint> |
5 | | -#include <cstring> |
6 | 4 | #include <random> |
7 | 5 | #include <type_traits> |
8 | 6 |
|
9 | 7 | // Portable random distributions that produce identical sequences across |
10 | 8 | // compilers and platforms (Windows/MSVC, Linux/GCC, macOS/Clang). |
11 | 9 | // |
12 | | -// std::uniform_int_distribution and std::normal_distribution are NOT portable: |
13 | | -// - MSVC STL, libstdc++, and libc++ each use different algorithm implementations |
14 | | -// - The same seed produces different sequences on different platforms |
15 | | -// - This makes regression tests and graph builds non-deterministic across CI runners |
16 | | -// |
17 | | -// Both distributions operate on std::mt19937 (which IS portable — same seed |
18 | | -// produces the same uniform sequence) and use deterministic algorithms: |
19 | | -// - DeterministicUniformIntDistribution: bounded rejection sampling |
20 | | -// - DeterministicNormalDistribution: Box-Muller transform with portable math |
21 | | -// |
22 | | -// Reference: Box, G.E.P. and Muller, M.E. (1958), "A Note on the Generation |
23 | | -// of Random Normal Deviates", Annals of Mathematical Statistics, 29(2), pp.610-611. |
| 10 | +// std::uniform_int_distribution is NOT portable across MSVC STL, libstdc++, and libc++. |
| 11 | +// DeterministicUniformIntDistribution operates on std::mt19937 using bounded |
| 12 | +// rejection sampling to guarantee identical integer sequences on all platforms. |
24 | 13 |
|
25 | 14 | namespace deglib::random { |
26 | 15 |
|
27 | | -namespace detail { |
28 | | - |
29 | | -// Portable float natural logarithm (log(x)) without platform libm dependencies. |
30 | | -// Operates strictly on 32-bit floats with 24-bit mantissa precision. |
31 | | -inline float portable_log(float x) { |
32 | | - uint32_t ix; |
33 | | - std::memcpy(&ix, &x, sizeof(float)); |
34 | | - int exp = static_cast<int>((ix >> 23) & 0xFF) - 127; |
35 | | - ix = (ix & 0x007FFFFF) | 0x3F800000; |
36 | | - float m; |
37 | | - std::memcpy(&m, &ix, sizeof(float)); |
38 | | - if (m > 1.41421356f) { |
39 | | - m *= 0.5f; |
40 | | - exp += 1; |
41 | | - } |
42 | | - float f = m - 1.0f; |
43 | | - float s = f / (2.0f + f); |
44 | | - float z = s * s; |
45 | | - float w = z * z; |
46 | | - float R = z * (0.6666666666666666f + w * (0.4f + w * (0.2857142857142857f + w * (0.2222222222222222f + w * 0.1818181818181818f)))); |
47 | | - return static_cast<float>(exp) * 0.6931471805599453f + s * (2.0f + R); |
48 | | -} |
49 | | - |
50 | | -// Portable float cosine (cos(x)) without platform libm dependencies. |
51 | | -// Input x is in [0, 2*PI). Uses symmetry reduction to [0, PI/2] and a 10th-degree polynomial. |
52 | | -inline float portable_cos(float x) { |
53 | | - constexpr float PI = 3.14159265358979323846f; |
54 | | - constexpr float HALF_PI = 1.57079632679489661923f; |
55 | | - constexpr float TWO_PI = 6.28318530717958647692f; |
56 | | - |
57 | | - // Symmetry reduction to [0, PI] |
58 | | - if (x > PI) { |
59 | | - x = TWO_PI - x; |
60 | | - } |
61 | | - |
62 | | - // Symmetry reduction to [0, PI/2] |
63 | | - bool sign = false; |
64 | | - if (x > HALF_PI) { |
65 | | - x = PI - x; |
66 | | - sign = true; |
67 | | - } |
68 | | - |
69 | | - // Polynomial for cos(x) on [0, PI/2] accurate to 1e-7 (full float precision) |
70 | | - float x2 = x * x; |
71 | | - float res = 1.0f - x2 * (0.5f - x2 * (0.041666666666666664f - x2 * (0.0013888888888888889f - x2 * (0.0000248015873015873f - x2 * 0.000000275573192239859f)))); |
72 | | - |
73 | | - return sign ? -res : res; |
74 | | -} |
75 | | - |
76 | | -} // namespace detail |
77 | | - |
78 | 16 | /// Portable uniform integer distribution using bounded rejection sampling. |
79 | 17 | /// Guarantees identical integer sequences across Windows, Linux, and macOS. |
80 | 18 | template <typename IntType = int> |
@@ -112,35 +50,4 @@ class DeterministicUniformIntDistribution { |
112 | 50 | uint64_t max_valid_{4294967295ULL}; |
113 | 51 | }; |
114 | 52 |
|
115 | | -/// Portable normal distribution using the Box-Muller transform with portable math. |
116 | | -/// Always consumes exactly 2 RNG values per call (no caching/spare state) |
117 | | -/// to ensure deterministic behavior regardless of how multiple instances |
118 | | -/// share the same RNG. |
119 | | -class DeterministicNormalDistribution { |
120 | | -public: |
121 | | - explicit DeterministicNormalDistribution(float mean, float stddev) |
122 | | - : mean_(mean), stddev_(stddev) {} |
123 | | - |
124 | | - // Generate a single normally-distributed float using the provided RNG. |
125 | | - // Operates strictly on 24-bit mantissa float division and portable math |
126 | | - // to guarantee 100% bit-identical sequences across Windows, Linux, and macOS. |
127 | | - float operator()(std::mt19937& rng) { |
128 | | - float u1, u2; |
129 | | - do { |
130 | | - uint32_t v1 = rng() >> 8; |
131 | | - u1 = static_cast<float>(v1) / 16777216.0f; |
132 | | - } while (u1 <= 0.0f); |
133 | | - |
134 | | - uint32_t v2 = rng() >> 8; |
135 | | - u2 = static_cast<float>(v2) / 16777216.0f; |
136 | | - |
137 | | - float mag = stddev_ * std::sqrt(-2.0f * detail::portable_log(u1)); |
138 | | - return mean_ + mag * detail::portable_cos(6.28318530717958647692f * u2); |
139 | | - } |
140 | | - |
141 | | -private: |
142 | | - float mean_; |
143 | | - float stddev_; |
144 | | -}; |
145 | | - |
146 | 53 | } // namespace deglib::random |
0 commit comments