dreamwell_math/

lib.rs

//! # dreamwell-math — Classical Mathematics for Simulation
//!
//! Pure mathematical primitives for the Dreamwell engine and research ecosystem.
//! No engine dependencies. No allocation on hot path. Standalone usable.
//!
//! ## Modules
//! - **Complex**: Complex number arithmetic (add, mul, conj, exp_i, norm)
//! - **Spatial**: Distance metrics, Gaussian kernels, AABB extraction
//! - **Wave**: Traveling wave, standing wave, wave packets, dispersion
//! - **Fibonacci**: Golden ratio, phyllotaxis disc, Fibonacci sphere
//! - **Interpolation**: Lerp, exponential decay, smoothstep
//! - **Random**: Deterministic splitmix64 PRNG
//!
//! - **Linear Algebra**: Complex matrix multiply (CGEMM), tiled blocking
//! - **Eigenvalue**: Jacobi rotation solver for Hermitian matrices
//! - **Matrix Exponential**: Scaling-and-squaring with Padé approximant
//!
//! ## Clean Compute Compliance
//! - Zero heap allocation in caller-provides-output functions
//! - All functions are `#[inline]` where beneficial
//! - No external dependencies
//! - Deterministic: same inputs always produce same outputs

// ══════════════════════════════════════════════════════════════
// Submodules — Linear Algebra, Eigenvalues, Matrix Exponential
// ══════════════════════════════════════════════════════════════

pub mod eigen;
pub mod linalg;
pub mod matrix_exp;

#[cfg(feature = "gpu-compute")]
pub mod gpu_linalg;

#[cfg(feature = "gpu-compute")]
pub mod gpu_training;

// ══════════════════════════════════════════════════════════════
// Complex Numbers
// ══════════════════════════════════════════════════════════════

/// Complex number z = re + i·im.
///
/// Used as the base type for density matrices, Hamiltonians, and
/// interference kernels in the quantum simulation layer.
/// Arithmetic follows standard complex field rules.
#[derive(Debug, Clone, Copy, Default, PartialEq)]
#[repr(C)]
pub struct Complex {
    pub re: f32,
    pub im: f32,
}

// SAFETY: Complex is #[repr(C)] with two f32 fields — trivially Pod and Zeroable.
#[cfg(feature = "gpu-compute")]
unsafe impl bytemuck::Pod for Complex {}
#[cfg(feature = "gpu-compute")]
unsafe impl bytemuck::Zeroable for Complex {}

impl Complex {
    pub const ZERO: Self = Self { re: 0.0, im: 0.0 };
    pub const ONE: Self = Self { re: 1.0, im: 0.0 };
    pub const I: Self = Self { re: 0.0, im: 1.0 };

    #[inline]
    pub fn new(re: f32, im: f32) -> Self {
        Self { re, im }
    }
    #[inline]
    pub fn from_real(r: f32) -> Self {
        Self { re: r, im: 0.0 }
    }
    #[inline]
    pub fn from_imag(i: f32) -> Self {
        Self { re: 0.0, im: i }
    }

    /// Complex conjugate: z* = re - i·im
    #[inline]
    pub fn conj(self) -> Self {
        Self {
            re: self.re,
            im: -self.im,
        }
    }

    /// Squared magnitude: |z|² = re² + im²
    #[inline]
    pub fn norm_sq(self) -> f32 {
        self.re * self.re + self.im * self.im
    }

    /// Magnitude: |z| = √(re² + im²)
    #[inline]
    pub fn norm(self) -> f32 {
        self.norm_sq().sqrt()
    }

    /// Complex multiplication: (a+bi)(c+di) = (ac-bd) + (ad+bc)i
    #[inline]
    pub fn mul(self, other: Self) -> Self {
        Self {
            re: self.re * other.re - self.im * other.im,
            im: self.re * other.im + self.im * other.re,
        }
    }

    #[inline]
    pub fn add(self, other: Self) -> Self {
        Self {
            re: self.re + other.re,
            im: self.im + other.im,
        }
    }

    #[inline]
    pub fn sub(self, other: Self) -> Self {
        Self {
            re: self.re - other.re,
            im: self.im - other.im,
        }
    }

    #[inline]
    pub fn scale(self, s: f32) -> Self {
        Self {
            re: self.re * s,
            im: self.im * s,
        }
    }

