1// Deterministic, stream-addressable RNG (xoshiro256++ seeded by splitmix64 from (seed, stream)).2// Each (seed, stream) pair gives an independent, reproducible stream suitable for parallel Monte Carlo.3#pragma once4#include <cstdint>5#include <cmath>6#include <array>7#include "raid/core/types.hpp"89namespace raid {1011class Rng {12 public:13 explicit Rng(u64 seed = 0x5eed'0000'0000'0001ULL, u64 stream = 0) { reseed(seed, stream); }14 void reseed(u64 seed, u64 stream) {15 seed_ = seed; stream_ = stream;16 u64 z = seed ^ (stream * 0x9E3779B97F4A7C15ULL) ^ 0xD1B54A32D192ED03ULL;17 for (auto& w : s_) w = splitmix(z);18 if (s_[0] == 0 && s_[1] == 0 && s_[2] == 0 && s_[3] == 0) s_[0] = 1;19 }20 u64 seed() const { return seed_; }21 u64 stream() const { return stream_; }22 u64 next_u64() {23 const u64 result = rotl(s_[0] + s_[3], 23) + s_[0];24 const u64 t = s_[1] << 17;25 s_[2] ^= s_[0]; s_[3] ^= s_[1]; s_[1] ^= s_[2]; s_[0] ^= s_[3]; s_[2] ^= t; s_[3] = rotl(s_[3], 45);26 return result;27 }28 /// Uniform on [0,1) with 53 random bits.29 f64 uniform() { return static_cast<f64>(next_u64() >> 11) * 0x1.0p-53; }30 /// Standard normal, polar Box–Muller (caches the second draw).31 f64 normal() {32 if (has_cached_) { has_cached_ = false; return cached_; }33 f64 u, v, s;34 do { u = 2.0 * uniform() - 1.0; v = 2.0 * uniform() - 1.0; s = u * u + v * v; } while (s >= 1.0 || s == 0.0);35 const f64 m = std::sqrt(-2.0 * std::log(s) / s);36 cached_ = v * m; has_cached_ = true;37 return u * m;38 }39 /// Gamma(shape k, scale 1), Marsaglia–Tsang.40 f64 gamma(f64 k) {41 if (k < 1.0) { const f64 u = uniform(); return gamma(k + 1.0) * std::pow(u, 1.0 / k); }42 const f64 d = k - 1.0 / 3.0, c = 1.0 / std::sqrt(9.0 * d);43 for (;;) {44 f64 x, v;45 do { x = normal(); v = 1.0 + c * x; } while (v <= 0.0);46 v = v * v * v;47 const f64 u = uniform();48 if (u < 1.0 - 0.0331 * x * x * x * x) return d * v;49 if (std::log(u) < 0.5 * x * x + d * (1.0 - v + std::log(v))) return d * v;50 }51 }52 /// Student-t with nu degrees of freedom, unit scale (variance nu/(nu-2)).53 f64 student_t(f64 nu) { return normal() / std::sqrt(gamma(nu / 2.0) * 2.0 / nu); }54 /// Integer in [0, n)55 u64 below(u64 n) { return next_u64() % n; }56 /// Geometric number of trials with success probability p (support 1,2,…) — for the stationary bootstrap.57 u64 geometric(f64 p) { return 1 + static_cast<u64>(std::floor(std::log(1.0 - uniform()) / std::log(1.0 - p))); }5859 private:60 static constexpr u64 rotl(u64 x, int k) { return (x << k) | (x >> (64 - k)); }61 static u64 splitmix(u64& z) {62 z += 0x9E3779B97F4A7C15ULL;63 u64 r = z; r = (r ^ (r >> 30)) * 0xBF58476D1CE4E5B9ULL; r = (r ^ (r >> 27)) * 0x94D049BB133111EBULL;64 return r ^ (r >> 31);65 }66 std::array<u64, 4> s_{};67 u64 seed_{0}, stream_{0};68 f64 cached_{0}; bool has_cached_{false};69};7071} // namespace raid72