vj-sim 0.3.0

//! The solver's randomness source (plan section 1.6).
//!
//! The solver never seeds anything: every draw advances a caller-supplied
//! `u64` state, which for `randomize()` is the object's own
//! `ObjEntry.prng_seed` (LRM 18.14: random stability is per object). Two
//! solves from equal states therefore produce equal results, and the draw
//! function is the *same* splitmix64 step `$urandom` and the process RNG use
//! (`sched::process_rng::splitmix64_next`), not a second generator.

use crate::sched::process_rng::splitmix64_next;
use num_bigint::BigUint;
use num_traits::Zero;

/// One 64-bit draw; advances `seed`.
pub fn draw(seed: &mut u64) -> u64 {
    splitmix64_next(seed)
}

/// A uniformly distributed value in `[0, bound)`; `bound` must be non-zero.
///
/// Rejection sampling over whole 64-bit draws with the top word masked to
/// `bound`'s bit length, so the result is exactly uniform (no modulo bias)
/// and at most one in two tries is rejected.
///
/// # Panics
/// If `bound` is zero (there is no value below zero; callers check first).
pub fn below(seed: &mut u64, bound: &BigUint) -> BigUint {
    assert!(!bound.is_zero(), "rng::below: empty range");
    let bits = bound.bits();
    let words = bits.div_ceil(64) as usize;
    let top_bits = bits - 64 * (words as u64 - 1);
    let top_mask = if top_bits == 64 {
        u64::MAX
    } else {
        (1u64 << top_bits) - 1
    };
    loop {
        let mut digits = Vec::with_capacity(words);
        for i in 0..words {
            let d = draw(seed);
            digits.push(if i == words - 1 { d & top_mask } else { d });
        }
        let v = BigUint::from_slice(
            &digits
                .iter()
                .flat_map(|d| [*d as u32, (*d >> 32) as u32])
                .collect::<Vec<u32>>(),
        );
        if &v < bound {
            return v;
        }
    }
}

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

    #[test]
    fn draw_is_the_process_rng_step() {
        let (mut a, mut b) = (7u64, 7u64);
        assert_eq!(draw(&mut a), splitmix64_next(&mut b));
        assert_eq!(a, b);
    }

    #[test]
    fn below_stays_in_range_and_is_deterministic() {
        for bound in [1u64, 2, 3, 5, 1000, (1 << 63) + 1] {
            let b = BigUint::from(bound);
            let (mut s1, mut s2) = (42u64, 42u64);
            for _ in 0..200 {
                let v = below(&mut s1, &b);
                assert!(v < b, "{v} must be below {bound}");
                assert_eq!(v, below(&mut s2, &b));
            }
        }
        // Wider than one word.
        let wide = (BigUint::from(1u8) << 130) + BigUint::from(12345u32);
        let mut s = 1u64;
        let mut high_seen = false;
        for _ in 0..200 {
            let v = below(&mut s, &wide);
            assert!(v < wide);
            high_seen |= v.bits() > 129;
        }
        assert!(high_seen, "values must reach the top of a wide range");
    }

    #[test]
    fn below_is_unbiased_on_a_small_range() {
        let b = BigUint::from(3u8);
        let mut s = 9u64;
        let mut counts = [0u32; 3];
        for _ in 0..30_000 {
            let v: u32 = u32::try_from(below(&mut s, &b)).expect("below 3");
            counts[v as usize] += 1;
        }
        for c in counts {
            assert!((9_500..10_500).contains(&c), "skewed: {counts:?}");
        }
    }
}