Skip to main content

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 next_u64(&mut self) -> u64 {
54        let old = self.state;
55        self.step();
56        Self::output(old)
57    }
58
59    fn seed_from_u64(seed: u64) -> Self {
60        let seed128 = {
61            let mut s = seed;
62            let a = super::splitmix64(&mut s);
63            let b = super::splitmix64(&mut s);
64            ((a as u128) << 64) | (b as u128)
65        };
66        let inc = {
67            let mut s = seed.wrapping_add(0xda3e_39cb_94b9_5bdb);
68            let a = super::splitmix64(&mut s);
69            let b = super::splitmix64(&mut s);
70            (((a as u128) << 64) | (b as u128)) | 1
71        };
72        let mut rng = Self { state: 0, inc };
73        rng.step();
74        rng.state = rng.state.wrapping_add(seed128);
75        rng.step();
76        rng
77    }
78
79    fn jump(&mut self) -> Option<()> {
80        None
81    }
82
83    fn stream(_seed: u64, _stream_id: u64) -> Option<Self> {
84        None
85    }
86}
87
88impl Clone for Pcg64Dxsm {
89    fn clone(&self) -> Self {
90        Self {
91            state: self.state,
92            inc: self.inc,
93        }
94    }
95}
96
97#[cfg(test)]
98mod tests {
99    use super::*;
100
101    #[test]
102    fn deterministic_output() {
103        let mut a = Pcg64Dxsm::seed_from_u64(7);
104        let mut b = Pcg64Dxsm::seed_from_u64(7);
105        for _ in 0..1000 {
106            assert_eq!(a.next_u64(), b.next_u64());
107        }
108    }
109
110    #[test]
111    fn different_seeds_differ() {
112        let mut a = Pcg64Dxsm::seed_from_u64(7);
113        let mut b = Pcg64Dxsm::seed_from_u64(8);
114        let mut diff = false;
115        for _ in 0..100 {
116            if a.next_u64() != b.next_u64() {
117                diff = true;
118                break;
119            }
120        }
121        assert!(diff);
122    }
123
124    #[test]
125    fn full_range() {
126        let mut rng = Pcg64Dxsm::seed_from_u64(0xfeed_face);
127        let mut hi = false;
128        let mut lo = false;
129        for _ in 0..10_000 {
130            let v = rng.next_u64();
131            if v > u64::MAX / 2 {
132                hi = true;
133            } else {
134                lo = true;
135            }
136            if hi && lo {
137                break;
138            }
139        }
140        assert!(hi && lo);
141    }
142}