oorandom 9.3.3

//! A tiny, robust PRNG implementation.
//!
//! More specifically, it implements a single GOOD PRNG algorithm,
//! which is currently a permuted congruential generator.  It has two
//! implementations, one that returns `u32` and one that returns
//! `u64`.  It also has functions that return floats or integer
//! ranges.  And that's it.  What more do you need?
//!
//! For more info on PCG generators, see http://www.pcg-random.org/
//!
//! This was designed as a minimalist utility for video games.  No
//! promises are made about its quality, and if you use it for
//! cryptography you will get what you deserve.
//!
//! Works with `#![no_std]`, has no global state, and is generally
//! neato.

#![no_std]
use core::ops::Range;

/// A PRNG producing a 32-bit output.
#[derive(Copy, Clone, Debug, PartialEq)]
pub struct Rand32 {
    state: u64,
    inc: u64,
}

impl Rand32 {
    /// The default value for `increment`.
    /// This is basically arbitrary, it comes from the
    /// PCG reference C implementation:
    /// https://github.com/imneme/pcg-c/blob/master/include/pcg_variants.h#L284
    pub const DEFAULT_INC: u64 = 1442695040888963407;

    /// This is the number that you have to Really Get Right.
    ///
    /// The value used here is from the PCG C implementation:
    /// https://github.com/imneme/pcg-c/blob/master/include/pcg_variants.h#L278
    pub(crate) const MULTIPLIER: u64 = 6364136223846793005;

    /// Creates a new PRNG with the given seed and a default increment.
    pub fn new(seed: u64) -> Self {
        Self::new_inc(seed, Self::DEFAULT_INC)
    }

    /// Creates a new PRNG.  The two inputs, `seed` and `increment`,
    /// determine what you get; `increment` basically selects which
    /// sequence of all those possible the PRNG will produce, and the
    /// `seed` selects where in that sequence you start.
    ///
    /// Both are arbitrary; increment must be an odd number but this
    /// handles that for you
    pub fn new_inc(seed: u64, increment: u64) -> Self {
        let mut rng = Self {
            state: 0,
            inc: (increment << 1) | 1,
        };
        // TODO This initialization song-and-dance is a little odd; verify.
        let _ = rng.rand();
        rng.state = rng.state.wrapping_add(seed);
        let _ = rng.rand();
        rng
    }

    /// Produces a random `u32` in the range `[0, u32::MAX]`.
    pub fn rand(&mut self) -> u32 {
        let oldstate: u64 = self.state;
        self.state = oldstate
            .wrapping_mul(Self::MULTIPLIER)
            .wrapping_add(self.inc);
        let xorshifted: u32 = (((oldstate >> 18) ^ oldstate) >> 27) as u32;
        let rot: u32 = (oldstate >> 59) as u32;
        xorshifted.rotate_right(rot)
    }

    /// Produces a random `f32` in the range `[0.0, 1.0)`.
    pub fn rand_float(&mut self) -> f32 {
        const TOTAL_BITS: u32 = 32;
        const MANTISSA_SCALE: f32 = 1.0 / ((1u32 << core::f32::MANTISSA_DIGITS) as f32);
        let mut u = self.rand();
        u >>= TOTAL_BITS - core::f32::MANTISSA_DIGITS;
        u as f32 * MANTISSA_SCALE
    }

    /// Produces a random within the given bounds.  Like any `Range`,
    /// it includes the lower bound and excludes the upper one.
    ///
    /// This should be faster than `Self::rand() % end + start`, but the
    /// real advantage is it's more convenient.  Requires that
    /// `range.end <= range.start`.
    pub fn rand_range(&mut self, range: Range<u32>) -> u32 {
        // This is harder to do well than it looks, it seems.  I don't
        // trust Lokathor's implementation 'cause I don't understand
        // it, so I went to numpy's, which points me to "Lemire's
        // rejection algorithm": http://arxiv.org/abs/1805.10941
        //
        // Algorithms 3, 4 and 5 in that paper all seem fine modulo
        // minor performance differences, so this is algorithm 5.
        // It uses numpy's implementation, `buffered_bounded_lemire_uint32()`

        // TODO: make shifts, and's, math, etc non-panicking.

        debug_assert!(range.start < range.end);
        let range_starting_from_zero = 0..(range.end - range.start);

        let s: u32 = range_starting_from_zero.end;
        let mut m: u64 = u64::from(self.rand()) * u64::from(s);
        let mut leftover: u32 = (m & 0xFFFFFFFF) as u32;

        if leftover < s {
            // TODO: verify the wrapping_neg() here
            let threshold: u32 = s.wrapping_neg() % s;
            while leftover < threshold {
                m = u64::from(self.rand()) * u64::from(s);
                leftover = (m & 0xFFFFFFFF) as u32;
            }
        }
        (m >> 32) as u32 + range.start
    }
}

/// A PRNG producing a 64-bit output.
#[derive(Copy, Clone, Debug, PartialEq)]
pub struct Rand64 {
    state: u128,
    inc: u128,
}

impl Rand64 {
    /// The default value for `increment`.
    ///
    /// The value used here is from the PCG default C implementation: http://www.pcg-random.org/download.html
    pub const DEFAULT_INC: u128 = 0x2FE0E169_FFBD06E3_5BC307BD_4D2F814F;

    /// This is the number that you have to Really Get Right.
    ///
    /// The value used here is from the PCG C implementation:
    /// https://github.com/imneme/pcg-c/blob/master/include/pcg_variants.h#L288
    pub(crate) const MULTIPLIER: u128 = 47026247687942121848144207491837523525;

    /// Creates a new PRNG with the given seed and a default increment.
    pub fn new(seed: u128) -> Self {
        Self::new_inc(seed, Self::DEFAULT_INC)
    }

