forge_core/

rng.rs

//! Deterministic RNG (SplitMix64), no external dependencies.
//!
//! Reproducibility is a design goal of the whole engine family: the same seed
//! must yield the same result, on any platform, with or without parallelism.
//! This is the exact generator already used by `anvil-core` and
//! `rainflow-core`, lifted here so every optimizer in the ecosystem shares one
//! certified source of randomness.
//!
//! For parallel runs (island models, ensemble restarts) derive an independent
//! stream per worker with [`Rng::split`]; distinct stream ids give
//! statistically independent, fully reproducible sub-sequences. (Formally the
//! derived states live on the same SplitMix64 orbit, so overlap is not
//! impossible — merely astronomically unlikely, ~`streams²·draws/2⁶⁴`.)

/// Derives a deterministic stream seed from a base `seed` and a stream `id`.
///
/// Mixes the id into the seed through the SplitMix64 finalizer so that even
/// adjacent ids (0, 1, 2, …) produce well-separated seeds. Used to give each
/// island/ensemble worker an independent, reproducible sub-stream.
pub fn mix_seed(seed: u64, id: u64) -> u64 {
    let mut z = seed ^ id.wrapping_mul(0x9E37_79B9_7F4A_7C15);
    z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
    z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
    z ^ (z >> 31)
}

/// A SplitMix64 generator. Fast, good statistical quality, single-word state.
#[derive(Clone, Debug)]
pub struct Rng {
    state: u64,
}

impl Rng {
    /// Creates an RNG from a seed.
    pub fn new(seed: u64) -> Self {
        Rng { state: seed }
    }

    /// Derives an independent generator for parallel stream `id`.
    ///
    /// The child seed is [`mix_seed`]`(seed, id)`, so different `(seed, id)`
    /// pairs give non-overlapping sub-streams deterministically — the basis for
    /// reproducible island/ensemble parallelism.
    pub fn split(seed: u64, id: u64) -> Self {
        Rng::new(mix_seed(seed, id))
    }

    #[inline]
    fn next_u64(&mut self) -> u64 {
        self.state = self.state.wrapping_add(0x9E37_79B9_7F4A_7C15);
        let mut z = self.state;
        z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
        z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
        z ^ (z >> 31)
    }

    /// Uniform `f64` in `[0, 1)` with 53 bits of mantissa.
    #[inline]
    pub fn uniform(&mut self) -> f64 {
        (self.next_u64() >> 11) as f64 / (1u64 << 53) as f64
    }

    /// Uniform `f64` in `[lo, hi]` (the top end is reachable: `lo + (hi−lo)·u`
    /// with `u` half an ulp below 1 can round up to `hi`, which is fine for
    /// the inclusive bounds used throughout the crate).
    #[inline]
    pub fn uniform_in(&mut self, lo: f64, hi: f64) -> f64 {
        lo + (hi - lo) * self.uniform()
    }

    /// Standard normal sample via the Box–Muller transform.
    #[inline]
    pub fn normal(&mut self) -> f64 {
        let u1 = (1.0 - self.uniform()).max(f64::MIN_POSITIVE); // avoid ln(0)
        let u2 = self.uniform();
        (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos()
    }

    /// `usize` uniform in `[0, bound)`.
    ///
    /// Uses a plain modulo reduction: the bias is ~`bound/2⁶⁴` (immeasurable
    /// for the population/dimension sizes used here) and, unlike rejection
    /// sampling, it consumes exactly one draw — which keeps RNG consumption,
    /// and therefore byte-level reproducibility, independent of `bound`.
    ///
    /// # Panics
    /// If `bound == 0`.
    #[inline]
    pub fn index(&mut self, bound: usize) -> usize {
        assert!(bound > 0, "Rng::index: bound must be > 0");
        (self.next_u64() % bound as u64) as usize
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn deterministic() {
        let mut a = Rng::new(123);
        let mut b = Rng::new(123);
        for _ in 0..1000 {
            assert_eq!(a.next_u64(), b.next_u64());
        }
    }

    #[test]
    fn uniform_in_range() {
        let mut r = Rng::new(7);
        for _ in 0..10_000 {
            let x = r.uniform();
            assert!((0.0..1.0).contains(&x));
        }
    }

    #[test]
    fn index_in_range() {
        let mut r = Rng::new(9);
        for _ in 0..10_000 {
            assert!(r.index(13) < 13);
        }
    }

    /// Split streams are reproducible and diverge from each other.
    #[test]
    fn split_streams_are_independent_and_reproducible() {
        let mut s0a = Rng::split(42, 0);
        let mut s0b = Rng::split(42, 0);
        let mut s1 = Rng::split(42, 1);
        // Same (seed, id) reproduces exactly.
        assert_eq!(s0a.next_u64(), s0b.next_u64());
        // Different ids diverge (overwhelmingly likely; fixed seeds make it sure).
        let a: Vec<u64> = (0..8).map(|_| s0a.next_u64()).collect();
        let b: Vec<u64> = (0..8).map(|_| s1.next_u64()).collect();
        assert_ne!(a, b);
    }

    /// The normal sampler is centered near zero with ~unit spread.
    #[test]
    fn normal_is_roughly_standard() {
        let mut r = Rng::new(2024);
        let n = 100_000;
        let mut mean = 0.0;
        let mut m2 = 0.0;
        for k in 1..=n {
            let x = r.normal();
            let d = x - mean;
            mean += d / k as f64;
            m2 += d * (x - mean);
        }
        let var = m2 / (n as f64 - 1.0);
        assert!(mean.abs() < 0.02, "mean {mean}");
        assert!((var - 1.0).abs() < 0.05, "var {var}");
    }
}