vrd/
quasirandom.rs

1// Copyright © 2023-2026 vrd. All rights reserved.
2// SPDX-License-Identifier: Apache-2.0 OR MIT
3
4//! Quasi-random (low-discrepancy) sequences.
5//!
6//! Quasi-random sequences cover the unit cube `[0, 1)^D` more evenly
7//! than uniform PRNG draws - variance for Monte Carlo integration
8//! scales `O((log n)^d / n)` rather than `O(1/√n)`. They're standard
9//! tools for ray-tracing, financial simulation, and high-dimensional
10//! numerical integration.
11//!
12//! Three constructions are shipped:
13//!
14//! - [`VanDerCorputSequence`](crate::quasirandom::VanDerCorputSequence)
15//!   - 1-D, any prime base.
16//! - [`HaltonSequence`](crate::quasirandom::HaltonSequence)
17//!   - multi-dim, Van der Corput across the first primes. Supports
18//!   up to 32 dimensions.
19//! - [`SobolSequence`](crate::quasirandom::SobolSequence)
20//!   - multi-dim, uses precomputed direction numbers. Supports up
21//!   to 6 dimensions out of the box; the Joe-Kuo D6 file would
22//!   extend this past 21 000.
23//!
24//! These are **not** PRNGs and live alongside, not inside,
25//! [`crate::Random`]. Use a PRNG when you want unpredictability;
26//! use a quasi-random sequence when you want even coverage.
27
28#[cfg(feature = "alloc")]
29use alloc::vec::Vec;
30
31// ---------------------------------------------------------------------------
32// Van der Corput (1-D)
33// ---------------------------------------------------------------------------
34
35/// 1-D Van der Corput sequence in a given prime base. Each step
36/// reverses the index's base-`b` digit expansion and reads it as
37/// a fraction in `[0, 1)`.
38///
39/// # Examples
40///
41/// ```
42/// use vrd::quasirandom::VanDerCorputSequence;
43///
44/// let mut vdc = VanDerCorputSequence::new(2);
45/// assert!((vdc.next_point() - 0.5).abs() < 1e-12);
46/// assert!((vdc.next_point() - 0.25).abs() < 1e-12);
47/// ```
48#[derive(Clone, Copy, Debug)]
49pub struct VanDerCorputSequence {
50    /// Prime base; typically 2 or 3.
51    base: u32,
52    /// 1-indexed position in the sequence.
53    i: u64,
54}
55
56impl VanDerCorputSequence {
57    /// Builds a Van der Corput sequence in base `base`. `base` must
58    /// be ≥ 2 (typically a prime for good distribution).
59    pub fn new(base: u32) -> Self {
60        assert!(base >= 2, "Van der Corput base must be at least 2");
61        Self { base, i: 1 }
62    }
63
64    /// Skips forward by `n` points. Cheap because each point's value
65    /// is determined entirely by the index.
66    pub fn skip(&mut self, n: u64) {
67        self.i = self.i.saturating_add(n);
68    }
69
70    /// Returns the next point in the sequence.
71    pub fn next_point(&mut self) -> f64 {
72        let p = radical_inverse(self.i, self.base);
73        self.i += 1;
74        p
75    }
76}
77
78impl Iterator for VanDerCorputSequence {
79    type Item = f64;
80
81    fn next(&mut self) -> Option<f64> {
82        Some(self.next_point())
83    }
84}
85
86/// Radical inverse - reverses `n`'s digit expansion in `base`
87/// and reads it as a fraction in `[0, 1)`. Powers the Van der
88/// Corput and Halton sequences.
89#[inline]
90fn radical_inverse(mut n: u64, base: u32) -> f64 {
91    let b = u64::from(base);
92    let mut q = 0.0_f64;
93    let mut bk = 1.0_f64 / (base as f64);
94    while n > 0 {
95        let digit = n % b;
96        q += (digit as f64) * bk;
97        n /= b;
98        bk /= base as f64;
99    }
100    q
101}
102
103// ---------------------------------------------------------------------------
104// Halton (multi-dim)
105// ---------------------------------------------------------------------------
106
107/// First 32 primes - bases for the Halton sequence's per-dimension
108/// Van der Corput components.
109const HALTON_PRIMES: [u32; 32] = [
110    2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61,
111    67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131,
112];
113
114/// Maximum number of dimensions supported by [`HaltonSequence`].
115pub const HALTON_MAX_DIM: usize = HALTON_PRIMES.len();
116
117/// Multi-dimensional Halton sequence - Van der Corput across the
118/// first [`HALTON_MAX_DIM`] primes. Discrepancy is excellent for
119/// dimensions up to ~10; beyond that consider [`SobolSequence`] or
120/// scrambled Halton.
121///
122/// # Examples
123///
124/// ```
125/// use vrd::quasirandom::HaltonSequence;
126///
127/// let mut h = HaltonSequence::new(2);
128/// let p = h.next_point::<2>();
129/// assert_eq!(p.len(), 2);
130/// assert!(p[0] >= 0.0 && p[0] < 1.0);
131/// ```
132#[derive(Clone, Copy, Debug)]
133pub struct HaltonSequence {
134    /// Number of dimensions; `1..=HALTON_MAX_DIM`.
135    d: usize,
136    /// 1-indexed position in the sequence.
137    i: u64,
138}
139
140impl HaltonSequence {
141    /// Builds an `d`-dimensional Halton sequence. Panics if
142    /// `d > HALTON_MAX_DIM`.
143    pub fn new(d: usize) -> Self {
144        assert!(
145            d > 0 && d <= HALTON_MAX_DIM,
146            "Halton dimension must be in 1..={HALTON_MAX_DIM}"
147        );
148        Self { d, i: 1 }
149    }
150
151    /// Skips forward by `n` points.
152    pub fn skip(&mut self, n: u64) {
153        self.i = self.i.saturating_add(n);
154    }
155
156    /// Returns the next `D`-dimensional point. Panics if
157    /// `D != self.dimensions()`.
158    #[rustfmt::skip]
159    pub fn next_point<const D: usize>(&mut self) -> [f64; D] {
160        assert_eq!(D, self.d, "const D must match dimensions configured at construction");
161        let mut out = [0.0; D];
162        for (k, slot) in out.iter_mut().enumerate() {
163            *slot = radical_inverse(self.i, HALTON_PRIMES[k]);
164        }
165        self.i += 1;
166        out
167    }
168
169    /// Returns the next point as an allocated [`Vec`]. Useful when
170    /// `D` isn't a compile-time constant. Requires `alloc`.
171    #[cfg(feature = "alloc")]
172    pub fn next_point_vec(&mut self) -> Vec<f64> {
173        let mut out = Vec::with_capacity(self.d);
174        for &base in &HALTON_PRIMES[..self.d] {
175            out.push(radical_inverse(self.i, base));
176        }
177        self.i += 1;
178        out
179    }
180
181    /// Returns the dimension count.
182    pub fn dimensions(&self) -> usize {
183        self.d
184    }
185}
186
187// ---------------------------------------------------------------------------
188// Sobol (multi-dim)
189// ---------------------------------------------------------------------------
190
191/// Maximum number of dimensions [`SobolSequence`] supports with the
192/// shipped direction-number table.
193pub const SOBOL_MAX_DIM: usize = 6;
194
195/// Direction numbers for the first six Sobol dimensions
196/// (Bratley-Fox 1988 starter set). Each row holds 32 direction
197/// numbers, sized for `u32` indexing.
198// Hand-derived from polynomial primitives & m_i initial values
199// listed in Bratley & Fox (1988) Table 1. Computed via:
200//   v_i = m_i << (32 - i)               for i = 1..=s
201//   v_i = v_{i-s} ^ (v_{i-s} >> s)
202//                ^ a_1 v_{i-1} ^ ... ^ a_{s-1} v_{i-s+1}   for i > s
203// where (a_1, ..., a_{s-1}) are the polynomial's middle bits
204// and `s` is its degree.
205const SOBOL_DIRECTIONS: [[u32; 32]; SOBOL_MAX_DIM] = [
206    [
207        0x8000_0000,
208        0x4000_0000,
209        0x2000_0000,
210        0x1000_0000,
211        0x0800_0000,
212        0x0400_0000,
213        0x0200_0000,
214        0x0100_0000,
215        0x0080_0000,
216        0x0040_0000,
217        0x0020_0000,
218        0x0010_0000,
219        0x0008_0000,
220        0x0004_0000,
221        0x0002_0000,
222        0x0001_0000,
223        0x0000_8000,
224        0x0000_4000,
225        0x0000_2000,
226        0x0000_1000,
227        0x0000_0800,
228        0x0000_0400,
229        0x0000_0200,
230        0x0000_0100,
231        0x0000_0080,
232        0x0000_0040,
233        0x0000_0020,
234        0x0000_0010,
235        0x0000_0008,
236        0x0000_0004,
237        0x0000_0002,
238        0x0000_0001,
239    ],
240    // Dim 2: polynomial x + 1, m_1 = 1
241    [
242        0x8000_0000,
243        0xC000_0000,
244        0xA000_0000,
245        0xF000_0000,
246        0x8800_0000,
247        0xCC00_0000,
248        0xAA00_0000,
249        0xFF00_0000,
250        0x8080_0000,
251        0xC0C0_0000,
252        0xA0A0_0000,
253        0xF0F0_0000,
254        0x8888_0000,
255        0xCCCC_0000,
256        0xAAAA_0000,
257        0xFFFF_0000,
258        0x8000_8000,
259        0xC000_C000,
260        0xA000_A000,
261        0xF000_F000,
262        0x8800_8800,
263        0xCC00_CC00,
264        0xAA00_AA00,
265        0xFF00_FF00,
266        0x8080_8080,
267        0xC0C0_C0C0,
268        0xA0A0_A0A0,
269        0xF0F0_F0F0,
270        0x8888_8888,
271        0xCCCC_CCCC,
272        0xAAAA_AAAA,
273        0xFFFF_FFFF,
274    ],
275    // Dim 3: polynomial x^2 + x + 1, m = [1, 3]
276    [
277        0x8000_0000,
278        0x4000_0000,
279        0xE000_0000,
280        0x5000_0000,
281        0xF800_0000,
282        0x2C00_0000,
283        0xE200_0000,
284        0x4900_0000,
285        0xF7C0_0000,
286        0x4A60_0000,
287        0xEF00_0000,
288        0x5DD0_0000,
289        0xF108_0000,
290        0x2528_0000,
291        0xCCC2_0000,
292        0xFFFB_0000,
293        0x8001_8000,
294        0x4002_4000,
295        0xE003_E000,
296        0x5004_5000,
297        0xF807_F800,
298        0x2C0C_2C00,
299        0xE20D_E200,
300        0x4914_4900,
301        0xF7CD_F7C0,
302        0x4A6F_4A60,
303        0xEF7B_EF00,
304        0x5DDF_5DD0,
305        0xF108_F108,
306        0x2528_2528,
307        0xCCC2_CCC2,
308        0xFFFB_FFFB,
309    ],
310    // Dim 4: polynomial x^3 + x + 1, m = [1, 3, 1]
311    [
312        0x8000_0000,
313        0x4000_0000,
314        0x6000_0000,
315        0xD000_0000,
316        0x6800_0000,
317        0x3400_0000,
318        0xC600_0000,
319        0x6D00_0000,
320        0x3580_0000,
321        0xC340_0000,
322        0x6FA0_0000,
323        0xB510_0000,
324        0xC628_0000,
325        0xA714_0000,
326        0xC32A_0000,
327        0x6431_0000,
328        0xB2E1_8000,
329        0xC55A_4000,
330        0xEFBA_6000,
331        0xF563_D000,
332        0xD27F_6800,
333        0xE3B1_3400,
334        0xC568_C600,
335        0xB78F_6D00,
336        0xC44B_3580,
337        0x8261_C340,
338        0x4123_6FA0,
339        0x60D2_B510,
340        0xD026_C628,
341        0x6817_A714,
342        0x3412_C32A,
343        0xC600_6431,
344    ],
345    // Dim 5: polynomial x^3 + x^2 + 1, m = [1, 1, 3]
346    [
347        0x8000_0000,
348        0xC000_0000,
349        0x2000_0000,
350        0xB000_0000,
351        0xC800_0000,
352        0x6C00_0000,
353        0x4600_0000,
354        0x6700_0000,
355        0xC580_0000,
356        0xC1C0_0000,
357        0x2820_0000,
358        0x9870_0000,
359        0x9C68_0000,
360        0x4CB4_0000,
361        0xC65A_0000,
362        0x69C7_0000,
363        0xC32E_8000,
364        0xC09B_C000,
365        0x2031_E000,
366        0xB008_F000,
367        0xC8C4_4800,
368        0x6C66_E400,
369        0x4647_2A00,
370        0x6726_F900,
371        0xC5A6_82C0,
372        0xC176_4860,
373        0x2832_BA90,
374        0x987F_7CB0,
375        0x9C7E_3C50,
376        0x4C24_CE40,
377        0xC60B_E5A0,
378        0x69CC_99C0,
379    ],
380    // Dim 6: polynomial x^4 + x + 1, m = [1, 1, 1, 3]
381    [
382        0x8000_0000,
383        0xC000_0000,
384        0xA000_0000,
385        0x9000_0000,
386        0xD800_0000,
387        0xFC00_0000,
388        0xCE00_0000,
389        0xD900_0000,
390        0xFD80_0000,
391        0xCFC0_0000,
392        0xD9E0_0000,
393        0xFCD0_0000,
394        0xCE38_0000,
395        0xD974_0000,
396        0xFD46_0000,
397        0xCFEF_0000,
398        0xD9E8_8000,
399        0xFCD0_4000,
400        0xCE38_E000,
401        0xD974_9000,
402        0xFD46_D800,
403        0xCFEF_FC00,
404        0xD9E8_CE00,
405        0xFCD0_D900,
406        0xCE38_FD80,
407        0xD974_CFC0,
408        0xFD46_D9E0,
409        0xCFEF_FCD0,
410        0xD9E8_CE38,
411        0xFCD0_D974,
412        0xCE38_FD46,
413        0xD974_CFEF,
414    ],
415];
416
417/// Multi-dimensional Sobol sequence (Bratley-Fox 1988 starter set,
418/// dimensions 1–6). For higher-dimensional uses, switch to a Joe-Kuo
419/// D6 direction-number table (out of scope here).
420///
421/// # Examples
422///
423/// ```
424/// use vrd::quasirandom::SobolSequence;
425///
426/// let mut s = SobolSequence::new(2);
427/// let p = s.next_point::<2>();
428/// assert!(p[0] >= 0.0 && p[0] < 1.0);
429/// ```
430#[derive(Clone, Copy, Debug)]
431pub struct SobolSequence {
432    /// Number of dimensions; `1..=SOBOL_MAX_DIM`.
433    d: usize,
434    /// 0-indexed position. `i == 0` returns the canonical
435    /// origin point `(0, 0, ..., 0)`; later steps XOR into `x`.
436    i: u64,
437    /// Per-dimension running state. XOR'd with the appropriate
438    /// direction number on every step; converted to `f64` via
439    /// the 2³² scale to yield the output point.
440    x: [u32; SOBOL_MAX_DIM],
441}
442
443impl SobolSequence {
444    /// Builds a `d`-dimensional Sobol sequence. Panics if `d` is 0
445    /// or exceeds [`SOBOL_MAX_DIM`].
446    pub fn new(d: usize) -> Self {
447        assert!(
448            d > 0 && d <= SOBOL_MAX_DIM,
449            "Sobol dimension must be in 1..={SOBOL_MAX_DIM}"
450        );
451        Self {
452            d,
453            i: 0,
454            x: [0; SOBOL_MAX_DIM],
455        }
456    }
457
458    /// Skips forward by `n` points.
459    pub fn skip(&mut self, n: u64) {
460        for _ in 0..n {
461            self.advance();
462        }
463    }
464
465    /// Returns the next `D`-dimensional point.
466    #[rustfmt::skip]
467    pub fn next_point<const D: usize>(&mut self) -> [f64; D] {
468        assert_eq!(D, self.d, "const D must match dimensions configured at construction");
469        let scale = 1.0_f64 / ((1u64 << 32) as f64);
470        let mut out = [0.0; D];
471        if self.i == 0 {
472            // Convention: first point is (0, 0, …, 0).
473            self.i = 1;
474            return out;
475        }
476        let bit = trailing_zero_bit(self.i);
477        assert!(bit < 32, "Sobol exhausted: only 2^32 points per dimension");
478        for k in 0..D {
479            self.x[k] ^= SOBOL_DIRECTIONS[k][bit];
480            out[k] = (self.x[k] as f64) * scale;
481        }
482        self.i += 1;
483        out
484    }
485
486    /// Returns the next point as an allocated [`Vec`]. Requires
487    /// `alloc`.
488    #[cfg(feature = "alloc")]
489    pub fn next_point_vec(&mut self) -> Vec<f64> {
490        let scale = 1.0_f64 / ((1u64 << 32) as f64);
491        if self.i == 0 {
492            self.i = 1;
493            return alloc::vec![0.0; self.d];
494        }
495        let bit = trailing_zero_bit(self.i);
496        assert!(
497            bit < 32,
498            "Sobol exhausted: only 2^32 points per dimension"
499        );
500        let mut out = Vec::with_capacity(self.d);
501        for (k, dir) in SOBOL_DIRECTIONS.iter().enumerate().take(self.d)
502        {
503            self.x[k] ^= dir[bit];
504            out.push((self.x[k] as f64) * scale);
505        }
506        self.i += 1;
507        out
508    }
509
510    /// Steps the internal state forward by one position
511    /// without producing an output point. Used by
512    /// [`Self::skip`] to amortise long forward jumps.
513    fn advance(&mut self) {
514        if self.i == 0 {
515            self.i = 1;
516            return;
517        }
518        let bit = trailing_zero_bit(self.i);
519        if bit >= 32 {
520            return;
521        }
522        for (k, dir) in SOBOL_DIRECTIONS.iter().enumerate().take(self.d)
523        {
524            self.x[k] ^= dir[bit];
525        }
526        self.i += 1;
527    }
528
529    /// Returns the dimension count.
530    pub fn dimensions(&self) -> usize {
531        self.d
532    }
533
534    /// Sets the internal index. Test-only - used to exercise the
535    /// `bit >= 32` defensive branch in [`Self::advance`] without
536    /// running 2³² actual iterations.
537    #[cfg(test)]
538    pub(crate) fn set_i_for_test(&mut self, i: u64) {
539        self.i = i;
540    }
541}
542
543/// Returns the number of trailing zero bits in `n`. Sobol's
544/// next-bit-to-toggle is `trailing_zeros(i + 1)`; this is the
545/// well-known O(1) closed form.
546#[inline]
547fn trailing_zero_bit(n: u64) -> usize {
548    n.trailing_zeros() as usize
549}
550
551#[cfg(test)]
552mod tests {
553    use super::*;
554
555    #[test]
556    fn vdc_base_2_matches_classical() {
557        let mut vdc = VanDerCorputSequence::new(2);
558        let expected = [0.5, 0.25, 0.75, 0.125, 0.625, 0.375, 0.875];
559        for &e in &expected {
560            let got = vdc.next_point();
561            assert!((got - e).abs() < 1e-12, "got {got}, want {e}");
562        }
563    }
564
565    #[test]
566    fn vdc_in_unit_interval() {
567        let mut vdc = VanDerCorputSequence::new(3);
568        for _ in 0..256 {
569            let p = vdc.next_point();
570            assert!((0.0..1.0).contains(&p), "point {p} out of [0, 1)");
571        }
572    }
573
574    #[test]
575    fn halton_2d_in_unit_square() {
576        let mut h = HaltonSequence::new(2);
577        for _ in 0..1024 {
578            let p = h.next_point::<2>();
579            assert!((0.0..1.0).contains(&p[0]));
580            assert!((0.0..1.0).contains(&p[1]));
581        }
582    }
583
584    #[test]
585    fn halton_dim_panic_on_overflow() {
586        // Should not panic at construction with the max dim.
587        let _ = HaltonSequence::new(HALTON_MAX_DIM);
588    }
589
590    #[test]
591    fn sobol_2d_first_few_points() {
592        let mut s = SobolSequence::new(2);
593        // (0, 0) is the canonical starting point.
594        assert_eq!(s.next_point::<2>(), [0.0, 0.0]);
595        let p = s.next_point::<2>();
596        assert!(p[0] > 0.0 && p[0] < 1.0);
597        assert!(p[1] > 0.0 && p[1] < 1.0);
598    }
599
600    #[test]
601    fn sobol_iterator_in_unit_square() {
602        let mut s = SobolSequence::new(2);
603        for _ in 0..512 {
604            let p = s.next_point::<2>();
605            assert!((0.0..1.0).contains(&p[0]));
606            assert!((0.0..1.0).contains(&p[1]));
607        }
608    }
609
610    /// `skip(n)` on Van der Corput advances by `n` positions.
611    /// `VanDerCorputSequence` also impls `Iterator`, so we must
612    /// disambiguate from `Iterator::skip` (which is consuming).
613    #[test]
614    fn vdc_skip_advances_correctly() {
615        let mut a = VanDerCorputSequence::new(2);
616        let mut b = VanDerCorputSequence::new(2);
617        for _ in 0..5 {
618            let _ = a.next_point();
619        }
620        VanDerCorputSequence::skip(&mut b, 5);
621        assert_eq!(a.next_point(), b.next_point());
622    }
623
624    /// Van der Corput's `Iterator` impl returns the same value as
625    /// `next_point()`.
626    #[test]
627    fn vdc_iterator_matches_next_point() {
628        let mut a = VanDerCorputSequence::new(2);
629        let mut b = VanDerCorputSequence::new(2);
630        for _ in 0..8 {
631            assert_eq!(a.next_point(), b.next().unwrap());
632        }
633    }
634
635    /// `HaltonSequence::skip` + `dimensions` accessors.
636    #[test]
637    fn halton_skip_and_dimensions() {
638        let mut a = HaltonSequence::new(3);
639        let mut b = HaltonSequence::new(3);
640        assert_eq!(a.dimensions(), 3);
641        for _ in 0..7 {
642            let _ = a.next_point::<3>();
643        }
644        b.skip(7);
645        assert_eq!(a.next_point::<3>(), b.next_point::<3>());
646    }
647
648    /// `HaltonSequence::next_point_vec` returns the heap-allocated
649    /// equivalent of `next_point::<D>()`.
650    #[test]
651    #[cfg(feature = "alloc")]
652    fn halton_next_point_vec_matches_const_generic() {
653        let mut a = HaltonSequence::new(3);
654        let mut b = HaltonSequence::new(3);
655        let v_const = a.next_point::<3>();
656        let v_heap = b.next_point_vec();
657        assert_eq!(v_heap.len(), 3);
658        for k in 0..3 {
659            assert!((v_const[k] - v_heap[k]).abs() < 1e-12);
660        }
661    }
662
663    /// `SobolSequence::skip(n)` + `dimensions`.
664    #[test]
665    fn sobol_skip_and_dimensions() {
666        let mut a = SobolSequence::new(3);
667        let mut b = SobolSequence::new(3);
668        assert_eq!(b.dimensions(), 3);
669        for _ in 0..7 {
670            let _ = a.next_point::<3>();
671        }
672        b.skip(7);
673        assert_eq!(a.next_point::<3>(), b.next_point::<3>());
674    }
675
676    /// `SobolSequence::next_point_vec`: covers both the `i == 0`
677    /// fast path (returns all-zero vec) and the normal sampling
678    /// path.
679    #[test]
680    #[cfg(feature = "alloc")]
681    fn sobol_next_point_vec_paths() {
682        let mut s = SobolSequence::new(3);
683        let first = s.next_point_vec();
684        assert_eq!(first, alloc::vec![0.0; 3]);
685        let second = s.next_point_vec();
686        assert_eq!(second.len(), 3);
687        assert!(second.iter().any(|&v| v != 0.0));
688    }
689
690    /// `VanDerCorputSequence::new(base)` panics if base < 2.
691    #[test]
692    #[should_panic(expected = "at least 2")]
693    fn vdc_base_below_two_panics() {
694        let _ = VanDerCorputSequence::new(1);
695    }
696
697    /// `HaltonSequence::new(d)` panics for d == 0.
698    #[test]
699    #[should_panic(expected = "must be in")]
700    fn halton_zero_dim_panics() {
701        let _ = HaltonSequence::new(0);
702    }
703
704    /// `HaltonSequence::new(d)` panics for d > HALTON_MAX_DIM.
705    #[test]
706    #[should_panic(expected = "must be in")]
707    fn halton_over_max_dim_panics() {
708        let _ = HaltonSequence::new(HALTON_MAX_DIM + 1);
709    }
710
711    /// `SobolSequence::new(d)` panics for d == 0.
712    #[test]
713    #[should_panic(expected = "must be in")]
714    fn sobol_zero_dim_panics() {
715        let _ = SobolSequence::new(0);
716    }
717
718    /// `SobolSequence::new(d)` panics for d > SOBOL_MAX_DIM.
719    #[test]
720    #[should_panic(expected = "must be in")]
721    fn sobol_over_max_dim_panics() {
722        let _ = SobolSequence::new(SOBOL_MAX_DIM + 1);
723    }
724
725    /// `Sobol::advance` early-returns when `i == 0` (first call).
726    /// Cover by calling `skip(1)` on a fresh sequence: that triggers
727    /// `advance()` with i==0, which sets i=1 and returns.
728    #[test]
729    fn sobol_skip_from_zero() {
730        let mut a = SobolSequence::new(2);
731        let mut b = SobolSequence::new(2);
732        // a does the canonical first call.
733        let _ = a.next_point::<2>();
734        // b skips one - should put it in the same state.
735        b.skip(1);
736        assert_eq!(a.next_point::<2>(), b.next_point::<2>());
737    }
738
739    /// Covers the `bit >= 32` defensive branch in `advance()`.
740    /// In practice this would only fire after 2³² advances; the
741    /// test harness sets `i` directly to `1 << 32` so the next
742    /// `advance()` triggers it.
743    #[test]
744    fn sobol_advance_handles_exhaustion() {
745        let mut s = SobolSequence::new(2);
746        s.set_i_for_test(1u64 << 32);
747        // Skip one - internally calls advance(), which hits the
748        // `bit >= 32` branch and returns without mutating x.
749        let x_before = s.x;
750        s.skip(1);
751        assert_eq!(
752            s.x, x_before,
753            "exhausted advance must not mutate state"
754        );
755    }
756
757    /// Monte Carlo π convergence: Sobol's error after N points
758    /// should shrink faster than Halton's, both faster than a fair
759    /// uniform PRNG would. With N=4096 we're not aiming for tight
760    /// thresholds - just verifying that all three converge into
761    /// the right neighbourhood.
762    #[test]
763    fn quasi_pi_estimates_converge() {
764        const N: usize = 4096;
765        let mut h = HaltonSequence::new(2);
766        let mut s = SobolSequence::new(2);
767
768        let mut h_inside = 0usize;
769        let mut s_inside = 0usize;
770        for _ in 0..N {
771            let p = h.next_point::<2>();
772            if p[0] * p[0] + p[1] * p[1] < 1.0 {
773                h_inside += 1;
774            }
775            let q = s.next_point::<2>();
776            if q[0] * q[0] + q[1] * q[1] < 1.0 {
777                s_inside += 1;
778            }
779        }
780        let pi_h = 4.0 * (h_inside as f64) / (N as f64);
781        let pi_s = 4.0 * (s_inside as f64) / (N as f64);
782        assert!(
783            (pi_h - core::f64::consts::PI).abs() < 0.05,
784            "Halton π = {pi_h}"
785        );
786        assert!(
787            (pi_s - core::f64::consts::PI).abs() < 0.05,
788            "Sobol π = {pi_s}"
789        );
790    }
791}
vrd/quasirandom.rs

vrd/
quasirandom.rs
