lds-rs 0.1.10

//! Low-Discrepancy Sequence (LDS) Generator
//!
//! This library implements a set of low-discrepancy sequence generators, which are used to create
//! sequences of numbers that are more evenly distributed than random numbers. These sequences are
//! particularly useful in various fields such as computer graphics, numerical integration, and
//! Monte Carlo simulations.
//!
//! The library defines several structs, each representing a different type of low-discrepancy sequence generator. The main types of sequences implemented are:
//!
//! 1. van der Corput sequence
//! 2. Halton sequence
//! 3. Circle sequence
//! 4. Disk sequence
//! 5. Sphere sequence
//! 6. 3-Sphere Hopf sequence
//! 7. N-dimensional Halton sequence
//!
//! Each generator takes specific inputs, usually in the form of base numbers or sequences of base numbers. These bases determine how the sequences are generated. The generators produce outputs in the form of floating-point numbers or vectors of floating-point numbers, depending on the dimensionality of the sequence.
//!
//! The core algorithm used in most of these generators is the van der Corput sequence. This sequence
//! is created by expressing integers in a given base, reversing the digits, and placing them after a
//! decimal point. For example, in base 2, the sequence would start: 1/2, 1/4, 3/4, 1/8, 5/8, and so on.
//!
//! The Halton sequence extends this concept to multiple dimensions by using a different base for each
//! dimension. The Circle and Sphere sequences use trigonometric functions to map these
//! low-discrepancy sequences onto circular or spherical surfaces.
//!
//! The library also includes utility functions and constants to support these generators. For instance, there's a list of prime numbers that can be used as bases for the sequences.
//!
//! Each generator struct has methods to produce the next value in the sequence (`pop()`) and to reset the sequence to a specific starting point (`reseed()`). This allows for flexible use of the generators in various applications.
//!
//! The purpose of this library is to provide a toolkit for generating well-distributed sequences of
//! numbers, which can be used in place of random numbers in many applications to achieve more uniform
//! coverage of a given space or surface. This can lead to more efficient and accurate results in
//! tasks like sampling, integration, and optimization.
//!
//! ![Halton 2D scatter plot][halton-2d-scatter]
//!
//! *Figure: 500 points of a 2D Halton sequence with bases 2 and 3.*
//!
//! The image above is embedded via `embed-doc-image`. To build docs with images locally:
//!
//! ```bash
//! cargo doc --features doc-images --open
//! ```
//!
//! On docs.rs, images are rendered automatically.
//!
//! # Visual examples
//!
//! | Generator | Visualisation |
//! |---|---|
//! | **Van der Corput** (base 2) — 64 values | ![VdC step sequence][vdc-sequence] |
//! | **Circle** (base 2) — 200 points | ![Circle points][circle-points] |
//! | **Disk** (bases 2, 3) — 200 points | ![Disk points][disk-points] |
//! | **Halton 2D** (bases 2, 3) — 500 points | ![Halton 2D scatter][halton-2d-scatter] |

// Embed images for docs (feature-gated for rustdoc compatibility).
// Image paths are relative to the crate root.
#![cfg_attr(
    feature = "doc-images",
    doc = embed_doc_image::embed_image!("vdc-sequence", "docs/images/vdc-sequence.png")
)]
#![cfg_attr(
    feature = "doc-images",
    doc = embed_doc_image::embed_image!("circle-points", "docs/images/circle-points.png")
)]
#![cfg_attr(
    feature = "doc-images",
    doc = embed_doc_image::embed_image!("disk-points", "docs/images/disk-points.png")
)]
#![cfg_attr(
    feature = "doc-images",
    doc = embed_doc_image::embed_image!("halton-2d-scatter", "docs/images/halton-2d-scatter.png")
)]

use std::f64::consts::PI;
use std::sync::atomic::{AtomicU64, Ordering};

/// Constant for 2π
pub const TWO_PI: f64 = 2.0 * PI;

/// Maximum number of digits for van der Corput sequence
pub const MAX_DIGITS: usize = 64;

/// van der Corput sequence function
///
/// Converts an integer to its radical inverse by reversing its base-$$b$$ expansion:
///
/// $$ n = \sum_{k=0}^{m} d_k b^k \quad\longrightarrow\quad \phi_b(n) = \sum_{k=0}^{m} \frac{d_k}{b^{k+1}} $$
///
/// where $$ d_k \in \{0, 1, \dots, b-1\} $$ are the base-$$b$$ digits of $$n$$.
///
/// # Arguments
///
/// * `count` - The number for which we want to calculate the van der Corput sequence value
/// * `base` - The base of the number system being used (defaults to 2)
///
/// # Examples
///
/// ```
/// use lds_rs::vdc;
/// assert_eq!(vdc(11, 2), 0.8125);
/// ```
pub fn vdc(count: u64, base: u64) -> f64 {
    let mut count = count;
    let mut reslt = 0.0;
    let mut denom = 1.0;
    let base_f64 = base as f64;

    while count != 0 {
        denom *= base_f64;
        let remainder = (count % base) as f64;
        count /= base;
        reslt += remainder / denom;
    }
    reslt
}

/// van der Corput sequence generator
///
/// Generates the van der Corput sequence, a low-discrepancy sequence commonly used in
/// quasi-Monte Carlo methods. The sequence is generated by iterating over a
/// base and calculating the fractional part of the number in that base.
///
/// The radical inverse process reverses the base-$$b$$ digits of the integer.
/// For example, $$n = 13 = 1101_2$$, reversed digits give $$\phi_2(13) = 0.1011_2 = 0.6875$$.
///
/// # Examples
///
/// ```
/// use lds_rs::VdCorput;
/// let mut vgen = VdCorput::new(2);
/// vgen.reseed(0);
/// assert_eq!(vgen.pop(), 0.5);
/// assert_eq!(vgen.pop(), 0.25);
/// assert_eq!(vgen.pop(), 0.75);
/// ```
#[cfg_attr(feature = "doc-images", doc = svgbobdoc::transform!(
/// ```svgbob
///      .──────────────.
///      │  n = 13₁₀    │
///      │  ≡ 1101₂     │
///      '──────┬───────'
///             │ reverse base-2 digits
///             ▼
///      .──────────────.
///      │  reversed:   │
///      │  1011₂       │
///      │  as decimal: │
///      │  0.1011₂     │
///      │  = 0.6875    │
///      '──────────────'
/// ```
))]
#[derive(Debug)]
pub struct VdCorput {
    count: AtomicU64,
    base: u64,
    rev_lst: Vec<f64>,
}

impl VdCorput {
    /// Creates a new van der Corput sequence generator with the given base
    ///
    /// # Arguments
    ///
    /// * `base` - The base of the number system (defaults to 2 if not specified)
    pub fn new(base: u64) -> Self {
        assert!(base >= 2, "base must be >= 2, got {}", base);
        let mut rev_lst = Vec::with_capacity(MAX_DIGITS);
        let mut reverse = 1.0;
        let base_f64 = base as f64;

        for _ in 0..MAX_DIGITS {
            reverse /= base_f64;
            rev_lst.push(reverse);
        }

        Self {
            count: AtomicU64::new(0),
            base,
            rev_lst,
        }
    }

    /// Generates the next value in the sequence
    ///
    /// $$ \phi_b(n) = \sum_{k=0}^{m} \frac{d_k}{b^{k+1}} $$
    ///
    /// Increments the count and calculates the van der Corput sequence value
    /// for that count and base.
    pub fn pop(&mut self) -> f64 {
        let count = self.count.fetch_add(1, Ordering::Relaxed) + 1; // ignore 0
        let mut count = count;
        let mut res = 0.0;
        let mut i = 0;

        while count != 0 {
            let remainder = (count % self.base) as f64;
            count /= self.base;
            if remainder != 0.0 {
                res += remainder * self.rev_lst[i];
            }
            i += 1;
        }
        res
    }

