// Deterministic, stream-addressable RNG (xoshiro256++ seeded by splitmix64 from (seed, stream)).
// Each (seed, stream) pair gives an independent, reproducible stream suitable for parallel Monte Carlo.
#pragma once
#include <cstdint>
#include <cmath>
#include <array>
#include "raid/core/types.hpp"

namespace raid {

class Rng {
 public:
  explicit Rng(u64 seed = 0x5eed'0000'0000'0001ULL, u64 stream = 0) { reseed(seed, stream); }
  void reseed(u64 seed, u64 stream) {
    seed_ = seed; stream_ = stream;
    u64 z = seed ^ (stream * 0x9E3779B97F4A7C15ULL) ^ 0xD1B54A32D192ED03ULL;
    for (auto& w : s_) w = splitmix(z);
    if (s_[0] == 0 && s_[1] == 0 && s_[2] == 0 && s_[3] == 0) s_[0] = 1;
  }
  u64 seed() const { return seed_; }
  u64 stream() const { return stream_; }
  u64 next_u64() {
    const u64 result = rotl(s_[0] + s_[3], 23) + s_[0];
    const u64 t = s_[1] << 17;
    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);
    return result;
  }
  /// Uniform on [0,1) with 53 random bits.
  f64 uniform() { return static_cast<f64>(next_u64() >> 11) * 0x1.0p-53; }
  /// Standard normal, polar Box–Muller (caches the second draw).
  f64 normal() {
    if (has_cached_) { has_cached_ = false; return cached_; }
    f64 u, v, s;
    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);
    const f64 m = std::sqrt(-2.0 * std::log(s) / s);
    cached_ = v * m; has_cached_ = true;
    return u * m;
  }
  /// Gamma(shape k, scale 1), Marsaglia–Tsang.
  f64 gamma(f64 k) {
    if (k < 1.0) { const f64 u = uniform(); return gamma(k + 1.0) * std::pow(u, 1.0 / k); }
    const f64 d = k - 1.0 / 3.0, c = 1.0 / std::sqrt(9.0 * d);
    for (;;) {
      f64 x, v;
      do { x = normal(); v = 1.0 + c * x; } while (v <= 0.0);
      v = v * v * v;
      const f64 u = uniform();
      if (u < 1.0 - 0.0331 * x * x * x * x) return d * v;
      if (std::log(u) < 0.5 * x * x + d * (1.0 - v + std::log(v))) return d * v;
    }
  }
  /// Student-t with nu degrees of freedom, unit scale (variance nu/(nu-2)).
  f64 student_t(f64 nu) { return normal() / std::sqrt(gamma(nu / 2.0) * 2.0 / nu); }
  /// Integer in [0, n)
  u64 below(u64 n) { return next_u64() % n; }
  /// Geometric number of trials with success probability p (support 1,2,…) — for the stationary bootstrap.
  u64 geometric(f64 p) { return 1 + static_cast<u64>(std::floor(std::log(1.0 - uniform()) / std::log(1.0 - p))); }

 private:
  static constexpr u64 rotl(u64 x, int k) { return (x << k) | (x >> (64 - k)); }
  static u64 splitmix(u64& z) {
    z += 0x9E3779B97F4A7C15ULL;
    u64 r = z; r = (r ^ (r >> 30)) * 0xBF58476D1CE4E5B9ULL; r = (r ^ (r >> 27)) * 0x94D049BB133111EBULL;
    return r ^ (r >> 31);
  }
  std::array<u64, 4> s_{};
  u64 seed_{0}, stream_{0};
  f64 cached_{0}; bool has_cached_{false};
};

}  // namespace raid