r2rs-nmath 0.1.1

use crate::traits::RNG;

static m1: u32 = 4294967087;
static m2: u32 = 4294944443;

#[derive(Clone, Debug)]
pub struct LEcuyerCMRG {
    state: [u32; 6],
}

impl Default for LEcuyerCMRG {
    fn default() -> Self {
        Self::new()
    }
}

impl LEcuyerCMRG {
    pub fn new() -> Self {
        let mut ret = Self { state: [0; 6] };
        ret.randomize();
        ret
    }

    pub fn fixup_seeds(&mut self) {
        let mut all_zero = (0..3).all(|j| {
            let temp = self.state[j];
            temp == 0
        });
        let mut all_ok = (0..3).all(|j| {
            let temp = self.state[j];
            temp < m1
        });
        if all_zero || !all_ok {
            self.randomize();
        }

        all_zero = (3..6).all(|j| {
            let temp = self.state[j];
            temp == 0
        });
        all_ok = (3..6).all(|j| {
            let temp = self.state[j];
            temp < m2
        });
        if all_zero || !all_ok {
            self.randomize();
        }
    }
}

impl RNG for LEcuyerCMRG {
    fn set_seed(&mut self, mut seed: u32) {
        for _ in 0..50 {
            seed = seed.wrapping_mul(69069).wrapping_add(1);
        }

        (0..6).for_each(|j| {
            seed = seed.wrapping_mul(69069).wrapping_add(1);
            while seed >= m2 {
                seed = seed.wrapping_mul(69069).wrapping_add(1);
            }
            self.state[j] = seed;
        });

        self.fixup_seeds()
    }

    fn get_state(&self) -> Vec<i32> {
        self.state.clone().iter().map(|&i| i as i32).collect()
    }

    fn unif_rand(&mut self) -> f64 {
        let mut p1 = 1403580 * self.state[1] as i64 - 810728 * self.state[0] as i64;
        /* p1 % m1 would surely do */
        let mut k = p1 / 4294967087;
        p1 -= k * 4294967087;
        if p1 < 0 {
            p1 += 4294967087
        }
        self.state[0] = self.state[1];
        self.state[1] = self.state[2];
        self.state[2] = p1 as u32;
        let mut p2 = 527612 * self.state[5] as i64 - 1370589 * self.state[3] as i64;
        k = p2 / 4294944443;
        p2 -= k as i64 * 4294944443;
        if p2 < 0 {
            p2 += 4294944443
        }
        self.state[3] = self.state[4];
        self.state[4] = self.state[5];
        self.state[5] = p2 as u32;
        return (if p1 > p2 {
            p1 - p2
        } else {
            (p1 - p2) + 4294967087
        }) as f64
            * 2.328306549295727688e-10f64;
    }
}