    /// Returns the next value without advancing the state (peek)
    ///
    /// $$ \phi_b(n) = \sum_{k=0}^{m} \frac{d_k}{b^{k+1}} $$
    ///
    /// Allows looking at the next value in the sequence without consuming it.
    pub fn peek(&self) -> f64 {
        let mut count = self.count.load(Ordering::Relaxed) + 1;
        let mut res = 0.0;
        let mut i = 0;

        while count != 0 {
            let remainder = (count % self.base) as f64;
            count /= self.base;
            if remainder != 0.0 {
                res += remainder * self.rev_lst[i];
            }
            i += 1;
        }
        res
    }

    /// Advances the sequence by `n` values without computing them
    ///
    /// # Arguments
    ///
    /// * `n` - The number of values to advance
    pub fn advance(&self, n: u64) {
        self.count.fetch_add(n, Ordering::Relaxed);
    }

    /// Returns the current index (number of values generated so far)
    pub fn get_index(&self) -> u64 {
        self.count.load(Ordering::Relaxed)
    }

    /// Resets the state of the sequence generator to a specific seed value
    ///
    /// # Arguments
    ///
    /// * `seed` - The seed value that determines the starting point of the sequence generation
    pub fn reseed(&mut self, seed: u64) {
        self.count.store(seed, Ordering::Relaxed);
    }
}

impl Iterator for VdCorput {
    type Item = f64;

    /// Returns the next value in the sequence
    ///
    /// This allows VdCorput to be used with iterator methods like `.take()`, `.collect()`, etc.
    fn next(&mut self) -> Option<Self::Item> {
        Some(self.pop())
    }
}

impl Default for VdCorput {
    fn default() -> Self {
        Self::new(2)
    }
}

impl Clone for VdCorput {
    /// Creates a deep copy of the van der Corput generator
    ///
    /// The cloned generator starts with the same internal state as the original,
    /// including the current count, base, and reversed digit list.
    fn clone(&self) -> Self {
        Self {
            count: AtomicU64::new(self.count.load(Ordering::Relaxed)),
            base: self.base,
            rev_lst: self.rev_lst.clone(),
        }
    }
}

/// Halton sequence generator
///
/// Generates points in a 2-dimensional space using the Halton sequence.
/// The Halton sequence is a low-discrepancy sequence that is often used in
/// quasi-Monte Carlo methods. It is generated by iterating over two different
/// bases and calculating the fractional parts of the numbers in those bases.
///
/// # Examples
///
/// ```
/// use lds_rs::Halton;
/// let mut hgen = Halton::new([2, 3]);
/// hgen.reseed(0);
/// let res = hgen.pop();
/// assert_eq!(res[0], 0.5);
/// assert!((res[1] - 1.0/3.0).abs() < 1e-10);
/// ```
#[cfg_attr(feature = "doc-images", doc = svgbobdoc::transform!(
/// ```svgbob
///  .───────────.     .───────────.
///  │ VdC(2, n) │────►│ x_coord   │
///  '───────────'     '───────────'
///  .───────────.     .───────────.
///  │ VdC(3, n) │────►│ y_coord   │
///  '───────────'     '───────────'
/// ```
))]
pub struct Halton {
    vdc0: VdCorput,
    vdc1: VdCorput,
}

impl Halton {
    /// Creates a new Halton sequence generator with the given bases
    ///
    /// # Arguments
    ///
    /// * `base` - An array of two integers used as bases for generating the Halton sequence
    pub fn new(base: [u64; 2]) -> Self {
        Self {
            vdc0: VdCorput::new(base[0]),
            vdc1: VdCorput::new(base[1]),
        }
    }

    /// Generates the next point in the Halton sequence
    ///
    /// Returns the next point in the Halton sequence as a `[f64; 2]`.
    pub fn pop(&mut self) -> [f64; 2] {
        [self.vdc0.pop(), self.vdc1.pop()]
    }

    /// Returns the next point without advancing the state (peek)
    pub fn peek(&self) -> [f64; 2] {
        [self.vdc0.peek(), self.vdc1.peek()]
    }

    /// Skips `n` points in the sequence without computing them
    ///
    /// # Arguments
    ///
    /// * `n` - The number of points to skip
    pub fn advance(&self, n: u64) {
        self.vdc0.advance(n);
        self.vdc1.advance(n);
    }

    /// Returns the current index (number of points generated so far)
    pub fn get_index(&self) -> u64 {
        self.vdc0.get_index()
    }

    /// Resets the state of the sequence generator to a specific seed value
    ///
    /// # Arguments
    ///
    /// * `seed` - The seed value that determines the starting point of the sequence generation
    pub fn reseed(&mut self, seed: u64) {
        self.vdc0.reseed(seed);
        self.vdc1.reseed(seed);
    }
}

impl Iterator for Halton {
    type Item = [f64; 2];

    /// Returns the next point in the Halton sequence
    ///
    /// This allows Halton to be used with iterator methods like `.take()`, `.collect()`, etc.
    fn next(&mut self) -> Option<Self::Item> {
        Some(self.pop())
    }
}

impl Clone for Halton {
    /// Creates a deep copy of the Halton generator
    ///
    /// The cloned generator contains cloned VdCorput instances for each dimension,
    /// preserving the same internal state.
    fn clone(&self) -> Self {
        Self {
            vdc0: self.vdc0.clone(),
            vdc1: self.vdc1.clone(),
        }
    }
}
/// Unit Circle sequence generator
///
/// Generates points on the unit circle using a low-discrepancy sequence.
///
/// Maps a van der Corput value to an angle $$\theta = 2\pi v$$ and returns
/// $$(\cos\theta,\; \sin\theta)$$ on the unit circle $$S^1$$.
///
/// # Examples
///
/// ```
/// use lds_rs::Circle;
/// let mut cgen = Circle::new(2);
/// cgen.reseed(0);
/// let res = cgen.pop();
/// // Should be approximately [-1.0, 0.0] (cos(π), sin(π))
/// assert!((res[0] + 1.0).abs() < 1e-10);
/// assert!(res[1].abs() < 1e-10);
/// ```
#[cfg_attr(feature = "doc-images", doc = svgbobdoc::transform!(
/// ```svgbob
///  .───────────.        .───────────────.
///  │ φ₂(n)     │───────►│ θ = 2π·φ₂(n) │
///  │ ∈ [0, 1)  │   θ    │               │
///  '───────────'        '───────┬───────'
///                               │
///                               ▼
///                        .───────────────.
///                        │ (cos θ, sin θ)│
///                        │  on unit S¹   │
///                        '───────────────'
/// ```
))]
pub struct Circle {
    vdc: VdCorput,
}

impl Circle {
    /// Creates a new Circle sequence generator with the given base
    ///
    /// # Arguments
    ///
    /// * `base` - The base of the van der Corput sequence
    pub fn new(base: u64) -> Self {
        assert!(base >= 2, "base must be >= 2, got {}", base);
        Self {
            vdc: VdCorput::new(base),
        }
    }

    /// Generates the next point on the unit circle
    ///
    /// $$ \theta = 2\pi v, \qquad (\cos\theta,\; \sin\theta) $$
    ///
    /// Maps the van der Corput value $$ v \in \[0,1\] $$ to $$ \[0, 2\pi\] $$.
    ///
    /// Returns the next point on the unit circle as a `[f64; 2]`.
    pub fn pop(&mut self) -> [f64; 2] {
        let theta = self.vdc.pop() * TWO_PI; // map to [0, 2π]
        [theta.cos(), theta.sin()]
    }

