conspire/math/random/
mod.rs

1use std::{
2    cell::Cell,
3    f64::consts::TAU,
4    time::{SystemTime, UNIX_EPOCH},
5};
6
7thread_local! {
8    static STATE: Cell<u64> = const { Cell::new(0) };
9}
10
11fn seed() -> u64 {
12    let now = SystemTime::now()
13        .duration_since(UNIX_EPOCH)
14        .unwrap_or_default();
15    let t = now.as_nanos() as u64;
16    let x = 0u8;
17    let addr = (&x as *const u8 as usize) as u64;
18    let mut s = t ^ addr.wrapping_mul(0x9E3779B97F4A7C15);
19    if s == 0 {
20        s = 1;
21    }
22    s
23}
24
25fn next_u64() -> u64 {
26    STATE.with(|st| {
27        let mut s = st.get();
28        if s == 0 {
29            s = seed();
30        }
31        s ^= s >> 12;
32        s ^= s << 25;
33        s ^= s >> 27;
34        st.set(s);
35        s.wrapping_mul(0x2545F4914F6CDD1D)
36    })
37}
38
39fn get_random() -> u8 {
40    (next_u64() >> 56) as u8
41}
42
43pub fn random_u8(max: u8) -> u8 {
44    if max == u8::MAX {
45        return get_random();
46    }
47    let bound = (max as u16) + 1;
48    let threshold = (256u16 / bound) * bound;
49    loop {
50        let v = get_random() as u16;
51        if v < threshold {
52            return (v % bound) as u8;
53        }
54    }
55}
56
57pub fn random_u64() -> u64 {
58    next_u64()
59}
60
61pub fn random_uniform() -> f64 {
62    let x = next_u64() >> 11;
63    (x as f64) * (1.0 / ((1u64 << 53) as f64))
64}
65
66thread_local! {
67    static NORMAL_SPARE: Cell<Option<f64>> = const { Cell::new(None) };
68}
69
70fn random_normal_standard() -> f64 {
71    NORMAL_SPARE.with(|spare| {
72        if let Some(z) = spare.take() {
73            return z;
74        }
75        let mut u1 = random_uniform();
76        while u1 <= 0.0 {
77            u1 = random_uniform();
78        }
79        let u2 = random_uniform();
80        let r = (-2.0 * u1.ln()).sqrt();
81        let (s, c) = (TAU * u2).sin_cos();
82        let z0 = r * c;
83        let z1 = r * s;
84        spare.set(Some(z1));
85        z0
86    })
87}
88
89pub fn random_normal(mean: f64, std: f64) -> f64 {
90    mean + std * random_normal_standard()
91}
92
93fn random_exp1() -> f64 {
94    let mut u = random_uniform();
95    while u <= 0.0 {
96        u = random_uniform();
97    }
98    -u.ln()
99}
100
101fn random_gamma_k3_scale1() -> f64 {
102    random_exp1() + random_exp1() + random_exp1()
103}
104
105pub fn random_x2_normal(mean: f64, std: f64) -> f64 {
106    let m = mean / std;
107    let z_star = if m >= -1.0 { m + 1.0 } else { 0.0 };
108    let h_min = 0.5 * (z_star - m).powi(2) - z_star;
109    loop {
110        let z = random_gamma_k3_scale1();
111        let h = 0.5 * (z - m).powi(2) - z;
112        let acceptance_probability = (-(h - h_min)).exp();
113        if random_uniform() < acceptance_probability {
114            return std * z;
115        }
116    }
117}
conspire/math/random/mod.rs

conspire/math/random/
mod.rs
