1#[derive(Debug)]
4pub struct SimRng {
5 state: u64,
6}
7const FNV_OFFSET: u64 = 14_695_981_039_346_656_037;
10const FNV_PRIME: u64 = 1_099_511_628_211;
11const BOX_MULLER_RADIUS_FACTOR: f64 = 2.0;
12
13fn fnv1a_bytes(mut hash: u64, bytes: &[u8]) -> u64 {
14 for &b in bytes {
15 hash ^= u64::from(b);
16 hash = hash.wrapping_mul(FNV_PRIME);
17 }
18
19 hash
20}
21
22fn derive_seed(seed_base: u64, chunk_id: &str, iteration_index: u32) -> u64 {
23 let mut h = FNV_OFFSET;
24
25 h = fnv1a_bytes(h, &seed_base.to_le_bytes());
26 h = fnv1a_bytes(h, chunk_id.as_bytes());
27 h = fnv1a_bytes(h, &iteration_index.to_le_bytes());
28
29 h
30}
31
32#[must_use]
35pub fn seed_prefix(seed_base: u64, chunk_id: &str) -> u64 {
36 let mut hash = FNV_OFFSET;
37 hash = fnv1a_bytes(hash, &seed_base.to_le_bytes());
38 fnv1a_bytes(hash, chunk_id.as_bytes())
39}
40impl SimRng {
43 fn new_state(seed: u64) -> Self {
44 Self {
45 state: if seed == 0 { 1 } else { seed },
46 }
47 }
48
49 #[must_use]
53 pub fn for_iteration(seed_base: u64, chunk_id: &str, iteration_index: u32) -> Self {
54 Self::new_state(derive_seed(seed_base, chunk_id, iteration_index))
55 }
56
57 #[must_use]
59 pub fn from_prefix(prefix: u64, iteration_index: u32) -> Self {
60 Self::new_state(fnv1a_bytes(prefix, &iteration_index.to_le_bytes()))
61 }
62
63 pub fn next_u64(&mut self) -> u64 {
65 const SHIFT_A: u32 = 13;
66 const SHIFT_B: u32 = 7;
67 const SHIFT_C: u32 = 17;
68 let mut x = self.state;
69
70 x ^= x << SHIFT_A;
71 x ^= x >> SHIFT_B;
72 x ^= x << SHIFT_C;
73 self.state = x;
74
75 x
76 }
77
78 pub fn next_f64(&mut self) -> f64 {
80 const DISCARD_BITS: u32 = 11;
81 const U53_SCALE: f64 = 1.0 / 9_007_199_254_740_992.0;
82 const U32_RANGE: f64 = 4_294_967_296.0;
83
84 let sample = self.next_u64() >> DISCARD_BITS;
85 let [low_0, low_1, low_2, low_3, high_0, high_1, high_2, high_3] = sample.to_le_bytes();
86 let low = u32::from_le_bytes([low_0, low_1, low_2, low_3]);
87 let high = u32::from_le_bytes([high_0, high_1, high_2, high_3]);
88
89 (f64::from(high) * U32_RANGE + f64::from(low)) * U53_SCALE
90 }
91
92 pub fn next_bool_at_prob(&mut self, p: f64) -> bool {
94 self.next_f64() < p.clamp(0.0, 1.0)
95 }
96
97 pub fn range_u32(&mut self, minimum: u32, maximum: u32) -> u32 {
99 if minimum >= maximum {
100 return minimum;
101 }
102
103 let width = u64::from(maximum - minimum);
104 let rejection_start = u64::MAX - u64::MAX % width;
105
106 loop {
107 let sample = self.next_u64();
108
109 if sample < rejection_start {
110 let Ok(offset) = u32::try_from(sample % width) else {
111 return minimum;
112 };
113
114 return minimum + offset;
115 }
116 }
117 }
118
119 pub fn gaussian(&mut self, mean: f64, standard_deviation: f64) -> f64 {
121 if standard_deviation <= 0.0 {
122 return mean;
123 }
124
125 let radius =
126 (-BOX_MULLER_RADIUS_FACTOR * self.next_f64().max(f64::MIN_POSITIVE).ln()).sqrt();
127 let angle = std::f64::consts::TAU * self.next_f64();
128
129 mean + standard_deviation * radius * angle.cos()
130 }
131}
132
133#[cfg(test)]
134mod tests {
135 use googletest::prelude::*;
136
137 use super::*;
138
139 #[gtest]
140 fn for_iteration_known_answer_sequence() -> Result<()> {
141 let mut rng = SimRng::for_iteration(42, "chunk-0", 0);
142 let u64s: Vec<u64> = (0..8).map(|_| rng.next_u64()).collect();
143 let f64s: Vec<f64> = (0..4).map(|_| rng.next_f64()).collect();
144
145 let golden_u64: Vec<u64> = vec![
146 7_908_305_473_855_710_578,
147 1_476_394_694_489_775_416,
148 9_147_714_132_149_237_762,
149 3_278_087_791_654_346_354,
150 246_832_694_459_627_214,
151 13_679_168_970_414_221_987,
152 12_652_401_039_501_625_782,
153 2_456_717_250_748_426_789,
154 ];
155
156 verify_that!(u64s, container_eq(golden_u64))?;
157
158 verify_that!(
159 f64s,
160 elements_are![
161 near(0.362_907_512_845_637_27, 1e-15),
162 near(0.141_259_249_568_543_85, 1e-15),
163 near(0.628_055_221_170_333_3, 1e-15),
164 near(0.246_435_189_170_538_75, 1e-15),
165 ]
166 )
167 }
168
169 #[gtest]
170 fn for_iteration_is_deterministic() -> Result<()> {
171 let mut a = SimRng::for_iteration(42, "chunk-0", 7);
172 let mut b = SimRng::for_iteration(42, "chunk-0", 7);
173 let stream_a: Vec<_> = (0..32).map(|_| a.next_u64()).collect();
174 let stream_b: Vec<_> = (0..32).map(|_| b.next_u64()).collect();
175
176 verify_that!(stream_a, container_eq(stream_b))
177 }
178
179 #[gtest]
180 fn different_iteration_indices_diverge() -> Result<()> {
181 let mut a = SimRng::for_iteration(42, "chunk-0", 0);
182 let mut b = SimRng::for_iteration(42, "chunk-0", 1);
183 let stream_a: Vec<u64> = (0..8).map(|_| a.next_u64()).collect();
184 let stream_b: Vec<u64> = (0..8).map(|_| b.next_u64()).collect();
185
186 verify_that!(stream_a, not(eq(&stream_b)))
187 }
188
189 #[gtest]
190 fn from_prefix_matches_for_iteration() -> Result<()> {
191 let prefix = seed_prefix(99, "alpha");
192 let mut from_prefix = SimRng::from_prefix(prefix, 5);
193 let mut full = SimRng::for_iteration(99, "alpha", 5);
194
195 let prefixed_stream: Vec<_> = (0..16).map(|_| from_prefix.next_u64()).collect();
196 let full_stream: Vec<_> = (0..16).map(|_| full.next_u64()).collect();
197
198 verify_that!(prefixed_stream, container_eq(full_stream))
199 }
200
201 #[gtest]
202 fn next_f64_is_in_unit_interval() -> Result<()> {
203 let mut rng = SimRng::for_iteration(1, "chunk", 0);
204
205 for _ in 0..10_000 {
206 let v = rng.next_f64();
207
208 verify_that!(v, all!(ge(0.0), lt(1.0)))?;
209 }
210
211 Ok(())
212 }
213
214 #[gtest]
215 fn next_bool_clamps_probability_and_consumes_one_draw() -> Result<()> {
216 let mut below_zero = SimRng::for_iteration(7, "draw-count", 3);
217 let mut at_zero = SimRng::for_iteration(7, "draw-count", 3);
218
219 verify_that!(below_zero.next_bool_at_prob(-1.0), eq(false))?;
220 verify_that!(at_zero.next_bool_at_prob(0.0), eq(false))?;
221 verify_that!(below_zero.next_u64(), eq(at_zero.next_u64()))?;
222
223 let mut above_one = SimRng::for_iteration(7, "draw-count", 3);
224 let mut at_one = SimRng::for_iteration(7, "draw-count", 3);
225
226 verify_that!(above_one.next_bool_at_prob(2.0), eq(true))?;
227 verify_that!(at_one.next_bool_at_prob(1.0), eq(true))?;
228
229 verify_that!(above_one.next_u64(), eq(at_one.next_u64()))
230 }
231
232 #[gtest]
233 fn integer_range_is_deterministic_and_respects_half_open_bounds() -> Result<()> {
234 let mut first = SimRng::for_iteration(42, "heartbeat", 7);
235 let mut second = SimRng::for_iteration(42, "heartbeat", 7);
236 let first_stream: Vec<_> = (0..100_000).map(|_| first.range_u32(1, 5_249)).collect();
237 let second_stream: Vec<_> = (0..100_000).map(|_| second.range_u32(1, 5_249)).collect();
238
239 verify_true!(
240 second_stream
241 .iter()
242 .all(|sample| (1..5_249).contains(sample))
243 )?;
244 verify_that!(first_stream, container_eq(second_stream))?;
245
246 verify_that!(first.range_u32(9, 9), eq(9))
247 }
248
249 #[gtest]
250 fn gaussian_stream_is_deterministic_with_expected_distribution() -> Result<()> {
251 const SAMPLE_COUNT: u32 = 100_000;
252 const EXPECTED_MEAN: f64 = 5_250.0;
253 const EXPECTED_STANDARD_DEVIATION: f64 = 25.0;
254
255 let mut first = SimRng::for_iteration(42, "heartbeat", 9);
256 let mut second = SimRng::for_iteration(42, "heartbeat", 9);
257 let samples: Vec<_> = (0..SAMPLE_COUNT)
258 .map(|_| first.gaussian(EXPECTED_MEAN, EXPECTED_STANDARD_DEVIATION))
259 .collect();
260 let repeated: Vec<_> = (0..SAMPLE_COUNT)
261 .map(|_| second.gaussian(EXPECTED_MEAN, EXPECTED_STANDARD_DEVIATION))
262 .collect();
263 let sample_count = f64::from(SAMPLE_COUNT);
264 let mean = samples.iter().sum::<f64>() / sample_count;
265 let variance = samples
266 .iter()
267 .map(|sample| (sample - mean).powi(2))
268 .sum::<f64>()
269 / sample_count;
270
271 verify_that!(samples, container_eq(repeated))?;
272 verify_that!(mean, near(EXPECTED_MEAN, 0.25))?;
273 verify_that!(variance.sqrt(), near(EXPECTED_STANDARD_DEVIATION, 0.25))?;
274
275 verify_that!(first.gaussian(123.0, 0.0), eq(123.0))
276 }
277
278 #[gtest]
279 fn zero_seed_falls_back_to_one() -> Result<()> {
280 let mut rng = SimRng::from_prefix(0, 0);
281 let _ = rng.next_u64();
282 let _ = rng.next_u64();
283
284 verify_that!(rng.next_u64(), not(eq(0)))
285 }
286}