ferray_random/bitgen/
pcg64dxsm.rs

1// ferray-random: PCG64DXSM BitGenerator implementation
2//
3// PCG-DXSM 128/64 (cheap-multiply variant of PCG-XSL-RR) — recommended
4// upgrade over PCG64 by Melissa O'Neill in 2022. NumPy added it as
5// `PCG64DXSM` in 1.21. Period 2^128.
6
7#![allow(clippy::unreadable_literal)]
8
9use super::BitGenerator;
10
11/// PCG64DXSM (DXSM output function) pseudo-random number generator.
12///
13/// 128-bit linear congruential generator with the DXSM (double-xorshift
14/// multiply) output function. Period 2^128. Recommended over the
15/// classic XSL-RR `Pcg64` for new applications because the DXSM
16/// permutation has stronger empirical bit-level randomness.
17pub struct Pcg64Dxsm {
18    state: u128,
19    inc: u128,
20}
21
22const PCG_CHEAP_MULTIPLIER: u64 = 0xda94_2042_e4dd_58b5;
23
24impl Pcg64Dxsm {
25    #[inline]
26    const fn step(&mut self) {
27        // Cheap-multiplier LCG advance: state = state * mult + inc, where
28        // mult is a 64-bit constant lifted to 128 bits.
29        let mult128 = PCG_CHEAP_MULTIPLIER as u128;
30        self.state = self.state.wrapping_mul(mult128).wrapping_add(self.inc);
31    }
32
33    #[inline]
34    const fn output(state: u128) -> u64 {
35        // DXSM permutation:
36        //   hi = state >> 64; lo = state | 1
37        //   hi ^= hi >> 32
38        //   hi *= cheap_mult
39        //   hi ^= hi >> 48
40        //   hi *= lo
41        let hi0 = (state >> 64) as u64;
42        let lo = (state as u64) | 1;
43        let mut hi = hi0;
44        hi ^= hi >> 32;
45        hi = hi.wrapping_mul(PCG_CHEAP_MULTIPLIER);
46        hi ^= hi >> 48;
47        hi = hi.wrapping_mul(lo);
48        hi
49    }
50}
51
52impl BitGenerator for Pcg64Dxsm {
53    fn state_bytes(&self) -> Result<Vec<u8>, ferray_core::FerrayError> {
54        let mut out = Vec::with_capacity(32);
55        out.extend_from_slice(&self.state.to_le_bytes());
56        out.extend_from_slice(&self.inc.to_le_bytes());
57        Ok(out)
58    }
59
60    fn set_state_bytes(&mut self, bytes: &[u8]) -> Result<(), ferray_core::FerrayError> {
61        if bytes.len() != 32 {
62            return Err(ferray_core::FerrayError::invalid_value(format!(
63                "Pcg64Dxsm state must be 32 bytes, got {}",
64                bytes.len()
65            )));
66        }
67        self.state = u128::from_le_bytes(bytes[0..16].try_into().unwrap());
68        let inc = u128::from_le_bytes(bytes[16..32].try_into().unwrap());
69        if inc & 1 == 0 {
70            return Err(ferray_core::FerrayError::invalid_value(
71                "Pcg64Dxsm inc must be odd",
72            ));
73        }
74        self.inc = inc;
75        Ok(())
76    }
77
78    fn next_u64(&mut self) -> u64 {
79        let old = self.state;
80        self.step();
81        Self::output(old)
82    }
83
84    fn seed_from_u64(seed: u64) -> Self {
85        let seed128 = {
86            let mut s = seed;
87            let a = super::splitmix64(&mut s);
88            let b = super::splitmix64(&mut s);
89            ((a as u128) << 64) | (b as u128)
90        };
91        let inc = {
92            let mut s = seed.wrapping_add(0xda3e_39cb_94b9_5bdb);
93            let a = super::splitmix64(&mut s);
94            let b = super::splitmix64(&mut s);
95            (((a as u128) << 64) | (b as u128)) | 1
96        };
97        let mut rng = Self { state: 0, inc };
98        rng.step();
99        rng.state = rng.state.wrapping_add(seed128);
100        rng.step();
101        rng
102    }
103
104    fn jump(&mut self) -> Option<()> {
105        None
106    }
107
108    fn stream(_seed: u64, _stream_id: u64) -> Option<Self> {
109        None
110    }
111}
112
113impl Clone for Pcg64Dxsm {
114    fn clone(&self) -> Self {
115        Self {
116            state: self.state,
117            inc: self.inc,
118        }
119    }
120}
121
122#[cfg(test)]
123mod tests {
124    use super::*;
125
126    #[test]
127    fn deterministic_output() {
128        let mut a = Pcg64Dxsm::seed_from_u64(7);
129        let mut b = Pcg64Dxsm::seed_from_u64(7);
130        for _ in 0..1000 {
131            assert_eq!(a.next_u64(), b.next_u64());
132        }
133    }
134
135    #[test]
136    fn different_seeds_differ() {
137        let mut a = Pcg64Dxsm::seed_from_u64(7);
138        let mut b = Pcg64Dxsm::seed_from_u64(8);
139        let mut diff = false;
140        for _ in 0..100 {
141            if a.next_u64() != b.next_u64() {
142                diff = true;
143                break;
144            }
145        }
146        assert!(diff);
147    }
148
149    #[test]
150    fn full_range() {
151        let mut rng = Pcg64Dxsm::seed_from_u64(0xfeed_face);
152        let mut hi = false;
153        let mut lo = false;
154        for _ in 0..10_000 {
155            let v = rng.next_u64();
156            if v > u64::MAX / 2 {
157                hi = true;
158            } else {
159                lo = true;
160            }
161            if hi && lo {
162                break;
163            }
164        }
165        assert!(hi && lo);
166    }
167}
ferray_random/bitgen/pcg64dxsm.rs

ferray_random/bitgen/
pcg64dxsm.rs
