roberts_sequence/lib.rs
1//! An implementation of the Roberts quasirandom sequence
2//!
3//! See <https://extremelearning.com.au/unreasonable-effectiveness-of-quasirandom-sequences/>, <https://www.martysmods.com/a-better-r2-sequence/>
4use std::array;
5
6/// Compute the value of the `DIM`-dimensional Roberts sequence at `index`
7///
8/// `shift` is a vector of elements in [0, 1) specifying a toroidal shift that can be used to
9/// generate distinct (though still highly correlated) sequences. This is useful to e.g. improve
10/// temporal accumulation without requiring large indexes.
11///
12/// Numerical limitations degrade the quality of output as `index` grows. For most applications,
13/// this effect should be insignificant for sequences that end before one million.
14///
15/// Guaranteed to be deterministic for DIM <= 4.
16pub fn r<const DIM: usize>(shift: [f64; DIM], index: u32) -> [f64; DIM] {
17 let inv_phi = 1.0 / phi(DIM);
18 let mut alpha_accum = inv_phi;
19 array::from_fn(|i| {
20 // Compute alpha by repeated multiplication to avoid nondeterminism in powi
21 let alpha = alpha_accum;
22 alpha_accum *= inv_phi;
23 // Improve precision when `index` is large per
24 // https://www.martysmods.com/a-better-r2-sequence/
25 let alpha = 1.0f64 - alpha;
26 // LLVM should optimize `alpha` into a constant in release builds
27 (shift[i] + alpha * (index + 1) as f64) % 1.0
28 })
29}
30
31/// Generalized golden ratio
32fn phi(dim: usize) -> f64 {
33 match dim {
34 1 => std::f64::consts::GOLDEN_RATIO,
35 2 => 1.324717957244746, // Plastic ratio
36 3 => 1.2207440846057596,
37 4 => 1.1673039782614187,
38 // Add more hard-coded cases to guarantee higher-dimensional determinism
39 _ => {
40 // General case, used to compute the above deterministic cases. Not guaranteed to be
41 // deterministic due to `powf`.
42 let mut x = 2.0;
43 for _ in 0..30 {
44 x = (x + 1.0f64).powf(1.0 / ((dim + 1) as f64))
45 }
46 x
47 }
48 }
49}
50
51#[cfg(test)]
52mod tests {
53 use super::*;
54
55 #[test]
56 fn determinism() {
57 assert_eq!(
58 array::from_fn::<_, 3, _>(|i| r::<1>([0.0], i as u32)[0]),
59 [0.3819660112501052, 0.7639320225002104, 0.14589803375031574]
60 );
61 }
62}