Source
stdlib/random/random.cpp
1
// Copyright (c) 2026 BigBrain LLC. MIT-licensed (see LICENSE).2
// Original work; see ACKNOWLEDGMENTS.md for the open-source ideas we build upon.3
#include "random.hpp"5
#include <cmath>6
#include <cstdint>7
#include <random>9
namespace cheatah::random {11
namespace {13
// WHY THE DISTRIBUTIONS ARE WRITTEN OUT HERE INSTEAD OF USING <random>'s14
//15
// `std::mt19937_64` is fully specified by the standard: the same seed gives the same 64-bit stream16
// on every implementation, forever. The DISTRIBUTIONS are not. `std::uniform_real_distribution`,17
// `uniform_int_distribution` and `normal_distribution` are all implementation-defined, so libstdc++18
// and libc++ turn one identical engine stream into DIFFERENT numbers.19
//20
// That quietly broke the promise `seed()` makes. A program seeded with 42 printed one pi estimate on21
// Linux and another on macOS — caught the first time the test suite ran on Apple Silicon, by a test22
// whose own comment claimed "seeded RNG makes the estimate fully reproducible". It was reproducible23
// per platform, which is not the same thing and is the more dangerous kind of almost-true.24
//25
// Defining the mappings here makes a seed mean one thing everywhere. All three are the standard26
// constructions, chosen so the arithmetic is exact and the consumption pattern is fixed:27
// each call takes a known number of engine draws, so two runs cannot drift apart.29
std::mt19937_64& engine() {30
// One engine PER THREAD: the `thread` module makes concurrent random() calls reachable, and a31
// single shared mt19937_64 would be a data race (torn state, lost advances). Each thread32
// self-seeds from std::random_device on first use; seed(s) seeds the CALLING thread only.33
thread_local std::mt19937_64 e{std::random_device{}()};34
return e;35
}37
/// A double in [0, 1) from ONE engine draw. The top 53 bits scaled by 2^-53: 53 is exactly the38
/// mantissa width, so every representable value is reachable and the multiply is exact.39
double canonical() {40
return static_cast<double>(engine()() >> 11) * 0x1.0p-53;41
}43
} // namespace45
void seed(unsigned long long s) { engine().seed(s); }47
double random() { return canonical(); }49
double uniform(double a, double b) { return a + (b - a) * canonical(); }51
long long randint(long long a, long long b) {52
if (b < a) {53
return a;54
}55
// Unbiased over the inclusive range by REJECTION rather than modulo: taking `draw % span`56
// directly would favour the low end whenever span does not divide 2^64. Computed in unsigned57
// arithmetic so a range spanning the whole of long long (where b - a overflows a signed type)58
// is still handled exactly.59
const std::uint64_t span = static_cast<std::uint64_t>(b) - static_cast<std::uint64_t>(a) + 1u;60
if (span == 0u) { // the full 64-bit range: every draw is already uniform61
return static_cast<long long>(engine()());62
}63
const std::uint64_t limit = UINT64_MAX - (UINT64_MAX % span) - 1u;64
std::uint64_t draw = engine()();65
while (draw > limit) { // discard the biased tail and redraw66
draw = engine()();67
}68
return static_cast<long long>(static_cast<std::uint64_t>(a) + (draw % span));69
}71
double gauss(double mu, double sigma) {72
// Box-Muller, using ONE of the two values it produces. Keeping the spare would be cheaper but73
// would make a call's output depend on how many gauss() calls preceded it, so the stream would74
// no longer be a pure function of the seed and the call sequence. Two draws per call, always.75
double u1 = canonical();76
while (u1 <= 0.0) { // log(0) is undefined; a zero draw is astronomically rare but not impossible77
u1 = canonical();78
}79
const double u2 = canonical();80
constexpr double kTwoPi = 6.283185307179586476925286766559;81
return mu + sigma * std::sqrt(-2.0 * std::log(u1)) * std::cos(kTwoPi * u2);82
}84
} // namespace cheatah::random