    /// Returns the next point without advancing the state (peek)
    ///
    /// $$ \theta = 2\pi v, \qquad (\cos\theta,\; \sin\theta) $$
    pub fn peek(&self) -> [f64; 2] {
        let theta = self.vdc.peek() * TWO_PI;
        [theta.cos(), theta.sin()]
    }

    /// Skips `n` points in the sequence without computing them
    ///
    /// # Arguments
    ///
    /// * `n` - The number of points to skip
    pub fn advance(&self, n: u64) {
        self.vdc.advance(n);
    }

    /// Returns the current index (number of points generated so far)
    pub fn get_index(&self) -> u64 {
        self.vdc.get_index()
    }

    /// Resets the state of the sequence generator to a specific seed value
    ///
    /// # Arguments
    ///
    /// * `seed` - The seed value that determines the starting point of the sequence generation
    pub fn reseed(&mut self, seed: u64) {
        self.vdc.reseed(seed);
    }
}

impl Iterator for Circle {
    type Item = [f64; 2];

    /// Returns the next point on the unit circle
    ///
    /// This allows Circle to be used with iterator methods like `.take()`, `.collect()`, etc.
    fn next(&mut self) -> Option<Self::Item> {
        Some(self.pop())
    }
}

impl Clone for Circle {
    fn clone(&self) -> Self {
        Self {
            vdc: self.vdc.clone(),
        }
    }
}

/// Unit Disk sequence generator
///
/// Generates points in the unit disk using a low-discrepancy sequence.
///
/// Maps two van der Corput values to polar coordinates:
/// $$\theta = 2\pi v_\theta$$ and $$r = \sqrt{v_r}$$, returning
/// $$(r\cos\theta,\; r\sin\theta)$$ inside the unit disk.
///
/// # Examples
///
/// ```
/// use lds_rs::Disk;
/// let mut dgen = Disk::new([2, 3]);
/// dgen.reseed(0);
/// let res = dgen.pop();
/// // First point should be on the unit disk
/// let radius_sq = res[0] * res[0] + res[1] * res[1];
/// assert!(radius_sq <= 1.0);
/// ```
#[cfg_attr(feature = "doc-images", doc = svgbobdoc::transform!(
/// ```svgbob
///  .───────────.   θ = 2π·φ₂(n)
///  │ φ₂(n)     │───► angle θ
///  '───────────'
///  .───────────.   r = √φ₃(n)
///  │ φ₃(n)     │───► radius r
///  '───────────'
///
///       combine:
///   (r·cosθ, r·sinθ)
///         │
///         ▼
///   .───────────────.
///   │ point inside  │
///   │ unit disk D²  │
///   '───────────────'
/// ```
))]
pub struct Disk {
    vdc0: VdCorput,
    vdc1: VdCorput,
}

impl Disk {
    /// Creates a new Disk sequence generator with the given bases
    ///
    /// # Arguments
    ///
    /// * `base` - An array of two integers used as bases for generating the sequence
    pub fn new(base: [u64; 2]) -> Self {
        Self {
            vdc0: VdCorput::new(base[0]),
            vdc1: VdCorput::new(base[1]),
        }
    }

    /// Generates the next point in the unit disk
    ///
    /// $$ \theta = 2\pi v_\theta, \qquad r = \sqrt{v_r}, \qquad (r\cos\theta,\; r\sin\theta) $$
    ///
    /// Returns the next point in the unit disk as a `[f64; 2]`.
    pub fn pop(&mut self) -> [f64; 2] {
        let theta = self.vdc0.pop() * TWO_PI; // map to [0, 2π]
        let radius = self.vdc1.pop().sqrt(); // map to [0, 1]
        [radius * theta.cos(), radius * theta.sin()]
    }

    /// Returns the next point without advancing the state (peek)
    ///
    /// $$ \theta = 2\pi v_\theta, \qquad r = \sqrt{v_r}, \qquad (r\cos\theta,\; r\sin\theta) $$
    pub fn peek(&self) -> [f64; 2] {
        let theta = self.vdc0.peek() * TWO_PI;
        let radius = self.vdc1.peek().sqrt();
        [radius * theta.cos(), radius * theta.sin()]
    }

    /// Skips `n` points in the sequence without computing them
    ///
    /// # Arguments
    ///
    /// * `n` - The number of points to skip
    pub fn advance(&self, n: u64) {
        self.vdc0.advance(n);
        self.vdc1.advance(n);
    }

    /// Returns the current index (number of points generated so far)
    pub fn get_index(&self) -> u64 {
        self.vdc0.get_index()
    }

    /// Resets the state of the sequence generator to a specific seed value
    ///
    /// # Arguments
    ///
    /// * `seed` - The seed value that determines the starting point of the sequence generation
    pub fn reseed(&mut self, seed: u64) {
        self.vdc0.reseed(seed);
        self.vdc1.reseed(seed);
    }
}

impl Iterator for Disk {
    type Item = [f64; 2];

    /// Returns the next point in the unit disk
    ///
    /// This allows Disk to be used with iterator methods like `.take()`, `.collect()`, etc.
    fn next(&mut self) -> Option<Self::Item> {
        Some(self.pop())
    }
}

impl Clone for Disk {
    fn clone(&self) -> Self {
        Self {
            vdc0: self.vdc0.clone(),
            vdc1: self.vdc1.clone(),
        }
    }
}

/// Unit Sphere sequence generator
///
/// Generates points on the unit sphere using a low-discrepancy sequence.
///
/// Uses a cylindrical equal-area projection from two van der Corput values:
/// $$z = 2v_z - 1$$ and $$\theta = 2\pi v_\theta$$, returning
/// $$(\sqrt{1-z^2}\cos\theta,\; \sqrt{1-z^2}\sin\theta,\; z)$$ on $$S^2$$.
///
/// # Examples
///
/// ```
/// use lds_rs::Sphere;
/// let mut sgen = Sphere::new([2, 3]);
/// sgen.reseed(0);
/// let res = sgen.pop();
/// // Should be on the unit sphere
/// let radius_sq = res[0] * res[0] + res[1] * res[1] + res[2] * res[2];
/// assert!((radius_sq - 1.0).abs() < 1e-10);
/// ```
#[cfg_attr(feature = "doc-images", doc = svgbobdoc::transform!(
/// ```svgbob
///  .───────────.   z = 2·φ₂(n) - 1
///  │ φ₂(n)     │───► z ∈ [-1, 1]
///  '───────────'
///  .───────────.   θ = 2π·φ₃(n)
///  │ φ₃(n)     │───► angle θ
///  '───────────'
///
///          │
///          ▼
///  .──────────────────────────.
///  │ ( √(1-z²)·cosθ,          │
///  │   √(1-z²)·sinθ,          │
///  │   z )                    │
///  │ point on unit sphere S²  │
///  '──────────────────────────'
/// ```
))]
pub struct Sphere {
    vdcgen: VdCorput,
    cirgen: Circle,
}

impl Sphere {
    /// Creates a new Sphere sequence generator with the given bases
    ///
    /// # Arguments
    ///
    /// * `base` - An array of two integers used as bases for generating the sequence
    pub fn new(base: [u64; 2]) -> Self {
        Self {
            vdcgen: VdCorput::new(base[0]),
            cirgen: Circle::new(base[1]),
        }
    }