    /// Euler's formula: e^{iθ} = cos(θ) + i·sin(θ)
    #[inline]
    pub fn exp_i(theta: f32) -> Self {
        Self {
            re: theta.cos(),
            im: theta.sin(),
        }
    }

    /// Complex exponential: e^z = e^re · (cos(im) + i·sin(im))
    #[inline]
    pub fn exp(self) -> Self {
        let r = self.re.exp();
        Self {
            re: r * self.im.cos(),
            im: r * self.im.sin(),
        }
    }
}

// ══════════════════════════════════════════════════════════════
// Spatial Mathematics
// ══════════════════════════════════════════════════════════════

/// Squared Euclidean distance between two 3D points.
/// Avoids sqrt for comparison (cheaper than full distance).
#[inline]
pub fn distance_sq_3d(a: &[f32; 3], b: &[f32; 3]) -> f32 {
    let dx = a[0] - b[0];
    let dy = a[1] - b[1];
    let dz = a[2] - b[2];
    dx * dx + dy * dy + dz * dz
}

/// Euclidean distance between two 3D points.
#[inline]
pub fn distance_3d(a: &[f32; 3], b: &[f32; 3]) -> f32 {
    distance_sq_3d(a, b).sqrt()
}

/// Gaussian radial basis function: exp(-d²/σ²).
/// Returns 1.0 at center, decays to ~0 at distance 3σ.
/// Used for spatial attention and relevance scoring.
#[inline]
pub fn gaussian_kernel(dist_sq: f32, sigma_sq: f32) -> f32 {
    (-dist_sq / sigma_sq).exp()
}

/// 3D dot product.
#[inline]
pub fn vec3_dot(a: [f32; 3], b: [f32; 3]) -> f32 {
    a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
}

/// 3D vector length.
#[inline]
pub fn vec3_len(v: [f32; 3]) -> f32 {
    vec3_dot(v, v).sqrt()
}

/// 3D vector scale.
#[inline]
pub fn vec3_scale(v: [f32; 3], s: f32) -> [f32; 3] {
    [v[0] * s, v[1] * s, v[2] * s]
}

/// 3D cross product.
#[inline]
pub fn vec3_cross(a: [f32; 3], b: [f32; 3]) -> [f32; 3] {
    [
        a[1] * b[2] - a[2] * b[1],
        a[2] * b[0] - a[0] * b[2],
        a[0] * b[1] - a[1] * b[0],
    ]
}

/// 3D vector normalize (returns zero vector if length < epsilon).
#[inline]
pub fn vec3_normalize(v: [f32; 3]) -> [f32; 3] {
    let len = vec3_len(v);
    if len < 1e-8 { [0.0; 3] } else { vec3_scale(v, 1.0 / len) }
}

// ══════════════════════════════════════════════════════════════
// Wave Mathematics
// ══════════════════════════════════════════════════════════════

/// Traveling wave: y = A·sin(kx - ωt).
/// k = wave number (2π/λ), ω = angular frequency (2πf).
/// Propagates in +x direction.
///
/// Reference: d'Alembert wave equation solution (1747).
#[inline]
pub fn traveling_wave(x: f32, t: f32, k: f32, omega: f32, amplitude: f32) -> f32 {
    amplitude * (k * x - omega * t).sin()
}

/// Standing wave: y = A·sin(kx)·cos(ωt).
/// Created by superposition of two counter-propagating traveling waves.
#[inline]
pub fn standing_wave(x: f32, t: f32, k: f32, omega: f32, amplitude: f32) -> f32 {
    amplitude * (k * x).sin() * (omega * t).cos()
}

/// 2D traveling wave with perpendicular tapering:
/// y = A·sin(k·along - ωt)·cos(k·perp·taper)
///
/// Used by the DreamWaveController for fabric-like particle motion.
#[inline]
pub fn traveling_wave_2d(along: f32, perp: f32, t: f32, k: f32, omega: f32, amplitude: f32, taper: f32) -> f32 {
    amplitude * (k * along - omega * t).sin() * (k * perp * taper).cos()
}

/// Wave number from wavelength: k = 2π/λ.
#[inline]
pub fn wave_number(wavelength: f32) -> f32 {
    std::f32::consts::TAU / wavelength
}

