Rust implementations of Xoshiro256**, PCG32, and a bonus SplitMix64

rand.rs

//! This is a direct port of these algorithms' respective
//! C implementations into Rust. I did this chiefly as a learning
//! exercise for myself, but after adopting the `wrapping_mul` and
//! `rotate_left` methods (which I discovered from looking at the
//! code in the `rand` crate) I ended up with essentially identical
//! implementations to what's in the popular `rand` crate. Also here
//! is a very quick-and-dirty main() that demonstrates some of the
//! functionality.
//! 
//! See https://payloadgame.dev/devlog/2019/02/28/random-generators.html
//! where I discuss the relative merits and benchmark results.

use std;

fn main() {
    let seed = 1880u64;
    let seq = 1880u64;
    
    let mut pcg = Pcg32::new(seed, seq);
    
    println!("PCG: {:?}", pcg);
    //for _x in 0..100 {
    //    let n = pcg.rand();
    //    println!("Random: {:?}", n);
    //}
    
    // Let's roll some dice!!
    let mut results = [0,0,0,0,0,0,0,0,0,0];
    for _x in 0..1_000_000 {
        let n = pcg.nroll(2, 4);
        results[n as usize] += 1;
    }
    println!("{:?}", results);
    
    let mut xo = Xoshiro256::new(seed);
    println!("{:?}", xo);
    for _x in 0..100 {
        let n = xo.rand();
        println!("Random: {:?}", n);
    }
}

/// Implementation in Rust of xoshiro256**
/// http://xoshiro.di.unimi.it/xoshiro256starstar.c
#[derive(Debug)]
struct Xoshiro256 ( u64, u64, u64, u64 );

impl Xoshiro256 {
    fn new(seed: u64) -> Xoshiro256 {
        //Per Blackman & Vigna, seed with output of splitmix64
        let mut seeder = SplitMix64::new(seed);
        Xoshiro256 ( seeder.next(), seeder.next(), seeder.next(), seeder.next() )
    }
    
    fn rand(&mut self) -> u64 {
        let result = self.1.wrapping_mul(5).rotate_left(7).wrapping_mul(9);
        
        let t = self.1 << 17;
        self.2 ^= self.0;
    	self.3 ^= self.1;
    	self.1 ^= self.2;
    	self.0 ^= self.3;
    
    	self.2 ^= t;
    
    	self.3 = self.3.rotate_left(45);
    	
    	result as u64
    }
}

/// http://xoshiro.di.unimi.it/splitmix64.c
#[derive(Debug)]
struct SplitMix64(u64);

impl SplitMix64 {
    fn new(seed: u64) -> SplitMix64 {
        SplitMix64(seed)
    }
    
    fn next(&mut self) -> u64 {
        // We need to do all our math as 128bits to avoid overflow
        // but we then need to drop all of the top 64 bits for correct
        // behavior
        let mut z = self.0.wrapping_add(0x9e3779b97f4a7c15);
        self.0 = z;
        
        z = (z ^ (z >> 30)).wrapping_mul(0xbf58476d1ce4e5b9);
        z = (z ^ (z >> 27)).wrapping_mul(0x94d049bb133111eb);
        z ^ (z >> 31)
    }
}

#[derive(Debug)]
struct Pcg32 {
    state: u64,
    stream: u64,
}

impl Pcg32 {
    fn new(seed: u64, stream: u64) -> Pcg32 {
        let mut pcg = Pcg32 { state: 0, stream: (stream << 1u8) | 1u64 };
        let mut _r = pcg.rand();
        pcg.state += seed;
        _r = pcg.rand();
        
        pcg
    }
    
    fn rand(&mut self) -> u32 {
        let oldstate = self.state;
        self.state = oldstate.wrapping_mul(6364136223846793005).wrapping_add(self.stream);
        
        // `oldstate` is u64, but after these shifts and xor we only want a u32
        let xorshifted = (((oldstate >> 18) ^ oldstate) >> 27) as u32;
        // By virtue of the 59-bit shift, rot is at most 5 bits; casting it to signed for the negation later
        let rot = (oldstate >> 59) as i32;
        
        (xorshifted >> rot) | (xorshifted << ( -rot & 31 ) )
    }
    
    fn rand_range(&mut self, min: u32, max: u32) -> u32 {
        let range = max - min;
        
        self.rand_bound(range) + min
    }
    
    fn rand_bound(&mut self, bound: u32) -> u32 {
        // This is magic math that gives us how many numbers
        // we have to throw out (from our overall range) in
        // order to produce an *unbiased* result in the
        // desired range:
        // (MAX - threshold) % bound == 0
        let threshold = (std::u32::MAX - bound) % bound;
        
        // Uniformity guarantees that this loop will eventually
        // terminate, and in fact in practice we'll be given
        // small bounds which will mean only a very small slice
        // of our total range will be being ignored; in other
        // words, most of the time this loop will exit on the
        // first iteration.
        let mut n: u32;
        loop {
            n = self.rand();
            if n > threshold {
                break;
            }
        }
        
        n % bound
    }
    
    fn roll(&mut self, sides: u32) -> u32 {
        self.rand_bound(sides) + 1
    }
    
    fn nroll(&mut self, dice: u32, sides: u32) -> u32 {
        let mut sum = 0;
        for _i in 0..dice {
            sum += self.roll(sides);
        }
        
        sum
    }
}