    /// Generates the next point on the unit sphere
    ///
    /// $$ \phi = 2v - 1,\quad v \in \[0,1\], \qquad (\sqrt{1-\phi^2}\cos\theta,\; \sqrt{1-\phi^2}\sin\theta,\; \phi) $$
    ///
    /// where $$ \theta = 2\pi v_\theta $$ comes from the Circle generator.
    /// This is a cylindrical equal-area projection.
    ///
    /// Returns the next point on the unit sphere as a `[f64; 3]`.
    pub fn pop(&mut self) -> [f64; 3] {
        let cosphi = 2.0 * self.vdcgen.pop() - 1.0; // map to [-1, 1]
        let sinphi = (1.0 - cosphi * cosphi).sqrt(); // cylindrical mapping
        let [cos, sin] = self.cirgen.pop();
        [sinphi * cos, sinphi * sin, cosphi]
    }

    /// Returns the next point without advancing the state (peek)
    ///
    /// $$ \phi = 2v - 1, \qquad (\sqrt{1-\phi^2}\cos\theta,\; \sqrt{1-\phi^2}\sin\theta,\; \phi) $$
    pub fn peek(&self) -> [f64; 3] {
        let cosphi = 2.0 * self.vdcgen.peek() - 1.0;
        let sinphi = (1.0 - cosphi * cosphi).sqrt();
        let [cos, sin] = self.cirgen.peek();
        [sinphi * cos, sinphi * sin, cosphi]
    }

    /// Skips `n` points in the sequence without computing them
    ///
    /// # Arguments
    ///
    /// * `n` - The number of points to skip
    pub fn advance(&self, n: u64) {
        self.cirgen.advance(n);
        self.vdcgen.advance(n);
    }

    /// Returns the current index (number of points generated so far)
    pub fn get_index(&self) -> u64 {
        self.vdcgen.get_index()
    }

    /// Resets the state of the sequence generator to a specific seed value
    ///
    /// # Arguments
    ///
    /// * `seed` - The seed value that determines the starting point of the sequence generation
    pub fn reseed(&mut self, seed: u64) {
        self.cirgen.reseed(seed);
        self.vdcgen.reseed(seed);
    }
}

impl Iterator for Sphere {
    type Item = [f64; 3];

    /// Returns the next point on the unit sphere
    ///
    /// This allows Sphere to be used with iterator methods like `.take()`, `.collect()`, etc.
    fn next(&mut self) -> Option<Self::Item> {
        Some(self.pop())
    }
}

impl Clone for Sphere {
    fn clone(&self) -> Self {
        Self {
            vdcgen: self.vdcgen.clone(),
            cirgen: self.cirgen.clone(),
        }
    }
}

/// Sphere-3 sequence generator using Hopf coordinates
///
/// Based on the paper:
/// Yershova, Anna, et al. "Generating uniform incremental grids on SO (3) using the Hopf fibration."
/// The International journal of robotics research 29.7 (2010): 801-812.
///
/// The Hopf fibration parametrises the 3-sphere $$S^3$$ using three VdC values:
/// $$\eta = \sqrt{v}$$, $$\psi = 2\pi v_\psi$$, $$\phi = 2\pi v_\phi$$,
/// producing $$(\cos\eta\cdot e^{i\psi},\; \sin\eta\cdot e^{i(\phi+\psi)})$$.
///
/// # Examples
///
/// ```
/// use lds_rs::Sphere3Hopf;
/// let mut sp3hgen = Sphere3Hopf::new([2, 3, 5]);
/// sp3hgen.reseed(0);
/// let res = sp3hgen.pop();
/// // Should be on the 3-sphere
/// let radius_sq = res[0] * res[0] + res[1] * res[1] + res[2] * res[2] + res[3] * res[3];
/// assert!((radius_sq - 1.0).abs() < 1e-10);
/// ```
#[cfg_attr(feature = "doc-images", doc = svgbobdoc::transform!(
/// ```svgbob
///  .───────────.  η = √φ₅(n)
///  │ φ₅(n)     │───► cos η, sin η (Hopf angle)
///  '───────────'
///  .───────────.  ψ = 2π·φ₂(n)
///  │ φ₂(n)     │───► ψ (Hopf phase)
///  '───────────'
///  .───────────.  φ = 2π·φ₃(n)
///  │ φ₃(n)     │───► φ (rotation)
///  '───────────'
///
///      combine via Hopf fibration:
///
///  S³ = (cosη·e^iψ,  sinη·e^i(φ+ψ))
///               │
///               ▼
///      .──────────────────.
///      │ 4D point on      │
///      │ unit 3-sphere S³ │
///      '──────────────────'
/// ```
))]
pub struct Sphere3Hopf {
    vdc0: VdCorput,
    vdc1: VdCorput,
    vdc2: VdCorput,
}

impl Sphere3Hopf {
    /// Creates a new Sphere3Hopf sequence generator with the given bases
    ///
    /// # Arguments
    ///
    /// * `base` - An array of three integers used as bases for generating the sequence
    pub fn new(base: [u64; 3]) -> Self {
        Self {
            vdc0: VdCorput::new(base[0]),
            vdc1: VdCorput::new(base[1]),
            vdc2: VdCorput::new(base[2]),
        }
    }

    /// Generates the next point on the 3-sphere using Hopf fibration
    ///
    /// The 3-sphere $$ S^3 $$ is parameterised by the Hopf fibration:
    ///
    /// $$ (\cos\eta\cdot e^{i\psi},\; \sin\eta\cdot e^{i(\phi+\psi)}) $$
    ///
    /// where $$ \phi,\psi \in [0, 2\pi) $$ and $$ \eta = \sqrt{v} $$ with
    /// $$ v \in \[0,1\] $$ (stratified sampling for uniform measure on $$ S^3 $$).
    ///
    /// Returns the next point on the 3-sphere as a `[f64; 4]`.
    pub fn pop(&mut self) -> [f64; 4] {
        let phi = self.vdc0.pop() * TWO_PI; // map to [0, 2π]
        let psy = self.vdc1.pop() * TWO_PI; // map to [0, 2π]
        let vdc = self.vdc2.pop();
        let cos_eta = vdc.sqrt();
        let sin_eta = (1.0 - vdc).sqrt();
        [
            cos_eta * psy.cos(),
            cos_eta * psy.sin(),
            sin_eta * (phi + psy).cos(),
            sin_eta * (phi + psy).sin(),
        ]
    }

    /// Returns the next point without advancing the state (peek)
    ///
    /// $$ (\cos\eta\cdot e^{i\psi},\; \sin\eta\cdot e^{i(\phi+\psi)}) $$
    pub fn peek(&self) -> [f64; 4] {
        let phi = self.vdc0.peek() * TWO_PI;
        let psy = self.vdc1.peek() * TWO_PI;
        let vdc = self.vdc2.peek();
        let cos_eta = vdc.sqrt();
        let sin_eta = (1.0 - vdc).sqrt();
        [
            cos_eta * psy.cos(),
            cos_eta * psy.sin(),
            sin_eta * (phi + psy).cos(),
            sin_eta * (phi + psy).sin(),
        ]
    }

    /// Skips `n` points in the sequence without computing them
    ///
    /// # Arguments
    ///
    /// * `n` - The number of points to skip
    pub fn advance(&self, n: u64) {
        self.vdc0.advance(n);
        self.vdc1.advance(n);
        self.vdc2.advance(n);
    }

    /// Returns the current index (number of points generated so far)
    pub fn get_index(&self) -> u64 {
        self.vdc0.get_index()
    }

    /// Resets the state of the sequence generator to a specific seed value
    ///
    /// # Arguments
    ///
    /// * `seed` - The seed value that determines the starting point of the sequence generation
    pub fn reseed(&mut self, seed: u64) {
        self.vdc0.reseed(seed);
        self.vdc1.reseed(seed);
        self.vdc2.reseed(seed);
    }
}