/// Angular frequency from frequency: ω = 2πf.
#[inline]
pub fn angular_frequency(frequency: f32) -> f32 {
    std::f32::consts::TAU * frequency
}

/// Linear dispersion relation: ω = c·k (non-dispersive medium).
#[inline]
pub fn dispersion_linear(k: f32, speed: f32) -> f32 {
    speed * k
}

/// Gaussian wave packet: ψ(x,t) = exp(-(x-x0-v_g·t)²/(4σ²)) · exp(i(k0·x - ω0·t))
/// Returns the real part (amplitude envelope × carrier).
#[inline]
pub fn wave_packet_gaussian(x: f32, t: f32, k0: f32, sigma: f32, group_velocity: f32) -> f32 {
    let center = group_velocity * t;
    let envelope = (-(x - center).powi(2) / (4.0 * sigma * sigma)).exp();
    let carrier = (k0 * x - dispersion_linear(k0, group_velocity) * t).cos();
    envelope * carrier
}

// ══════════════════════════════════════════════════════════════
// Fibonacci / Golden Ratio
// ══════════════════════════════════════════════════════════════

/// The golden ratio φ = (1 + √5) / 2 ≈ 1.6180339887.
#[inline]
pub fn golden_ratio() -> f32 {
    (1.0 + 5.0f32.sqrt()) / 2.0
}

/// The golden angle θ = 2π/φ² ≈ 2.3999632297 radians ≈ 137.508°.
/// Used in phyllotaxis for optimal packing on discs and spheres.
///
/// Reference: Vogel, H. (1979). "A better model for sunflower phyllotaxis."
#[inline]
pub fn golden_angle() -> f32 {
    let phi = golden_ratio();
    std::f32::consts::TAU / (phi * phi)
}

/// Fibonacci sphere: N points uniformly distributed on a unit sphere
/// using the golden ratio spiral.
///
/// Returns Vec of [x, y, z] unit vectors.
/// Deterministic: same count always produces the same distribution.
pub fn fibonacci_sphere(count: u32) -> Vec<[f32; 3]> {
    let phi = golden_ratio();
    let angle_inc = std::f32::consts::TAU * phi;
    let mut points = Vec::with_capacity(count as usize);
    for i in 0..count {
        let t = (i as f32 + 0.5) / count as f32;
        let polar = (1.0 - 2.0 * t).acos();
        let azimuth = angle_inc * i as f32;
        points.push([polar.sin() * azimuth.cos(), polar.cos(), polar.sin() * azimuth.sin()]);
    }
    points
}

/// Phyllotaxis disc: N points on a flat disc with golden angle spacing.
/// Radius grows as √(index) for uniform area density.
///
/// `rotation`: additional rotation offset (radians), e.g. for animation.
pub fn phyllotaxis_disc(count: u32, radius: f32, rotation: f32) -> Vec<[f32; 3]> {
    let ga = golden_angle();
    let mut points = Vec::with_capacity(count as usize);
    for i in 0..count {
        let frac = (i as f32 + 0.5) / count as f32;
        let r = frac.sqrt() * radius;
        let theta = i as f32 * ga + rotation;
        points.push([r * theta.cos(), 0.0, r * theta.sin()]);
    }
    points
}

// ══════════════════════════════════════════════════════════════
// Interpolation
// ══════════════════════════════════════════════════════════════

/// Linear interpolation: a + (b - a) · t.
#[inline]
pub fn lerp(a: f32, b: f32, t: f32) -> f32 {
    a + (b - a) * t
}

/// 3D linear interpolation.
#[inline]
pub fn lerp3(a: [f32; 3], b: [f32; 3], t: f32) -> [f32; 3] {
    [lerp(a[0], b[0], t), lerp(a[1], b[1], t), lerp(a[2], b[2], t)]
}

/// Exponential decay toward target: current + (target - current) · (1 - e^{-rate·dt}).
///
/// This is the discrete-time solution to the first-order ODE dx/dt = -rate·(x - target).
/// Also known as the Lindblad decoherence kernel when applied to off-diagonal elements.
///
/// Properties:
/// - rate = 0: no change
/// - rate -> infinity: instant snap to target
/// - Typical game rates: 2.0 (gentle, ~1.5s), 8.0 (snappy, ~0.25s)
#[inline]
pub fn exp_decay(current: f32, target: f32, rate: f32, dt: f32) -> f32 {
    current + (target - current) * (1.0 - (-rate * dt).exp())
}

