quantre/random/
rng_uniform.rs

1//! RNG from Uniform Distribution
2
3use std::cmp;
4use std::sync::atomic::AtomicU32;
5use std::sync::atomic::Ordering;
6use std::time::SystemTime;
7use std::time::UNIX_EPOCH;
8
9/// Global seed
10static SEED: AtomicU32 = AtomicU32::new(0);
11
12/// Pseudorandom seed using the current time
13fn rand_seed() -> u32 {
14    let seed = SystemTime::now()
15        .duration_since(UNIX_EPOCH)
16        .unwrap()
17        .as_nanos() as u32;
18    return seed;
19}
20
21/// `n` pseudorandom numbers from [0, 1)
22///
23/// Using the Linear Congruential Generator (LCG)
24/// Lewis, Goodman, and Miller in 1969
25///
26/// C-like Parameters from
27/// [wikipedia](https://en.wikipedia.org/wiki/Linear_congruential_generator)
28///
29pub fn rand(n: u64) -> Vec<f64> {
30    // If the first time and the SEED is not set
31    let _ = SEED.compare_exchange(0, rand_seed(), Ordering::Relaxed, Ordering::Relaxed);
32
33    let a: u32 = 22_695_477;
34    let c: u32 = 1;
35
36    // 32-bit SEED converted to 31-bit
37    let mut x: u32 = SEED.load(Ordering::Relaxed) & 0x7FFFFFFF;
38    let mut r: Vec<f64> = Vec::new();
39    for _ in 0..n {
40        // next seed
41        x = x.wrapping_mul(a).wrapping_add(c) & 0x7FFFFFFF;
42        SEED.store(x, Ordering::Relaxed);
43
44        // map to [0 1]
45        // 1. shift the 31-bit by 16 => 15-bit left
46        // 2. max = 2^15 = 32,768
47        // 3. x / max to map to [0 1]
48        r.push((x >> 16) as f64 / 32768.0);
49    }
50    return r;
51}
52
53/// `One` pseudorandom number from [0, 1)
54pub fn rand1() -> f64 {
55    let r = rand(1);
56    return r[0];
57}
58
59/// `n` pseudorandom integers from [a, b]
60pub fn randi(n: u64, start: i64, end: i64) -> Vec<i64> {
61    let min = cmp::min(start, end) as i64;
62    let max = cmp::max(start, end) as i64 + 1;
63    let length = (max - min) as f64;
64    let r = rand(n);
65    return r
66        .into_iter()
67        .map(|x| min + (length * x).floor() as i64)
68        .collect();
69}
70
71/// `One` pseudorandom integer from [a, b]
72pub fn randi1(start: i64, end: i64) -> i64 {
73    return randi(1, start, end)[0];
74}
75
76#[cfg(test)]
77mod tests {
78    use super::*;
79    use std::collections::HashSet;
80
81    /// Check one pseudorandom number is between 0 (inclusive) and 1 (exlusive)
82    #[test]
83    fn check_rand1() {
84        let n = 1000_000;
85        for _ in 0..n {
86            let r = rand1();
87            assert!(r >= 0.);
88            assert!(r < 1.);
89        }
90    }
91
92    /// Check the mean (mu) and variance (var) of pseudorandom numbers
93    #[test]
94    fn check_rand() {
95        let r: Vec<f64> = rand(1000_000);
96        let count: usize = r.iter().len();
97        let mu: f64 = r.iter().sum::<f64>() / count as f64;
98        let var = r
99            .iter()
100            .map(|v| {
101                let distance = *v - mu;
102                distance * distance
103            })
104            .sum::<f64>()
105            / count as f64;
106
107        // mean error
108        let er_mu = (0.5 - (mu * 1e4).round() / 1e4).abs();
109        assert!(er_mu < 0.01);
110
111        // variance error (theoretical variance: 1/12 ~ 0.0833)
112        let er_var = (1. / 12. - (var * 1e4).round() / 1e4).abs();
113        assert!(er_var < 0.001);
114    }
115
116    /// Check the set of pseudorandom integers contains every integer in [a, b], i.e. inclusive of a and b
117    #[test]
118    fn check_randi() {
119        let n = 1000_000;
120        let a = -2;
121        let b = 3;
122        let r = randi(n, a, b);
123        let set_test: HashSet<i64> = r.iter().clone().into_iter().map(|x| *x).collect();
124        let set_true: HashSet<i64> = HashSet::from([-2, -1, 0, 1, 2, 3]);
125        assert_eq!(set_test, set_true);
126    }
127
128    /// Check one pseudorandom integer is between the start and end args
129    #[test]
130    fn check_randi1() {
131        let a = -10;
132        let b = 10;
133        let n = 10_000;
134        for _ in 0..n {
135            let r = randi1(a, b);
136            assert!((a..=b).contains(&r));
137        }
138    }
139}
quantre/random/rng_uniform.rs

quantre/random/
rng_uniform.rs