impl Iterator for Sphere3Hopf {
    type Item = [f64; 4];

    /// Returns the next point on the 3-sphere
    ///
    /// This allows Sphere3Hopf to be used with iterator methods like `.take()`, `.collect()`, etc.
    fn next(&mut self) -> Option<Self::Item> {
        Some(self.pop())
    }
}

impl Clone for Sphere3Hopf {
    fn clone(&self) -> Self {
        Self {
            vdc0: self.vdc0.clone(),
            vdc1: self.vdc1.clone(),
            vdc2: self.vdc2.clone(),
        }
    }
}

/// N-dimensional Halton sequence generator
///
/// Generates points in N-dimensional space using the Halton sequence.
///
/// # Examples
///
/// ```
/// use lds_rs::HaltonN;
/// let mut hgen = HaltonN::new(&[2, 3, 5]);
/// hgen.reseed(0);
/// let res = hgen.pop();
/// assert_eq!(res[0], 0.5);
/// assert!((res[1] - 1.0/3.0).abs() < 1e-10);
/// assert!((res[2] - 0.2).abs() < 1e-10);
/// ```
pub struct HaltonN {
    vdcs: Vec<VdCorput>,
}

impl HaltonN {
    /// Creates a new N-dimensional Halton sequence generator with the given bases
    ///
    /// # Arguments
    ///
    /// * `base` - A slice of integers used as bases for each dimension
    pub fn new(base: &[u64]) -> Self {
        let vdcs = base.iter().map(|&b| VdCorput::new(b)).collect();
        Self { vdcs }
    }

    /// Generates the next point in the N-dimensional Halton sequence
    ///
    /// Returns the next point as a `Vec<f64>`.
    pub fn pop(&mut self) -> Vec<f64> {
        let mut result = Vec::with_capacity(self.vdcs.len());
        for vdc in &mut self.vdcs {
            result.push(vdc.pop());
        }
        result
    }

    /// Resets the state of the sequence generator to a specific seed value
    ///
    /// # Arguments
    ///
    /// * `seed` - The seed value that determines the starting point of the sequence generation
    pub fn reseed(&mut self, seed: u64) {
        for vdc in &mut self.vdcs {
            vdc.reseed(seed);
        }
    }
}

impl Iterator for HaltonN {
    type Item = Vec<f64>;

    /// Returns the next point in the N-dimensional Halton sequence
    ///
    /// This allows HaltonN to be used with iterator methods like `.take()`, `.collect()`, etc.
    fn next(&mut self) -> Option<Self::Item> {
        Some(self.pop())
    }
}

impl Clone for HaltonN {
    fn clone(&self) -> Self {
        Self {
            vdcs: self.vdcs.clone(),
        }
    }
}

/// First 1000 prime numbers
///
/// Can be used as bases for low-discrepancy sequences.
pub const PRIME_TABLE: [u64; 1000] = [
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97,
    101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193,
    197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307,
    311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421,
    431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547,
    557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659,
    661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797,
    809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929,
    937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033, 1039,
    1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153,
    1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279,
    1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399, 1409,
    1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499,
    1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609, 1613,
    1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741,
    1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873,
    1877, 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999,
    2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113,
    2129, 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251,
    2267, 2269, 2273, 2281, 2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, 2371,
    2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, 2437, 2441, 2447, 2459, 2467, 2473, 2477,
    2503, 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633, 2647,
    2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731,
    2741, 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851, 2857,
    2861, 2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, 3001,
    3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137, 3163,
    3167, 3169, 3181, 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257, 3259, 3271, 3299,
    3301, 3307, 3313, 3319, 3323, 3329, 3331, 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407,
    3413, 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533, 3539,
    3541, 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, 3659,
    3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779, 3793,
    3797, 3803, 3821, 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917, 3919,
    3923, 3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051,
    4057, 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, 4153, 4157, 4159, 4177, 4201,
    4211, 4217, 4219, 4229, 4231, 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, 4327,
    4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, 4421, 4423, 4441, 4447, 4451, 4457, 4463,
    4481, 4483, 4493, 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597, 4603,
    4621, 4637, 4639, 4643, 4649, 4651, 4657, 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733,
    4751, 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889, 4903,
    4909, 4919, 4931, 4933, 4937, 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, 5009,
    5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, 5099, 5101, 5107, 5113, 5119, 5147, 5153,
    5167, 5171, 5179, 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279, 5281, 5297, 5303,
    5309, 5323, 5333, 5347, 5351, 5381, 5387, 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441,
    5443, 5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, 5527, 5531, 5557, 5563, 5569,
    5573, 5581, 5591, 5623, 5639, 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, 5701,
    5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791, 5801, 5807, 5813, 5821, 5827, 5839, 5843,
    5849, 5851, 5857, 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939, 5953, 5981, 5987,
    6007, 6011, 6029, 6037, 6043, 6047, 6053, 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131,
    6133, 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, 6229, 6247, 6257, 6263, 6269,
    6271, 6277, 6287, 6299, 6301, 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, 6373,
    6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, 6481, 6491, 6521, 6529, 6547, 6551, 6553,
    6563, 6569, 6571, 6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673, 6679, 6689, 6691,
    6701, 6703, 6709, 6719, 6733, 6737, 6761, 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829,
    6833, 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, 6947, 6949, 6959, 6961, 6967,
    6971, 6977, 6983, 6991, 6997, 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103, 7109,
    7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, 7211, 7213, 7219, 7229, 7237, 7243, 7247,
    7253, 7283, 7297, 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, 7417, 7433, 7451,
    7457, 7459, 7477, 7481, 7487, 7489, 7499, 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559,
    7561, 7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, 7649, 7669, 7673, 7681, 7687,
    7691, 7699, 7703, 7717, 7723, 7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829, 7841,
    7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919,
];

/// Integer low-discrepancy sequence generators
pub mod ilds;

/// N-dimensional sphere sequence generators
pub mod sphere_n;

#[cfg(test)]
mod tests {
    use super::*;
    use approx::assert_relative_eq;

    #[test]
    fn test_vdc_function() {
        assert_eq!(vdc(11, 2), 0.8125);
        assert_eq!(vdc(0, 2), 0.0);
        assert_eq!(vdc(1, 2), 0.5);
        assert_eq!(vdc(2, 2), 0.25);
        assert_eq!(vdc(3, 2), 0.75);
    }

    #[test]
    fn test_vdcorput_pop() {
        let mut vgen = VdCorput::new(2);
        vgen.reseed(0);
        assert_eq!(vgen.pop(), 0.5);
        assert_eq!(vgen.pop(), 0.25);
        assert_eq!(vgen.pop(), 0.75);
        assert_eq!(vgen.pop(), 0.125);
    }

    #[test]
    fn test_vdcorput_reseed() {
        let mut vgen = VdCorput::new(2);
        vgen.reseed(5);
        assert_eq!(vgen.pop(), 0.375);
        vgen.reseed(0);
        assert_eq!(vgen.pop(), 0.5);
    }

    #[test]
    fn test_vdcorput_default() {
        let mut vgen = VdCorput::default();
        vgen.reseed(0);
        assert_eq!(vgen.pop(), 0.5);
        assert_eq!(vgen.pop(), 0.25);
    }

    #[test]
    fn test_halton_pop() {
        let mut hgen = Halton::new([2, 3]);
        hgen.reseed(0);
        let res = hgen.pop();
        assert_eq!(res[0], 0.5);
        assert_relative_eq!(res[1], 1.0 / 3.0, epsilon = 1e-10);

        let res = hgen.pop();
        assert_eq!(res[0], 0.25);
        assert_relative_eq!(res[1], 2.0 / 3.0, epsilon = 1e-10);
    }