/// Smoothstep: Hermite interpolation between edge0 and edge1.
/// Returns 0 below edge0, 1 above edge1, smooth curve between.
#[inline]
pub fn smoothstep(edge0: f32, edge1: f32, x: f32) -> f32 {
    let t = ((x - edge0) / (edge1 - edge0)).clamp(0.0, 1.0);
    t * t * (3.0 - 2.0 * t)
}

// ══════════════════════════════════════════════════════════════
// Deterministic Random
// ══════════════════════════════════════════════════════════════

/// Splitmix64: deterministic PRNG. Returns (value, next_state).
///
/// Used throughout the engine for reproducible benchmarks, seed-based
/// procedural generation, and deterministic Born rule measurement.
///
/// Reference: Vigna, S. (2015). "An experimental exploration of Marsaglia's
/// xorshift generators, scrambled."
#[inline]
pub fn splitmix64(state: u64) -> (u64, u64) {
    let s = state.wrapping_add(0x9E3779B97F4A7C15);
    let mut z = s;
    z = (z ^ (z >> 30)).wrapping_mul(0xBF58476D1CE4E5B9);
    z = (z ^ (z >> 27)).wrapping_mul(0x94D049BB133111EB);
    z = z ^ (z >> 31);
    (z, s)
}

/// Splitmix64 -> f32 in [0, 1). Deterministic.
#[inline]
pub fn splitmix64_f32(state: u64) -> (f32, u64) {
    let (val, next) = splitmix64(state);
    ((val >> 40) as f32 / (1u64 << 24) as f32, next)
}

// ══════════════════════════════════════════════════════════════
// Atomic Operations (for parallel solvers)
// ══════════════════════════════════════════════════════════════

/// Atomic f32 addition via compare-exchange on AtomicU32.
///
/// Used by the Parallel Jacobi solver for lock-free impulse accumulation
/// across rayon thread boundaries. Zero allocation. Wait-free under low contention.
#[inline]
pub fn atomic_f32_add(atom: &std::sync::atomic::AtomicU32, val: f32) {
    use std::sync::atomic::Ordering;
    loop {
        let bits = atom.load(Ordering::Relaxed);
        let current = f32::from_bits(bits);
        let new = current + val;
        match atom.compare_exchange_weak(bits, new.to_bits(), Ordering::Relaxed, Ordering::Relaxed) {
            Ok(_) => break,
            Err(_) => {} // CAS retry on contention
        }
    }
}

// ══════════════════════════════════════════════════════════════
// Spatial Hashing & Broadphase Primitives
// ══════════════════════════════════════════════════════════════

/// Forward half-neighborhood offsets for broadphase pair generation.
/// 14 cells: self + 13 lexicographically forward neighbors.
/// Produces zero-duplicate pairs by construction (Causal Spatial Mask).
///
/// This is the spatial analog of the causal attention mask in transformers:
/// only compute the upper triangle of the spatial relevance matrix.
///
/// Reference: Allen & Tildesley (1987), half-neighbor Verlet list.
/// Dreamwell innovation: applied to hash-grid broadphase (novel for game engines).
pub const HALF_NEIGHBORHOOD: [(i32, i32, i32); 14] = [
    (0, 0, 0), // self (j > i guard)
    (0, 0, 1),
    (0, 1, -1),
    (0, 1, 0),
    (0, 1, 1),
    (1, -1, -1),
    (1, -1, 0),
    (1, -1, 1),
    (1, 0, -1),
    (1, 0, 0),
    (1, 0, 1),
    (1, 1, -1),
    (1, 1, 0),
    (1, 1, 1),
];

/// FNV-1a hash of a 3D integer cell coordinate.
/// Used for spatial hash grid broadphase. Deterministic, fast, low collision rate.
#[inline]
pub fn fnv1a_3d(ix: i32, iy: i32, iz: i32) -> u64 {
    let mut h: u64 = 0xcbf29ce484222325;
    for byte in ix
        .to_le_bytes()
        .iter()
        .chain(iy.to_le_bytes().iter())
        .chain(iz.to_le_bytes().iter())
    {
        h ^= *byte as u64;
        h = h.wrapping_mul(0x100000001b3);
    }
    h
}

