Skip to main content

sci_form/distgeom/embedding/
rng.rs

1/// Trait for RNG sources used in distance geometry embedding.
2/// Both MinstdRand and Mt19937 implement this.
3pub trait DistGeomRng {
4    fn next_double(&mut self) -> f64;
5}
6
7pub struct MinstdRand {
8    state: u64,
9}
10
11impl MinstdRand {
12    const A: u64 = 48271;
13    const M: u64 = 2147483647; // 2^31 - 1
14
15    pub fn new(seed: u32) -> Self {
16        Self { state: seed as u64 }
17    }
18
19    pub fn get_state(&self) -> u64 {
20        self.state
21    }
22
23    /// Returns next integer in [1, M-1]
24    fn next_int(&mut self) -> u64 {
25        self.state = (self.state * Self::A) % Self::M;
26        self.state
27    }
28
29    /// Returns next double in [0, 1) matching boost::uniform_real_distribution<>(0.0, 1.0)
30    /// Formula: (val - min()) / (max() - min() + 1) = (val - 1) / (M - 1)
31    /// where min()=1, max()=M-1, M=2^31-1
32    pub fn next_double(&mut self) -> f64 {
33        let val = self.next_int();
34        (val as f64 - 1.0) / (Self::M as f64 - 1.0)
35    }
36}
37
38/// Mersenne Twister 19937 matching boost::mt19937.
39/// Used as the main RNG for distance picking and coordinate generation.
40pub struct Mt19937 {
41    state: [u32; 624],
42    index: usize,
43}
44
45impl Mt19937 {
46    const N: usize = 624;
47    const M: usize = 397;
48
49    pub fn new(seed: u32) -> Self {
50        let mut mt = Mt19937 {
51            state: [0u32; 624],
52            index: 624,
53        };
54        mt.state[0] = seed;
55        for i in 1..Self::N {
56            mt.state[i] = 1812433253u32
57                .wrapping_mul(mt.state[i - 1] ^ (mt.state[i - 1] >> 30))
58                .wrapping_add(i as u32);
59        }
60        mt
61    }
62
63    fn twist(&mut self) {
64        for i in 0..Self::N {
65            let y = (self.state[i] & 0x80000000) | (self.state[(i + 1) % Self::N] & 0x7fffffff);
66            self.state[i] = self.state[(i + Self::M) % Self::N] ^ (y >> 1);
67            if y & 1 != 0 {
68                self.state[i] ^= 0x9908b0df;
69            }
70        }
71        self.index = 0;
72    }
73
74    /// Generate next u32 value
75    pub fn next_u32(&mut self) -> u32 {
76        if self.index >= Self::N {
77            self.twist();
78        }
79        let mut y = self.state[self.index];
80        self.index += 1;
81
82        y ^= y >> 11;
83        y ^= (y << 7) & 0x9d2c5680;
84        y ^= (y << 15) & 0xefc60000;
85        y ^= y >> 18;
86        y
87    }
88
89    /// Generate next double in [0, 1) matching boost::uniform_real<>(0, 1)
90    /// Formula: mt_output / (mt_max + 1.0) = mt_output / 4294967296.0
91    pub fn next_double(&mut self) -> f64 {
92        self.next_u32() as f64 / 4294967296.0
93    }
94}
95
96impl DistGeomRng for MinstdRand {
97    fn next_double(&mut self) -> f64 {
98        self.next_double()
99    }
100}
101
102impl DistGeomRng for Mt19937 {
103    fn next_double(&mut self) -> f64 {
104        self.next_double()
105    }
106}