    #[test]
    fn test_circle_pop() {
        let mut cgen = Circle::new(2);
        cgen.reseed(0);
        let res = cgen.pop();
        // First point should be at angle π (180 degrees)
        assert_relative_eq!(res[0], -1.0, epsilon = 1e-10);
        assert_relative_eq!(res[1], 0.0, epsilon = 1e-10);

        let res = cgen.pop();
        // Second point should be at angle π/2 (90 degrees)
        assert_relative_eq!(res[0], 0.0, epsilon = 1e-10);
        assert_relative_eq!(res[1], 1.0, epsilon = 1e-10);
    }

    #[test]
    fn test_disk_pop() {
        let mut dgen = Disk::new([2, 3]);
        dgen.reseed(0);
        let res = dgen.pop();

        // Check that point is within unit disk
        let radius_sq = res[0] * res[0] + res[1] * res[1];
        assert!(radius_sq <= 1.0);

        // Generate a few more points and check they're all within unit disk
        for _ in 0..10 {
            let res = dgen.pop();
            let radius_sq = res[0] * res[0] + res[1] * res[1];
            assert!(radius_sq <= 1.0);
        }
    }

    #[test]
    fn test_sphere_pop() {
        let mut sgen = Sphere::new([2, 3]);
        sgen.reseed(0);
        let res = sgen.pop();

        // Check that point is on unit sphere
        let radius_sq = res[0] * res[0] + res[1] * res[1] + res[2] * res[2];
        assert_relative_eq!(radius_sq, 1.0, epsilon = 1e-10);

        // Generate a few more points and check they're all on unit sphere
        for _ in 0..5 {
            let res = sgen.pop();
            let radius_sq = res[0] * res[0] + res[1] * res[1] + res[2] * res[2];
            assert_relative_eq!(radius_sq, 1.0, epsilon = 1e-10);
        }
    }

    #[test]
    fn test_sphere3hopf_pop() {
        let mut sp3hgen = Sphere3Hopf::new([2, 3, 5]);
        sp3hgen.reseed(0);
        let res = sp3hgen.pop();

        // Check that point is on 3-sphere
        let radius_sq = res[0] * res[0] + res[1] * res[1] + res[2] * res[2] + res[3] * res[3];
        assert_relative_eq!(radius_sq, 1.0, epsilon = 1e-10);

        // Generate a few more points and check they're all on 3-sphere
        for _ in 0..5 {
            let res = sp3hgen.pop();
            let radius_sq = res[0] * res[0] + res[1] * res[1] + res[2] * res[2] + res[3] * res[3];
            assert_relative_eq!(radius_sq, 1.0, epsilon = 1e-10);
        }
    }

    #[test]
    fn test_haltonn_pop() {
        let mut hgen = HaltonN::new(&[2, 3, 5]);
        hgen.reseed(0);
        let res = hgen.pop();

        assert_eq!(res[0], 0.5);
        assert_relative_eq!(res[1], 1.0 / 3.0, epsilon = 1e-10);
        assert_relative_eq!(res[2], 0.2, epsilon = 1e-10);

        let res = hgen.pop();
        assert_eq!(res[0], 0.25);
        assert_relative_eq!(res[1], 2.0 / 3.0, epsilon = 1e-10);
        assert_relative_eq!(res[2], 0.4, epsilon = 1e-10);
    }

    #[test]
    fn test_prime_table() {
        // Check first few primes
        assert_eq!(PRIME_TABLE[0], 2);
        assert_eq!(PRIME_TABLE[1], 3);
        assert_eq!(PRIME_TABLE[2], 5);
        assert_eq!(PRIME_TABLE[3], 7);
        assert_eq!(PRIME_TABLE[4], 11);

        // Check table length
        assert_eq!(PRIME_TABLE.len(), 1000);

        // Check a prime near the end
        assert_eq!(PRIME_TABLE[999], 7919);
    }

    // Additional comprehensive tests for edge cases and different bases

    #[test]
    fn test_vdc_function_edge_cases() {
        // Test with count=0
        assert_eq!(vdc(0, 2), 0.0);
        assert_eq!(vdc(0, 3), 0.0);
        assert_eq!(vdc(0, 5), 0.0);

        // Test with count=1
        assert_eq!(vdc(1, 2), 0.5);
        assert_eq!(vdc(1, 3), 1.0 / 3.0);
        assert_eq!(vdc(1, 5), 0.2);

        // Test with different bases
        assert_eq!(vdc(10, 2), 0.3125); // 1010 in binary -> 0.0101
        assert_eq!(vdc(10, 3), 0.37037037037037035); // 101 in base 3 -> 0.101
        assert_eq!(vdc(10, 5), 0.08); // 20 in base 5 -> 0.02

        // Test with larger numbers
        let result = vdc(1000, 2);
        assert!((0.0..1.0).contains(&result));

        // Test with prime bases
        assert_eq!(vdc(6, 7), 6.0 / 7.0);
        assert_eq!(vdc(12, 11), 0.09917355371900827);
    }

    #[test]
    fn test_vdcorput_different_bases() {
        // Test with base 3
        let mut vgen = VdCorput::new(3);
        vgen.reseed(0);
        assert_relative_eq!(vgen.pop(), 1.0 / 3.0, epsilon = 1e-10);
        assert_relative_eq!(vgen.pop(), 2.0 / 3.0, epsilon = 1e-10);
        assert_relative_eq!(vgen.pop(), 1.0 / 9.0, epsilon = 1e-10);

        // Test with base 5
        let mut vgen = VdCorput::new(5);
        vgen.reseed(0);
        assert_relative_eq!(vgen.pop(), 0.2, epsilon = 1e-10);
        assert_relative_eq!(vgen.pop(), 0.4, epsilon = 1e-10);
        assert_relative_eq!(vgen.pop(), 0.6, epsilon = 1e-10);
        assert_relative_eq!(vgen.pop(), 0.8, epsilon = 1e-10);
        assert_relative_eq!(vgen.pop(), 0.04, epsilon = 1e-10);

        // Test with base 7
        let mut vgen = VdCorput::new(7);
        vgen.reseed(0);
        assert_relative_eq!(vgen.pop(), 1.0 / 7.0, epsilon = 1e-10);
        assert_relative_eq!(vgen.pop(), 2.0 / 7.0, epsilon = 1e-10);
    }

    #[test]
    fn test_vdcorput_large_values() {
        let mut vgen = VdCorput::new(2);
        vgen.reseed(1000);

        // Generate several values and ensure they're valid
        for _ in 0..10 {
            let value = vgen.pop();
            assert!((0.0..1.0).contains(&value));
        }
    }

    #[test]
    fn test_halton_different_bases() {
        // Test with bases 3 and 5
        let mut hgen = Halton::new([3, 5]);
        hgen.reseed(0);
        let res = hgen.pop();
        assert_relative_eq!(res[0], 1.0 / 3.0, epsilon = 1e-10);
        assert_relative_eq!(res[1], 0.2, epsilon = 1e-10);

        // Test with bases 5 and 7
        let mut hgen = Halton::new([5, 7]);
        hgen.reseed(0);
        let res = hgen.pop();
        assert_relative_eq!(res[0], 0.2, epsilon = 1e-10);
        assert_relative_eq!(res[1], 1.0 / 7.0, epsilon = 1e-10);
    }