/// Adiabatic force distribution: returns the fractional force to apply
/// on the current step within an adiabatic window.
///
/// Based on the Quantum Adiabatic Theorem (Born & Fock, 1928):
/// slow Hamiltonian changes keep the system near its ground state.
/// Applied to game physics: spreading impulses over N steps preserves
/// broadphase temporal coherence and eliminates p99 frame spikes.
///
/// Returns `force / window_size` if within the window, 0.0 otherwise.
#[inline]
pub fn adiabatic_fraction(step: u32, interval: u32, window: u32) -> f32 {
    let phase = step % interval;
    if step > 0 && phase < window {
        1.0 / window as f32
    } else {
        0.0
    }
}

// ══════════════════════════════════════════════════════════════
// Physics Primitives (extracted from Parallel Jacobi solver)
// ══════════════════════════════════════════════════════════════

/// Contact impulse magnitude (sequential impulse formula).
/// j = -(1 + e) * v_along_normal / inv_mass_sum
///
/// Reference: Catto, E. (2005). "Iterative Dynamics with Temporal Coherence." GDC.
#[inline]
pub fn contact_impulse(v_along_normal: f32, restitution: f32, inv_mass_sum: f32) -> f32 {
    if inv_mass_sum < 1e-8 {
        return 0.0;
    }
    -(1.0 + restitution) * v_along_normal / inv_mass_sum
}

/// Baumgarte positional correction for penetration resolution.
/// correction = max(0, depth - slop) * pct / inv_mass_sum
///
/// Prevents sinking by directly adjusting positions after impulse resolution.
/// Standard technique in game physics (Baumgarte stabilization, 1972).
#[inline]
pub fn baumgarte_correction(depth: f32, slop: f32, pct: f32, inv_mass_sum: f32) -> f32 {
    if inv_mass_sum < 1e-8 {
        return 0.0;
    }
    (depth - slop).max(0.0) * pct / inv_mass_sum
}

/// Coulomb friction clamp: limits tangential impulse to mu * |normal impulse|.
/// jt_clamped = clamp(jt, -|j_normal| * mu, |j_normal| * mu)
#[inline]
pub fn coulomb_friction_clamp(jt: f32, j_normal_abs: f32, mu: f32) -> f32 {
    jt.clamp(-j_normal_abs * mu, j_normal_abs * mu)
}

/// Union-find root with path compression. O(alpha(n)) amortized.
/// Used for island splitting in the Parallel Jacobi solver.
pub fn union_find_root(parent: &mut [usize], mut x: usize) -> usize {
    while parent[x] != x {
        parent[x] = parent[parent[x]]; // path halving
        x = parent[x];
    }
    x
}

/// Union-find merge by rank. O(alpha(n)) amortized.
pub fn union_find_merge(parent: &mut [usize], rank: &mut [u8], a: usize, b: usize) {
    let ra = union_find_root(parent, a);
    let rb = union_find_root(parent, b);
    if ra == rb {
        return;
    }
    if rank[ra] < rank[rb] {
        parent[ra] = rb;
    } else if rank[ra] > rank[rb] {
        parent[rb] = ra;
    } else {
        parent[rb] = ra;
        rank[ra] += 1;
    }
}

// ══════════════════════════════════════════════════════════════
// Additional Interpolation & Easing
// ══════════════════════════════════════════════════════════════

/// Quaternion spherical linear interpolation (slerp).
/// Interpolates between two unit quaternions [x, y, z, w] along the great arc.
pub fn slerp(a: [f32; 4], b: [f32; 4], t: f32) -> [f32; 4] {
    let mut dot = a[0] * b[0] + a[1] * b[1] + a[2] * b[2] + a[3] * b[3];
    let mut b = b;
    if dot < 0.0 {
        // shortest path
        b = [-b[0], -b[1], -b[2], -b[3]];
        dot = -dot;
    }
    if dot > 0.9995 {
        // near-parallel: lerp fallback
        let mut r = [0.0f32; 4];
        for i in 0..4 {
            r[i] = a[i] + (b[i] - a[i]) * t;
        }
        let len = (r[0] * r[0] + r[1] * r[1] + r[2] * r[2] + r[3] * r[3]).sqrt();
        if len > 1e-8 {
            for i in 0..4 {
                r[i] /= len;
            }
        }
        return r;
    }
    let theta = dot.clamp(-1.0, 1.0).acos();
    let sin_theta = theta.sin();
    let wa = ((1.0 - t) * theta).sin() / sin_theta;
    let wb = (t * theta).sin() / sin_theta;
    [
        a[0] * wa + b[0] * wb,
        a[1] * wa + b[1] * wb,
        a[2] * wa + b[2] * wb,
        a[3] * wa + b[3] * wb,
    ]
}