    /// Same as `Rand32::new_inc()`
    fn new_inc(seed: u128, increment: u128) -> Self {
        let mut rng = Self {
            state: 0,
            inc: (increment << 1) | 1,
        };
        // TODO: This initialization song-and-dance is a little odd; verify.
        let _ = rng.rand();
        rng.state = rng.state.wrapping_add(seed);
        let _ = rng.rand();
        rng
    }

    /// Produces a random `u64` in the range`[0, u64::MAX]`.
    pub fn rand(&mut self) -> u64 {
        let oldstate: u128 = self.state;
        self.state = oldstate
            .wrapping_mul(Self::MULTIPLIER)
            .wrapping_add(self.inc);
        let xorshifted: u64 = (((oldstate >> 29) ^ oldstate) >> 58) as u64;
        let rot: u32 = (oldstate >> 122) as u32;
        xorshifted.rotate_right(rot)
    }

    /// Produces a random `f64` in the range `[0.0, 1.0)`.
    pub fn rand_float(&mut self) -> f64 {
        const TOTAL_BITS: u32 = 64;
        const MANTISSA_SCALE: f64 = 1.0 / ((1u64 << core::f64::MANTISSA_DIGITS) as f64);
        let mut u = self.rand();
        u >>= TOTAL_BITS - core::f64::MANTISSA_DIGITS;
        u as f64 * MANTISSA_SCALE
    }

    /// Produces a random within the given bounds.  Like any `Range`,
    /// it includes the lower bound and excludes the upper one.
    ///
    /// This should be faster than `Self::rand() % end + start`, but the
    /// real advantage is it's more convenient.  Requires that
    /// `range.end <= range.start`.
    pub fn rand_range(&mut self, range: Range<u64>) -> u64 {
        // Same as `Rand32::rand_range()`
        // TODO: make shifts, and's, math, etc non-panicking.

        debug_assert!(range.start < range.end);
        let range_starting_from_zero = 0..(range.end - range.start);

        let s: u64 = range_starting_from_zero.end;
        let mut m: u128 = u128::from(self.rand()) * u128::from(s);
        let mut leftover: u64 = (m & 0xFFFFFFFF_FFFFFFFF) as u64;

        if leftover < s {
            // TODO: Verify the wrapping_negate() here
            let threshold: u64 = s.wrapping_neg() % s;
            while leftover < threshold {
                m = u128::from(self.rand()) * u128::from(s);
                leftover = (m & 0xFFFFFFFF_FFFFFFFF) as u64;
            }
        }
        (m >> 64) as u64
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use randomize::{self, PCG32, PCG64};

    #[test]
    fn test_rand32() {
        // Generate some random numbers and validate them against
        // a known-good generator.
        {
            let seed = 54321;
            let mut r1 = Rand32::new(seed);
            let mut r2 = PCG32::seed(seed, Rand32::DEFAULT_INC);
            for _ in 0..1000 {
                assert_eq!(r1.rand(), r2.next_u32());
            }
        }

        {
            let seed = 3141592653;
            let inc = 0xDEADBEEF;
            let mut r1 = Rand32::new_inc(seed, inc);
            let mut r2 = PCG32::seed(seed, inc);
            for _ in 0..1000 {
                assert_eq!(r1.rand(), r2.next_u32());
            }
        }
    }

    #[test]
    fn test_rand64() {
        // Generate some random numbers and validate them against
        // a known-good generator.
        {
            let seed = 54321;
            let mut r1 = Rand64::new(seed);
            let mut r2 = PCG64::seed(seed, Rand64::DEFAULT_INC);
            for _ in 0..1000 {
                assert_eq!(r1.rand(), r2.next_u64());
            }
        }

        {
            let seed = 3141592653;
            let inc = 0xDEADBEEF;
            let mut r1 = Rand64::new_inc(seed, inc);
            let mut r2 = PCG64::seed(seed, inc);
            for _ in 0..1000 {
                assert_eq!(r1.rand(), r2.next_u64());
            }
        }
    }

    #[test]
    fn test_float32() {
        {
            let seed = 2718281828;
            let mut r1 = Rand32::new(seed);
            let mut r2 = PCG32::seed(seed, Rand32::DEFAULT_INC);
            for _ in 0..1000 {
                // First just make sure they both work with randomize's
                // f32 conversion function -- sanity checks.
                let i1 = r1.rand();
                let i2 = r2.next_u32();
                assert_eq!(i1, i2);
                let f1 = randomize::f32_half_open_right(i1);
                let f2 = randomize::f32_half_open_right(i2);
                // We can directly compare floats 'cause we do no math, it's
                // literally the same bitwise algorithm with the same inputs.
                assert_eq!(f1, f2);

                // Make sure result is in [0.0, 1.0)
                assert!(f1 >= 0.0);
                assert!(f1 < 1.0);
            }

            for _ in 0..1000 {
                // Now make sure our own float conversion function works.
                let f1 = r1.rand_float();
                let f2 = randomize::f32_half_open_right(r2.next_u32());
                assert_eq!(f1, f2);
                assert!(f1 >= 0.0);
                assert!(f1 < 1.0);
            }
        }
    }

    #[test]
    fn test_randrange() {
        // Make sure ranges are valid and in the given range
        let seed = 2342_4223;
        let mut r1 = Rand32::new(seed);
        for _ in 0..1000 {
            // Generate our bounds at random
            let a = r1.rand();
            let b = r1.rand();
            if a == b {
                continue;
            }
            let (low, high) = if a < b { (a, b) } else { (b, a) };

            // Get a number in that range
            let in_range = r1.rand_range(low..high);
            assert!(in_range >= low);
            assert!(in_range < high);
        }
    }
}