    #[test]
    fn test_circle_different_bases() {
        // Test with base 3
        let mut cgen = Circle::new(3);
        cgen.reseed(0);
        let res = cgen.pop();
        // Should be on unit circle
        let radius_sq = res[0] * res[0] + res[1] * res[1];
        assert_relative_eq!(radius_sq, 1.0, epsilon = 1e-10);

        // Test with base 5
        let mut cgen = Circle::new(5);
        cgen.reseed(0);
        let res = cgen.pop();
        let radius_sq = res[0] * res[0] + res[1] * res[1];
        assert_relative_eq!(radius_sq, 1.0, epsilon = 1e-10);
    }

    #[test]
    fn test_disk_different_bases() {
        // Test with bases 3 and 5
        let mut dgen = Disk::new([3, 5]);
        dgen.reseed(0);
        let res = dgen.pop();
        let radius_sq = res[0] * res[0] + res[1] * res[1];
        assert!(radius_sq <= 1.0);

        // Test with bases 5 and 7
        let mut dgen = Disk::new([5, 7]);
        dgen.reseed(0);
        let res = dgen.pop();
        let radius_sq = res[0] * res[0] + res[1] * res[1];
        assert!(radius_sq <= 1.0);
    }

    #[test]
    fn test_sphere_different_bases() {
        // Test with bases 3 and 5
        let mut sgen = Sphere::new([3, 5]);
        sgen.reseed(0);
        let res = sgen.pop();
        let radius_sq = res[0] * res[0] + res[1] * res[1] + res[2] * res[2];
        assert_relative_eq!(radius_sq, 1.0, epsilon = 1e-10);

        // Test with bases 5 and 7
        let mut sgen = Sphere::new([5, 7]);
        sgen.reseed(0);
        let res = sgen.pop();
        let radius_sq = res[0] * res[0] + res[1] * res[1] + res[2] * res[2];
        assert_relative_eq!(radius_sq, 1.0, epsilon = 1e-10);
    }

    #[test]
    fn test_sphere3hopf_different_bases() {
        // Test with bases 3, 5, 7
        let mut sgen = Sphere3Hopf::new([3, 5, 7]);
        sgen.reseed(0);
        let res = sgen.pop();
        let radius_sq = res[0] * res[0] + res[1] * res[1] + res[2] * res[2] + res[3] * res[3];
        assert_relative_eq!(radius_sq, 1.0, epsilon = 1e-10);

        // Test with bases 5, 7, 11
        let mut sgen = Sphere3Hopf::new([5, 7, 11]);
        sgen.reseed(0);
        let res = sgen.pop();
        let radius_sq = res[0] * res[0] + res[1] * res[1] + res[2] * res[2] + res[3] * res[3];
        assert_relative_eq!(radius_sq, 1.0, epsilon = 1e-10);
    }

    #[test]
    fn test_haltonn_different_bases() {
        // Test with 4 dimensions
        let mut hgen = HaltonN::new(&[3, 5, 7, 11]);
        hgen.reseed(0);
        let res = hgen.pop();
        assert_eq!(res.len(), 4);
        assert_relative_eq!(res[0], 1.0 / 3.0, epsilon = 1e-10);
        assert_relative_eq!(res[1], 0.2, epsilon = 1e-10);
        assert_relative_eq!(res[2], 1.0 / 7.0, epsilon = 1e-10);
        assert_relative_eq!(res[3], 1.0 / 11.0, epsilon = 1e-10);

        // Test with 5 dimensions
        let mut hgen = HaltonN::new(&[2, 3, 5, 7, 11]);
        hgen.reseed(0);
        let res = hgen.pop();
        assert_eq!(res.len(), 5);
        assert_relative_eq!(res[0], 0.5, epsilon = 1e-10);
        assert_relative_eq!(res[1], 1.0 / 3.0, epsilon = 1e-10);
        assert_relative_eq!(res[2], 0.2, epsilon = 1e-10);
        assert_relative_eq!(res[3], 1.0 / 7.0, epsilon = 1e-10);
        assert_relative_eq!(res[4], 1.0 / 11.0, epsilon = 1e-10);
    }

    #[test]
    fn test_default_implementations() {
        // Test Default for VdCorput
        let mut vgen = VdCorput::default();
        vgen.reseed(0);
        assert_eq!(vgen.pop(), 0.5);

        // Test that default base is 2
        let vgen_default = VdCorput::default();
        let vgen_explicit = VdCorput::new(2);
        assert_eq!(vgen_default.base, vgen_explicit.base);
    }

    #[test]
    fn test_sequence_properties() {
        // Test that VdCorput sequence values are always in [0, 1)
        let mut vgen = VdCorput::new(2);
        for _ in 0..100 {
            let value = vgen.pop();
            assert!((0.0..1.0).contains(&value));
        }

        // Test that Halton sequence values are always in [0, 1) for each dimension
        let mut hgen = Halton::new([2, 3]);
        for _ in 0..100 {
            let res = hgen.pop();
            assert!(res[0] >= 0.0 && res[0] < 1.0);
            assert!(res[1] >= 0.0 && res[1] < 1.0);
        }

        // Test that Circle sequence points are always on unit circle
        let mut cgen = Circle::new(2);
        for _ in 0..100 {
            let res = cgen.pop();
            let radius_sq = res[0] * res[0] + res[1] * res[1];
            assert_relative_eq!(radius_sq, 1.0, epsilon = 1e-10);
        }

        // Test that Disk sequence points are always within unit disk
        let mut dgen = Disk::new([2, 3]);
        for _ in 0..100 {
            let res = dgen.pop();
            let radius_sq = res[0] * res[0] + res[1] * res[1];
            assert!(radius_sq <= 1.0);
        }

        // Test that Sphere sequence points are always on unit sphere
        let mut sgen = Sphere::new([2, 3]);
        for _ in 0..100 {
            let res = sgen.pop();
            let radius_sq = res[0] * res[0] + res[1] * res[1] + res[2] * res[2];
            assert_relative_eq!(radius_sq, 1.0, epsilon = 1e-10);
        }

        // Test that Sphere3Hopf sequence points are always on 3-sphere
        let mut sgen = Sphere3Hopf::new([2, 3, 5]);
        for _ in 0..100 {
            let res = sgen.pop();
            let radius_sq = res.iter().map(|&x| x * x).sum::<f64>();
            assert_relative_eq!(radius_sq, 1.0, epsilon = 1e-10);
        }
    }

    #[test]
    fn test_reseed_consistency() {
        // Test that reseed with the same value produces the same sequence
        let mut vgen = VdCorput::new(2);

        vgen.reseed(10);
        let seq1: Vec<_> = (0..5).map(|_| vgen.pop()).collect();

        vgen.reseed(10);
        let seq2: Vec<_> = (0..5).map(|_| vgen.pop()).collect();

        for i in 0..5 {
            assert_relative_eq!(seq1[i], seq2[i], epsilon = 1e-10);
        }

        // Test that reseed with different values produces different sequences
        vgen.reseed(10);
        let seq3: Vec<_> = (0..5).map(|_| vgen.pop()).collect();

        vgen.reseed(20);
        let seq4: Vec<_> = (0..5).map(|_| vgen.pop()).collect();

        let mut different = false;
        for i in 0..5 {
            if (seq3[i] - seq4[i]).abs() > 1e-10 {
                different = true;
                break;
            }
        }
        assert!(
            different,
            "Sequences with different seeds should be different"
        );
    }

    #[test]
    fn test_vdcorput_peek() {
        let mut vgen = VdCorput::new(2);
        vgen.reseed(0);
        let peeked = vgen.peek();
        assert_relative_eq!(peeked, 0.5, epsilon = 1e-10);
        let popped = vgen.pop();
        assert_relative_eq!(popped, 0.5, epsilon = 1e-10); // peek doesn't advance
        assert_relative_eq!(vgen.peek(), 0.25, epsilon = 1e-10);
    }

