Skip to main content

ferray_random/bitgen/
sfc64.rs

1// ferray-random: SFC64 (Small Fast Counting) BitGenerator implementation
2//
3// Sebastian Vigna's SFC64 — a 4-state 64-bit chaotic counting generator.
4// Period >= 2^64 (typical: 2^255). Used by NumPy as the `SFC64`
5// BitGenerator since 1.17.
6
7#![allow(clippy::unreadable_literal)]
8
9use super::BitGenerator;
10
11/// SFC64 (Small Fast Counting) pseudo-random number generator.
12///
13/// Four-word 64-bit state generator by Sebastian Vigna. Very fast,
14/// statistically robust on standard test suites. Period at least
15/// 2^64 from any starting state and typically 2^255.
16pub struct Sfc64 {
17    a: u64,
18    b: u64,
19    c: u64,
20    counter: u64,
21}
22
23impl Sfc64 {
24    #[inline]
25    fn advance(&mut self) -> u64 {
26        let tmp = self.a.wrapping_add(self.b).wrapping_add(self.counter);
27        self.counter = self.counter.wrapping_add(1);
28        self.a = self.b ^ (self.b >> 11);
29        self.b = self.c.wrapping_add(self.c << 3);
30        self.c = self.c.rotate_left(24).wrapping_add(tmp);
31        tmp
32    }
33}
34
35impl BitGenerator for Sfc64 {
36    fn state_bytes(&self) -> Result<Vec<u8>, ferray_core::FerrayError> {
37        let mut out = Vec::with_capacity(32);
38        out.extend_from_slice(&self.a.to_le_bytes());
39        out.extend_from_slice(&self.b.to_le_bytes());
40        out.extend_from_slice(&self.c.to_le_bytes());
41        out.extend_from_slice(&self.counter.to_le_bytes());
42        Ok(out)
43    }
44
45    fn set_state_bytes(&mut self, bytes: &[u8]) -> Result<(), ferray_core::FerrayError> {
46        if bytes.len() != 32 {
47            return Err(ferray_core::FerrayError::invalid_value(format!(
48                "Sfc64 state must be 32 bytes, got {}",
49                bytes.len()
50            )));
51        }
52        self.a = u64::from_le_bytes(bytes[0..8].try_into().unwrap());
53        self.b = u64::from_le_bytes(bytes[8..16].try_into().unwrap());
54        self.c = u64::from_le_bytes(bytes[16..24].try_into().unwrap());
55        self.counter = u64::from_le_bytes(bytes[24..32].try_into().unwrap());
56        Ok(())
57    }
58
59    fn next_u64(&mut self) -> u64 {
60        self.advance()
61    }
62
63    fn seed_from_u64(seed: u64) -> Self {
64        // NumPy's SFC64 init: set a = b = c = seed, counter = 1, then
65        // advance 12 times to "warm up" the state. We follow the
66        // upstream reference exactly.
67        let mut rng = Self {
68            a: seed,
69            b: seed,
70            c: seed,
71            counter: 1,
72        };
73        for _ in 0..12 {
74            let _ = rng.advance();
75        }
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 Sfc64 {
89    fn clone(&self) -> Self {
90        Self {
91            a: self.a,
92            b: self.b,
93            c: self.c,
94            counter: self.counter,
95        }
96    }
97}
98
99#[cfg(test)]
100mod tests {
101    use super::*;
102
103    #[test]
104    fn deterministic_output() {
105        let mut a = Sfc64::seed_from_u64(0xc0ffee);
106        let mut b = Sfc64::seed_from_u64(0xc0ffee);
107        for _ in 0..1000 {
108            assert_eq!(a.next_u64(), b.next_u64());
109        }
110    }
111
112    #[test]
113    fn different_seeds_differ() {
114        let mut a = Sfc64::seed_from_u64(1);
115        let mut b = Sfc64::seed_from_u64(2);
116        let mut diff = false;
117        for _ in 0..100 {
118            if a.next_u64() != b.next_u64() {
119                diff = true;
120                break;
121            }
122        }
123        assert!(diff);
124    }
125
126    #[test]
127    fn full_range() {
128        let mut rng = Sfc64::seed_from_u64(0x1234_5678_9abc_def0);
129        let mut hi = false;
130        let mut lo = false;
131        for _ in 0..10_000 {
132            let v = rng.next_u64();
133            if v > u64::MAX / 2 {
134                hi = true;
135            } else {
136                lo = true;
137            }
138            if hi && lo {
139                break;
140            }
141        }
142        assert!(hi && lo);
143    }
144}