hoomd_rand/
sfc.rs

1// Copyright (c) 2024-2026 The Regents of the University of Michigan.
2// Part of hoomd-rs, released under the BSD 3-Clause License.
3
4use core::convert::Infallible;
5
6use rand::{
7    SeedableRng,
8    rand_core::{TryRng, utils},
9};
10use serde::{Deserialize, Serialize};
11
12use crate::util::read_le_u64;
13
14/// The "Small Fast Chaotic" PRNG, originally designed by Chris Doty-Humphrey.
15///
16/// This PRNG holds 256 bits of state (including a 64-bit counter) and generates 64 bits
17/// of output with each step. The minimum cycle length is $`2^{64}`$, and the expected
18/// period is $`~2^{255}`$. Independent seeds are guaranteed to not collide within the
19/// first $`2^{64}`$ steps, although the actual time to collision may be much longer.
20///
21/// This specific implementation passes [PractRand] with >2TB of output for each of
22/// several low-entropy seeds, and >1TB for sets of 8 and 64 parallel streams whose seed
23/// differs by only a single bit.
24///
25/// ## Seeding (construction)
26///
27/// This generator implements the [`SeedableRng`] trait. Any method may be used,
28/// and the results are guaranteed to be [portable].
29///
30/// Seeding the generator with `seed_from_u64` is often the most convenient and
31/// guarantees good pseudorandom statistics.
32/// ```
33/// use hoomd_rand::SFC64;
34/// use rand::SeedableRng;
35/// let rng = SFC64::seed_from_u64(42);
36/// ```
37/// See also [Seeding RNGs] in the Rust Rand book.
38///
39/// ## Generation
40///
41/// The generators implements [`TryRng`] and thus also `Rng`.
42/// See also the [Random Values] chapter in the Rust Rand book.
43///
44/// [portable]: https://rust-random.github.io/book/crate-reprod.html
45/// [Seeding RNGs]: https://rust-random.github.io/book/guide-seeding.html
46/// [Random Values]: https://rust-random.github.io/book/guide-values.html
47/// [PractRand]: https://pracrand.sourceforge.net
48
49#[derive(Clone, Debug, PartialEq, Eq, Serialize, Deserialize)]
50pub struct SFC64 {
51    /// The internal state of the PRNG.
52    state: [u64; 3],
53    /// The current step of the PRNG. This value is incremented with each output, but
54    /// can be set independently without issues.
55    counter: u64,
56}
57
58impl SFC64 {
59    /// Set up the PRNG from four u64 values, discarding the first 12 outputs to avoid
60    /// correlations between similar seeds.
61    #[inline]
62    pub(crate) fn initialize(a: u64, b: u64, c: u64, counter: u64) -> Self {
63        let mut rng = Self {
64            state: [a, b, c],
65            counter,
66        };
67        (0..12).for_each(|_| {
68            rng.step();
69        });
70        rng
71    }
72    /// Increment the RNG forward, returning the generated result.
73    #[inline]
74    fn step(&mut self) -> u64 {
75        const BARREL_SHIFT: u32 = 24;
76        const RSHIFT: u64 = 11;
77        const LSHIFT: u64 = 3;
78        let out = self.state[0]
79            .wrapping_add(self.state[1])
80            .wrapping_add(self.counter);
81        self.state[0] = self.state[1] ^ (self.state[1] >> RSHIFT);
82        self.state[1] = self.state[2].wrapping_add(self.state[2] << LSHIFT);
83        self.state[2] = self.state[2].rotate_left(BARREL_SHIFT).wrapping_add(out);
84        self.counter = self.counter.wrapping_add(1); // Weyl increment of 1
85        out
86    }
87}
88
89impl TryRng for SFC64 {
90    type Error = Infallible;
91
92    #[inline]
93    fn try_next_u64(&mut self) -> Result<u64, Self::Error> {
94        Ok(self.step())
95    }
96    #[inline]
97    #[expect(
98        clippy::cast_possible_truncation,
99        reason = "the truncation is intended"
100    )]
101    fn try_next_u32(&mut self) -> Result<u32, Self::Error> {
102        Ok(self.step() as u32)
103    }
104    #[inline]
105    fn try_fill_bytes(&mut self, dst: &mut [u8]) -> Result<(), Self::Error> {
106        utils::fill_bytes_via_next_word(dst, || Ok(self.step()))
107    }
108}
109impl SeedableRng for SFC64 {
110    type Seed = [u8; 24];
111    #[inline]
112    fn from_seed(seed: Self::Seed) -> Self {
113        let seed = &mut seed.as_slice();
114        Self::initialize(read_le_u64(seed), read_le_u64(seed), read_le_u64(seed), 0)
115    }
116    #[inline]
117    fn seed_from_u64(state: u64) -> Self {
118        Self::initialize(state, 0, 0, 1)
119    }
120}
121
122#[cfg(test)]
123mod tests {
124    use super::*;
125    use rand::{Rng, SeedableRng};
126    use rstest::rstest;
127
128    // To generate test data:
129    // ```python
130    // from numpy.random import SFC64
131    // rng = SFC64()
132    // rng.state = {
133    // "bit_generator": "SFC64",
134    // "state": {"state":np.array([0,0,0,1], dtype=np.uint64)},
135    // "has_uint32": 0,
136    // "uinteger": 0
137    // }
138    // Generator(rng).integers(
139    //     np.uint64(2**64-1), size=30, dtype=np.uint64, endpoint=True
140    // )
141    //
142    // array([                   1,                    2,                   12,
143    //                   150994975,     2533275243454595,     7601064794802825,
144    //           81067317779357880,   939904496648135862,  8358283255545778679,
145    //         3766502869832372394,  3749490695777847966, 10585829324311968283,
146    // # Done warming up =====================================================
147    //         4237781876154851393, 17705428440413258140,  1322197197711907681,
148    //          822724228132957142,  2474202602039083746,  5912426283212852001,
149    //        15821317571115833757, 10375476962501160791, 16772721532102950986,
150    //        10344472337605944266, 17980158625826311727,  5068054200821024954,
151    //         7387740087731720141,  3233965506304800125,   805889469043559033,
152    //         5406730423683535818, 15071469935640508508,  4156074611580516626],
153    //       dtype=uint64)
154
155    // ```
156
157    #[rustfmt::skip]
158    const SFC64_REFERENCE_OUTPUT: [u64; 17] = [
159        4_237_781_876_154_851_393,  17_705_428_440_413_258_140,
160        1_322_197_197_711_907_681,     822_724_228_132_957_142,
161        2_474_202_602_039_083_746,   5_912_426_283_212_852_001,
162        15_821_317_571_115_833_757, 10_375_476_962_501_160_791,
163        16_772_721_532_102_950_986, 10_344_472_337_605_944_266,
164        17_980_158_625_826_311_727,  5_068_054_200_821_024_954,
165        7_387_740_087_731_720_141,   3_233_965_506_304_800_125,
166        805_889_469_043_559_033,     5_406_730_423_683_535_818,
167        15_071_469_935_640_508_508
168    ];
169
170    #[rstest::fixture]
171    fn large_uniform_sample() -> Vec<u64> {
172        const N: u32 = 2_u32.pow(23);
173        let mut rng = SFC64::initialize(456_981, 0xcafe, 9_345_663_908, 123_456_789);
174        (0..N).map(|_| rng.next_u64()).collect::<Vec<_>>()
175    }
176
177    #[rstest::fixture]
178    fn zeros(large_uniform_sample: Vec<u64>) -> u32 {
179        large_uniform_sample
180            .iter()
181            .fold(0, |acc, x| acc + x.count_zeros())
182    }
183    #[rstest::fixture]
184    fn ones(large_uniform_sample: Vec<u64>) -> u32 {
185        large_uniform_sample
186            .iter()
187            .fold(0, |acc, x| acc + x.count_ones())
188    }
189
190    #[test]
191    fn test_sfc64_against_numpy() {
192        let mut rng = SFC64::seed_from_u64(0);
193        (0..17).for_each(|i| {
194            assert_eq!(
195                rng.next_u64(),
196                SFC64_REFERENCE_OUTPUT[i],
197                "failed at index {i}"
198            );
199        });
200    }
201
202    /// This test is quite slow, but tests that our generation of the ~2^17th value
203    /// matches the reference impl
204    ///
205    /// NOTE: set state to [0,0,0,1], then
206    /// `rg.integers(np.uint64(2**64-1), size=2**(17)+12, dtype=np.uint64, endpoint=True)[-1]`
207    #[test]
208    fn test_sfc64_deep() {
209        let mut rng = SFC64::seed_from_u64(0);
210        (0..(2_u64.pow(17) - 1)).for_each(|_| {
211            rng.next_u64();
212        });
213        assert_eq!(rng.next_u64(), 4_977_758_738_274_538_201);
214    }
215
216    #[rstest]
217    fn test_uniformity_mean(large_uniform_sample: Vec<u64>, zeros: u32, ones: u32) {
218        let standard_error = ((4 * 64 * large_uniform_sample.len()) as f64)
219            .sqrt()
220            .recip();
221        approxim::assert_abs_diff_eq!(
222            f64::from(zeros) / f64::from(ones),
223            1.0,
224            epsilon = standard_error
225        );
226    }
227    #[rstest]
228    fn test_uniformity_variance(large_uniform_sample: Vec<u64>, ones: u32) {
229        let n_bits: u128 = large_uniform_sample.len() as u128 * 64;
230
231        let variance = (u128::from(ones) * (n_bits - u128::from(ones))) as f64
232            / (n_bits * (n_bits - 1)) as f64;
233        let epsilon = (n_bits as f64).sqrt().recip(); // approximate
234        approxim::assert_abs_diff_eq!(variance, 0.25, epsilon = epsilon);
235    }
236
237    #[rstest]
238    fn test_autocorrelation(
239        large_uniform_sample: Vec<u64>,
240        #[values(1, 2, 4, 8, 64, 256, 65536)] lag: usize,
241    ) {
242        let n = large_uniform_sample.len();
243        let sample_f64: Vec<f64> = large_uniform_sample.iter().map(|&x| x as f64).collect();
244
245        let mean = sample_f64.iter().sum::<f64>() / n as f64;
246        let var = sample_f64.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / n as f64;
247        let cov = (0..n - lag)
248            .map(|i| (sample_f64[i] - mean) * (sample_f64[i + lag] - mean))
249            .sum::<f64>()
250            / (n - lag) as f64;
251        let autocorr = cov / var;
252
253        let epsilon = 3.0 / (n as f64).sqrt();
254        approxim::assert_abs_diff_eq!(autocorr, 0.0, epsilon = epsilon);
255    }
256}
hoomd_rand/sfc.rs

hoomd_rand/
sfc.rs