    #[test]
    fn test_vdcorput_advance() {
        let mut vgen = VdCorput::new(2);
        vgen.reseed(0);
        vgen.advance(3);
        assert_relative_eq!(vgen.pop(), 0.125, epsilon = 1e-10);
        vgen.reseed(0);
        vgen.advance(4);
        // vdc(5, 2) = binary 101 reversed = 0.101 = 0.625
        assert_relative_eq!(vgen.pop(), 0.625, epsilon = 1e-10);
    }

    #[test]
    fn test_vdcorput_get_index() {
        let mut vgen = VdCorput::new(2);
        assert_eq!(vgen.get_index(), 0);
        vgen.pop();
        assert_eq!(vgen.get_index(), 1);
        vgen.pop();
        assert_eq!(vgen.get_index(), 2);
        vgen.reseed(5);
        assert_eq!(vgen.get_index(), 5);
    }

    #[test]
    fn test_halton_peek() {
        let mut hgen = Halton::new([2, 3]);
        hgen.reseed(0);
        let peeked = hgen.peek();
        assert_relative_eq!(peeked[0], 0.5, epsilon = 1e-10);
        assert_relative_eq!(peeked[1], 1.0 / 3.0, epsilon = 1e-10);
        let popped = hgen.pop();
        assert_relative_eq!(popped[0], 0.5, epsilon = 1e-10);
        assert_relative_eq!(popped[1], 1.0 / 3.0, epsilon = 1e-10);
    }

    #[test]
    fn test_halton_advance() {
        let mut hgen = Halton::new([2, 3]);
        hgen.reseed(0);
        hgen.advance(2);
        let popped = hgen.pop();
        assert_relative_eq!(popped[0], 0.75, epsilon = 1e-10);
        assert_relative_eq!(popped[1], 1.0 / 9.0, epsilon = 1e-10);
    }

    #[test]
    fn test_halton_iterator() {
        let mut hgen = Halton::new([2, 3]);
        hgen.reseed(0);
        let values: Vec<[f64; 2]> = hgen.take(3).collect();
        assert_eq!(values.len(), 3);
        assert_relative_eq!(values[0][0], 0.5, epsilon = 1e-10);
        assert_relative_eq!(values[0][1], 1.0 / 3.0, epsilon = 1e-10);
        assert_relative_eq!(values[1][0], 0.25, epsilon = 1e-10);
        assert_relative_eq!(values[1][1], 2.0 / 3.0, epsilon = 1e-10);
    }

    #[test]
    fn test_circle_peek() {
        let mut cgen = Circle::new(2);
        cgen.reseed(0);
        let peeked = cgen.peek();
        let popped = cgen.pop();
        assert_relative_eq!(peeked[0], popped[0], epsilon = 1e-10);
        assert_relative_eq!(peeked[1], popped[1], epsilon = 1e-10);
    }

    #[test]
    fn test_circle_advance() {
        let mut cgen = Circle::new(2);
        cgen.reseed(0);
        cgen.advance(5);
        let res = cgen.pop();
        let radius_sq = res[0] * res[0] + res[1] * res[1];
        assert_relative_eq!(radius_sq, 1.0, epsilon = 1e-10);
    }

    #[test]
    fn test_circle_get_index() {
        let mut cgen = Circle::new(2);
        assert_eq!(cgen.get_index(), 0);
        cgen.pop();
        assert_eq!(cgen.get_index(), 1);
    }

    #[test]
    fn test_disk_peek() {
        let mut dgen = Disk::new([2, 3]);
        dgen.reseed(0);
        let peeked = dgen.peek();
        let popped = dgen.pop();
        assert_relative_eq!(peeked[0], popped[0], epsilon = 1e-10);
        assert_relative_eq!(peeked[1], popped[1], epsilon = 1e-10);
    }

    #[test]
    fn test_disk_advance() {
        let mut dgen = Disk::new([2, 3]);
        dgen.reseed(0);
        dgen.advance(5);
        let res = dgen.pop();
        let radius_sq = res[0] * res[0] + res[1] * res[1];
        assert!(radius_sq <= 1.0);
    }

    #[test]
    fn test_sphere_peek() {
        let mut sgen = Sphere::new([2, 3]);
        sgen.reseed(0);
        let peeked = sgen.peek();
        let popped = sgen.pop();
        assert_relative_eq!(peeked[0], popped[0], epsilon = 1e-10);
        assert_relative_eq!(peeked[1], popped[1], epsilon = 1e-10);
        assert_relative_eq!(peeked[2], popped[2], epsilon = 1e-10);
    }

    #[test]
    fn test_sphere_advance() {
        let mut sgen = Sphere::new([2, 3]);
        sgen.reseed(0);
        sgen.advance(5);
        let res = sgen.pop();
        let radius_sq = res[0] * res[0] + res[1] * res[1] + res[2] * res[2];
        assert_relative_eq!(radius_sq, 1.0, epsilon = 1e-10);
    }

    #[test]
    fn test_sphere3hopf_peek() {
        let mut sgen = Sphere3Hopf::new([2, 3, 5]);
        sgen.reseed(0);
        let peeked = sgen.peek();
        let popped = sgen.pop();
        for i in 0..4 {
            assert_relative_eq!(peeked[i], popped[i], epsilon = 1e-10);
        }
    }

    #[test]
    fn test_sphere3hopf_advance() {
        let mut sgen = Sphere3Hopf::new([2, 3, 5]);
        sgen.reseed(0);
        sgen.advance(5);
        let res = sgen.pop();
        let radius_sq = res.iter().map(|&x| x * x).sum::<f64>();
        assert_relative_eq!(radius_sq, 1.0, epsilon = 1e-10);
    }

    #[test]
    fn test_sphere3hopf_get_index() {
        let mut sgen = Sphere3Hopf::new([2, 3, 5]);
        assert_eq!(sgen.get_index(), 0);
        sgen.pop();
        assert_eq!(sgen.get_index(), 1);
    }

    #[test]
    fn test_halton_get_index() {
        let mut hgen = Halton::new([2, 3]);
        assert_eq!(hgen.get_index(), 0);
        hgen.pop();
        assert_eq!(hgen.get_index(), 1);
    }

    #[test]
    fn test_circle_clone() {
        let mut cgen = Circle::new(2);
        cgen.reseed(5);
        let mut cloned = cgen.clone();
        assert_eq!(cloned.pop(), cgen.pop());
    }

    #[test]
    fn test_disk_clone() {
        let mut dgen = Disk::new([2, 3]);
        dgen.reseed(5);
        let mut cloned = dgen.clone();
        assert_eq!(cloned.pop(), dgen.pop());
    }

    #[test]
    fn test_sphere_clone() {
        let mut sgen = Sphere::new([2, 3]);
        sgen.reseed(5);
        let mut cloned = sgen.clone();
        assert_eq!(cloned.pop(), sgen.pop());
    }

    #[test]
    fn test_sphere3hopf_clone() {
        let mut sgen = Sphere3Hopf::new([2, 3, 5]);
        sgen.reseed(5);
        let mut cloned = sgen.clone();
        assert_eq!(cloned.pop(), sgen.pop());
    }

    #[test]
    fn test_haltonn_clone() {
        let mut hgen = HaltonN::new(&[2, 3, 5]);
        hgen.reseed(5);
        let mut cloned = hgen.clone();
        assert_eq!(cloned.pop(), hgen.pop());
    }
}