/// Range remapping: maps value from [in_min, in_max] to [out_min, out_max].
#[inline]
pub fn remap(value: f32, in_min: f32, in_max: f32, out_min: f32, out_max: f32) -> f32 {
    let t = (value - in_min) / (in_max - in_min);
    out_min + (out_max - out_min) * t
}

/// Cubic ease-in-out: smooth acceleration and deceleration.
/// t in [0, 1]. Returns smooth curve from 0 to 1.
#[inline]
pub fn ease_in_out_cubic(t: f32) -> f32 {
    if t < 0.5 {
        4.0 * t * t * t
    } else {
        1.0 - (-2.0 * t + 2.0).powi(3) / 2.0
    }
}

/// Saturating clamp to [0, 1].
#[inline]
pub fn clamp01(v: f32) -> f32 {
    v.clamp(0.0, 1.0)
}

// ══════════════════════════════════════════════════════════════
// Tests
// ══════════════════════════════════════════════════════════════

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

    #[test]
    fn complex_arithmetic() {
        let a = Complex::new(1.0, 2.0);
        let b = Complex::new(3.0, -1.0);
        let c = a.mul(b); // (1+2i)(3-i) = 3-i+6i-2i² = 5+5i
        assert!((c.re - 5.0).abs() < 1e-6);
        assert!((c.im - 5.0).abs() < 1e-6);
    }

    #[test]
    fn complex_euler() {
        let z = Complex::exp_i(std::f32::consts::PI); // e^{iπ} = -1
        assert!((z.re - (-1.0)).abs() < 1e-5);
        assert!(z.im.abs() < 1e-5);
    }

    #[test]
    fn complex_conjugate() {
        let z = Complex::new(3.0, 4.0);
        let zc = z.conj();
        assert_eq!(zc.re, 3.0);
        assert_eq!(zc.im, -4.0);
        assert!((z.mul(zc).im).abs() < 1e-6); // z·z* is real
    }

    #[test]
    fn distance_sq() {
        let a = [1.0, 2.0, 3.0];
        let b = [4.0, 6.0, 3.0];
        assert!((distance_sq_3d(&a, &b) - 25.0).abs() < 1e-6); // 9+16+0
    }

    #[test]
    fn gaussian_at_center() {
        assert!((gaussian_kernel(0.0, 1.0) - 1.0).abs() < 1e-6);
    }

    #[test]
    fn gaussian_decays() {
        let near = gaussian_kernel(1.0, 4.0);
        let far = gaussian_kernel(9.0, 4.0);
        assert!(near > far);
    }

    #[test]
    fn traveling_wave_oscillates() {
        let y0 = traveling_wave(0.0, 0.0, 1.0, 1.0, 1.0);
        let y1 = traveling_wave(std::f32::consts::PI, 0.0, 1.0, 1.0, 1.0);
        assert!(y0.abs() < 1e-6); // sin(0) = 0
        assert!(y1.abs() < 1e-6); // sin(π) ≈ 0
    }

    #[test]
    fn wave_number_from_wavelength() {
        let k = wave_number(1.0);
        assert!((k - std::f32::consts::TAU).abs() < 1e-5);
    }

    #[test]
    fn golden_ratio_value() {
        assert!((golden_ratio() - 1.6180339).abs() < 1e-5);
    }

    #[test]
    fn golden_angle_value() {
        assert!((golden_angle() - 2.3999632).abs() < 1e-4);
    }

    #[test]
    fn fibonacci_sphere_count() {
        let pts = fibonacci_sphere(128);
        assert_eq!(pts.len(), 128);
        for p in &pts {
            let len = vec3_len(*p);
            assert!((len - 1.0).abs() < 0.01);
        }
    }

    #[test]
    fn phyllotaxis_disc_count() {
        let pts = phyllotaxis_disc(64, 1.0, 0.0);
        assert_eq!(pts.len(), 64);
    }

    #[test]
    fn lerp_endpoints() {
        assert_eq!(lerp(0.0, 10.0, 0.0), 0.0);
        assert_eq!(lerp(0.0, 10.0, 1.0), 10.0);
        assert!((lerp(0.0, 10.0, 0.5) - 5.0).abs() < 1e-6);
    }

    #[test]
    fn exp_decay_converges() {
        let mut val = 0.0f32;
        for _ in 0..100 {
            val = exp_decay(val, 1.0, 8.0, 1.0 / 60.0);
        }
        assert!((val - 1.0).abs() < 0.01);
    }

    #[test]
    fn smoothstep_boundaries() {
        assert!((smoothstep(0.0, 1.0, -0.5)).abs() < 1e-6);
        assert!((smoothstep(0.0, 1.0, 1.5) - 1.0).abs() < 1e-6);
        assert!((smoothstep(0.0, 1.0, 0.5) - 0.5).abs() < 0.1); // ~0.5 at midpoint
    }

    #[test]
    fn splitmix64_deterministic() {
        let (a, _) = splitmix64(42);
        let (b, _) = splitmix64(42);
        assert_eq!(a, b);
    }

    #[test]
    fn splitmix64_f32_range() {
        let (val, _) = splitmix64_f32(42);
        assert!(val >= 0.0 && val < 1.0);
    }

    #[test]
    fn vec3_operations() {
        let a = [1.0, 0.0, 0.0];
        let b = [0.0, 1.0, 0.0];
        let cross = vec3_cross(a, b);
        assert!((cross[2] - 1.0).abs() < 1e-6); // x × y = z
    }

    #[test]
    fn half_neighborhood_14_cells() {
        assert_eq!(HALF_NEIGHBORHOOD.len(), 14);
    }

    #[test]
    fn fnv1a_deterministic() {
        assert_eq!(fnv1a_3d(1, 2, 3), fnv1a_3d(1, 2, 3));
        assert_ne!(fnv1a_3d(1, 2, 3), fnv1a_3d(3, 2, 1));
    }

    #[test]
    fn adiabatic_in_window() {
        assert!((adiabatic_fraction(1, 60, 10) - 0.1).abs() < 1e-6);
        assert_eq!(adiabatic_fraction(15, 60, 10), 0.0); // outside window
    }

    #[test]
    fn contact_impulse_separating() {
        // Bodies separating (v_along_n > 0): impulse should be negative (or zero for game use)
        let j = contact_impulse(1.0, 0.5, 2.0);
        assert!(j < 0.0); // -(1+0.5)*1.0/2.0 = -0.75
    }

    #[test]
    fn union_find_basic() {
        let mut parent: Vec<usize> = (0..5).collect();
        let mut rank = vec![0u8; 5];
        union_find_merge(&mut parent, &mut rank, 0, 1);
        union_find_merge(&mut parent, &mut rank, 2, 3);
        union_find_merge(&mut parent, &mut rank, 1, 3);
        assert_eq!(union_find_root(&mut parent, 0), union_find_root(&mut parent, 3));
    }

    #[test]
    fn slerp_endpoints() {
        let a = [0.0, 0.0, 0.0, 1.0]; // identity quaternion
        let b = [0.0, 1.0, 0.0, 0.0]; // 180° around Y
        let mid = slerp(a, b, 0.0);
        assert!((mid[3] - 1.0).abs() < 0.01); // at t=0, should be near a
    }

    #[test]
    fn remap_basic() {
        assert!((remap(5.0, 0.0, 10.0, 0.0, 100.0) - 50.0).abs() < 1e-6);
    }

    #[test]
    fn ease_cubic_endpoints() {
        assert!(ease_in_out_cubic(0.0).abs() < 1e-6);
        assert!((ease_in_out_cubic(1.0) - 1.0).abs() < 1e-6);
    }

    #[test]
    fn clamp01_bounds() {
        assert_eq!(clamp01(-0.5), 0.0);
        assert_eq!(clamp01(1.5), 1.0);
        assert!((clamp01(0.5) - 0.5).abs() < 1e-6);
    }
}