alice-sdf 1.11.0

//! Interval Arithmetic for SDF evaluation
//!
//! Evaluates SDF trees over spatial regions (AABBs) instead of single points.
//! Returns conservative distance bounds [lo, hi] guaranteed to contain
//! all possible distance values within the region.
//!
//! - If `result.lo > 0`: entire region is outside the surface (safe to skip)
//! - If `result.hi < 0`: entire region is inside the surface
//! - Otherwise: surface may cross through this region
//!
//! Author: Moroya Sakamoto

use crate::crispy::round_half_up;
use crate::types::SdfNode;
use glam::{Quat, Vec3};
use std::ops::{Add, Mul, Neg, Sub};

/// A closed interval [lo, hi] representing a range of possible values
#[derive(Clone, Copy, Debug)]
pub struct Interval {
    /// Lower bound
    pub lo: f32,
    /// Upper bound
    pub hi: f32,
}

impl Interval {
    /// Create a new interval
    #[inline(always)]
    pub fn new(lo: f32, hi: f32) -> Self {
        debug_assert!(lo <= hi + 1e-6, "lo ({}) > hi ({})", lo, hi);
        Self { lo, hi }
    }

    /// Create a point interval [v, v]
    #[inline(always)]
    pub const fn point(v: f32) -> Self {
        Self { lo: v, hi: v }
    }

    /// The entire real line
    pub const EVERYTHING: Self = Self {
        lo: f32::NEG_INFINITY,
        hi: f32::INFINITY,
    };

    /// Zero interval
    pub const ZERO: Self = Self { lo: 0.0, hi: 0.0 };

    /// Check if the interval is entirely positive
    #[inline(always)]
    pub fn is_positive(self) -> bool {
        self.lo > 0.0
    }

    /// Check if the interval is entirely negative
    #[inline(always)]
    pub fn is_negative(self) -> bool {
        self.hi < 0.0
    }

    /// Width of the interval (hi - lo)
    #[inline(always)]
    pub fn width(self) -> f32 {
        self.hi - self.lo
    }

    /// Midpoint of the interval
    #[inline(always)]
    pub fn midpoint(self) -> f32 {
        (self.lo + self.hi) * 0.5
    }

    /// Check if a scalar value is within the interval
    #[inline(always)]
    pub fn contains(self, v: f32) -> bool {
        v >= self.lo && v <= self.hi
    }

    /// Check if two intervals overlap
    #[inline(always)]
    pub fn overlaps(self, other: Self) -> bool {
        self.lo <= other.hi && other.lo <= self.hi
    }

    /// Intersection of two intervals (narrowest overlap)
    #[inline(always)]
    pub fn intersect(self, other: Self) -> Self {
        Self {
            lo: self.lo.max(other.lo),
            hi: self.hi.min(other.hi),
        }
    }

    /// Hull (union) of two intervals (widest encompassing)
    #[inline(always)]
    pub fn hull(self, other: Self) -> Self {
        Self {
            lo: self.lo.min(other.lo),
            hi: self.hi.max(other.hi),
        }
    }

    /// Absolute value of an interval
    #[inline(always)]
    pub fn abs(self) -> Self {
        if self.lo >= 0.0 {
            self
        } else if self.hi <= 0.0 {
            Self {
                lo: -self.hi,
                hi: -self.lo,
            }
        } else {
            Self {
                lo: 0.0,
                hi: self.hi.max(-self.lo),
            }
        }
    }

    /// Square root (clamped to non-negative)
    #[inline(always)]
    pub fn sqrt(self) -> Self {
        Self {
            lo: self.lo.max(0.0).sqrt(),
            hi: self.hi.max(0.0).sqrt(),
        }
    }

    /// Square of an interval
    #[inline(always)]
    pub fn sqr(self) -> Self {
        if self.lo >= 0.0 {
            Self {
                lo: self.lo * self.lo,
                hi: self.hi * self.hi,
            }
        } else if self.hi <= 0.0 {
            Self {
                lo: self.hi * self.hi,
                hi: self.lo * self.lo,
            }
        } else {
            Self {
                lo: 0.0,
                hi: (self.lo * self.lo).max(self.hi * self.hi),
            }
        }
    }

    /// Minimum of two intervals
    #[inline(always)]
    pub fn min(self, other: Self) -> Self {
        Self {
            lo: self.lo.min(other.lo),
            hi: self.hi.min(other.hi),
        }
    }

    /// Maximum of two intervals
    #[inline(always)]
    pub fn max(self, other: Self) -> Self {
        Self {
            lo: self.lo.max(other.lo),
            hi: self.hi.max(other.hi),
        }
    }

    /// Clamp to scalar range
    #[inline(always)]
    pub fn clamp(self, lo: f32, hi: f32) -> Self {
        Self {
            lo: self.lo.clamp(lo, hi),
            hi: self.hi.clamp(lo, hi),
        }
    }

    /// Expand interval by a constant in both directions
    #[inline(always)]
    pub fn expand(self, amount: f32) -> Self {
        Self {
            lo: self.lo - amount,
            hi: self.hi + amount,
        }
    }
}

impl Add for Interval {
    type Output = Self;
    #[inline(always)]
    fn add(self, rhs: Self) -> Self {
        Self {
            lo: self.lo + rhs.lo,
            hi: self.hi + rhs.hi,
        }
    }
}

impl Sub for Interval {
    type Output = Self;
    #[inline(always)]
    fn sub(self, rhs: Self) -> Self {
        Self {
            lo: self.lo - rhs.hi,
            hi: self.hi - rhs.lo,
        }
    }
}

impl Mul for Interval {
    type Output = Self;
    #[inline(always)]
    fn mul(self, rhs: Self) -> Self {
        let a = self.lo * rhs.lo;
        let b = self.lo * rhs.hi;
        let c = self.hi * rhs.lo;
        let d = self.hi * rhs.hi;
        Self {
            lo: a.min(b).min(c).min(d),
            hi: a.max(b).max(c).max(d),
        }
    }
}

impl Neg for Interval {
    type Output = Self;
    #[inline(always)]
    fn neg(self) -> Self {
        Self {
            lo: -self.hi,
            hi: -self.lo,
        }
    }
}

impl Mul<f32> for Interval {
    type Output = Self;
    #[inline(always)]
    fn mul(self, rhs: f32) -> Self {
        if rhs >= 0.0 {
            Self {
                lo: self.lo * rhs,
                hi: self.hi * rhs,
            }
        } else {
            Self {
                lo: self.hi * rhs,
                hi: self.lo * rhs,
            }
        }
    }
}

impl Sub<f32> for Interval {
    type Output = Self;
    #[inline(always)]
    fn sub(self, rhs: f32) -> Self {
        Self {
            lo: self.lo - rhs,
            hi: self.hi - rhs,
        }
    }
}

// ============================================================
// Vec3Interval: 3D interval (axis-aligned box)
// ============================================================

/// 3D interval representing an axis-aligned box region
#[derive(Clone, Copy, Debug)]
pub struct Vec3Interval {
    /// X interval
    pub x: Interval,
    /// Y interval
    pub y: Interval,
    /// Z interval
    pub z: Interval,
}

impl Vec3Interval {
    /// Create from min/max corners
    #[inline]
    pub fn from_bounds(min: Vec3, max: Vec3) -> Self {
        Self {
            x: Interval::new(min.x, max.x),
            y: Interval::new(min.y, max.y),
            z: Interval::new(min.z, max.z),
        }
    }

    /// Translate by subtracting offset
    #[inline(always)]
    pub fn translate(self, offset: Vec3) -> Self {
        Self {
            x: self.x - Interval::point(offset.x),
            y: self.y - Interval::point(offset.y),
            z: self.z - Interval::point(offset.z),
        }
    }

    /// Rotate by quaternion (conservative: bounding box of rotated corners)
    #[inline]
    pub fn rotate(self, rotation: Quat) -> Self {
        let inv_rot = rotation.conjugate();
        let corners = [
            Vec3::new(self.x.lo, self.y.lo, self.z.lo),
            Vec3::new(self.x.lo, self.y.lo, self.z.hi),
            Vec3::new(self.x.lo, self.y.hi, self.z.lo),
            Vec3::new(self.x.lo, self.y.hi, self.z.hi),
            Vec3::new(self.x.hi, self.y.lo, self.z.lo),
            Vec3::new(self.x.hi, self.y.lo, self.z.hi),
            Vec3::new(self.x.hi, self.y.hi, self.z.lo),
            Vec3::new(self.x.hi, self.y.hi, self.z.hi),
        ];
        let first = inv_rot * corners[0];
        let mut min = first;
        let mut max = first;
        for &c in &corners[1..] {
            let r = inv_rot * c;
            min = min.min(r);
            max = max.max(r);
        }
        Self::from_bounds(min, max)
    }

    /// 3D interval length: sqrt(x² + y² + z²)
    #[inline]
    pub fn length(self) -> Interval {
        (self.x.sqr() + self.y.sqr() + self.z.sqr()).sqrt()
    }

    /// 2D length in XZ: sqrt(x² + z²)
    #[inline]
    pub fn length_xz(self) -> Interval {
        (self.x.sqr() + self.z.sqr()).sqrt()
    }

    /// Component-wise mirror (abs specified axes)
    #[inline]
    pub fn mirror(self, axes: Vec3) -> Self {
        Self {
            x: if axes.x != 0.0 { self.x.abs() } else { self.x },
            y: if axes.y != 0.0 { self.y.abs() } else { self.y },
            z: if axes.z != 0.0 { self.z.abs() } else { self.z },
        }
    }
}

// ============================================================
// Interval SDF evaluation
// ============================================================

/// Evaluate an SDF tree over a spatial region, returning distance bounds.
///
/// For all points p in the region: `result.lo <= sdf(p) <= result.hi`
pub fn eval_interval(node: &SdfNode, bounds: Vec3Interval) -> Interval {
    match node {
        // ============ Exact-IA Primitives ============
        SdfNode::Sphere { radius } => bounds.length() - Interval::point(*radius),

        SdfNode::Box3d { half_extents } => {
            let h = *half_extents;
            let qx = bounds.x.abs() - Interval::point(h.x);
            let qy = bounds.y.abs() - Interval::point(h.y);
            let qz = bounds.z.abs() - Interval::point(h.z);
            let outer = Vec3Interval {
                x: qx.max(Interval::ZERO),
                y: qy.max(Interval::ZERO),
                z: qz.max(Interval::ZERO),
            }
            .length();
            let inner = qx.max(qy).max(qz).min(Interval::ZERO);
            outer + inner
        }

        SdfNode::Cylinder {
            radius,
            half_height,
        } => {
            let dx = bounds.length_xz() - Interval::point(*radius);
            let dy = bounds.y.abs() - Interval::point(*half_height);
            let ox = dx.max(Interval::ZERO);
            let oy = dy.max(Interval::ZERO);
            (ox.sqr() + oy.sqr()).sqrt() + dx.max(dy).min(Interval::ZERO)
        }

        SdfNode::Torus {
            major_radius,
            minor_radius,
        } => {
            let qx = bounds.length_xz() - Interval::point(*major_radius);
            (qx.sqr() + bounds.y.sqr()).sqrt() - Interval::point(*minor_radius)
        }

        SdfNode::Plane { normal, distance } => {
            bounds.x * Interval::point(normal.x)
                + bounds.y * Interval::point(normal.y)
                + bounds.z * Interval::point(normal.z)
                - Interval::point(*distance)
        }

        SdfNode::Capsule {
            point_a,
            point_b,
            radius,
        } => {
            let ab = *point_b - *point_a;
            let ab_sq = ab.length_squared();
            if ab_sq < 1e-10 {
                let p = bounds.translate(*point_a);
                return p.length() - Interval::point(*radius);
            }
            let pa = Vec3Interval {
                x: bounds.x - Interval::point(point_a.x),
                y: bounds.y - Interval::point(point_a.y),
                z: bounds.z - Interval::point(point_a.z),
            };
            let dot = pa.x * Interval::point(ab.x)
                + pa.y * Interval::point(ab.y)
                + pa.z * Interval::point(ab.z);
            let t = (dot * (1.0 / ab_sq)).clamp(0.0, 1.0);
            let cx = pa.x - Interval::point(ab.x) * t;
            let cy = pa.y - Interval::point(ab.y) * t;
            let cz = pa.z - Interval::point(ab.z) * t;
            Vec3Interval {
                x: cx,
                y: cy,
                z: cz,
            }
            .length()
                - Interval::point(*radius)
        }

        SdfNode::InfiniteCylinder { radius } => bounds.length_xz() - Interval::point(*radius),

        SdfNode::RoundedBox {
            half_extents,
            round_radius,
        } => {
            let h = *half_extents;
            let qx = bounds.x.abs() - Interval::point(h.x);
            let qy = bounds.y.abs() - Interval::point(h.y);
            let qz = bounds.z.abs() - Interval::point(h.z);
            let outer = Vec3Interval {
                x: qx.max(Interval::ZERO),
                y: qy.max(Interval::ZERO),
                z: qz.max(Interval::ZERO),
            }
            .length();
            let inner = qx.max(qy).max(qz).min(Interval::ZERO);
            outer + inner - Interval::point(*round_radius)
        }

        // ============ Conservative-IA Primitives (bounding sphere) ============
        SdfNode::Cone {
            radius,
            half_height,
        } => ia_bsphere(bounds, (2.0 * half_height).hypot(*radius)),
        // iq's ellipsoid bound `k0 (k0 - 1) / k1` scales the unit-sphere distance
        // by at most max(r) / min(r), so its gradient is bounded by that ratio.
        // Exact SDF since 1.11.0 (Eberly nearest point)
        SdfNode::Ellipsoid { .. } => ia_lipschitz(node, bounds, 1.0),
        SdfNode::RoundedCone {
            r1,
            r2,
            half_height,
        } => ia_bsphere(bounds, (2.0 * half_height).hypot(r1.max(*r2))),
        SdfNode::Pyramid { half_height } => ia_bsphere(
            bounds,
            (2.0 * half_height).mul_add(2.0 * half_height, 0.25).sqrt(),
        ),
        SdfNode::Octahedron { size } => ia_bsphere(bounds, *size),
        SdfNode::HexPrism { .. } => ia_lipschitz(node, bounds, 1.0),
        SdfNode::Link {
            half_length,
            r1,
            r2,
        } => ia_bsphere(bounds, half_length + r1 + r2),
        SdfNode::Triangle { .. } | SdfNode::Bezier { .. } => Interval::EVERYTHING,
        SdfNode::CappedCone {
            half_height,
            r1,
            r2,
        } => ia_bsphere(bounds, (2.0 * half_height).hypot(r1.max(*r2))),
        SdfNode::CappedTorus {
            major_radius,
            minor_radius,
            ..
        } => ia_bsphere(bounds, major_radius + minor_radius),
        SdfNode::RoundedCylinder { .. } => ia_lipschitz(node, bounds, 1.0),
        // max of unit-normal plane distances: 1-Lipschitz lower bound
        SdfNode::TriangularPrism { .. } => ia_lipschitz(node, bounds, 1.0),
        SdfNode::CutSphere { radius, .. } => ia_bsphere(bounds, *radius),
        SdfNode::CutHollowSphere { .. } => ia_lipschitz(node, bounds, 1.0),
        SdfNode::DeathStar { ra, rb, .. } => ia_bsphere(bounds, ra.max(*rb)),
        SdfNode::SolidAngle { radius, .. } => ia_bsphere(bounds, *radius),
        SdfNode::Rhombus {
            la,
            lb,
            half_height,
            ..
        } => ia_bsphere(bounds, la.max(*lb).hypot(*half_height)),
        SdfNode::Horseshoe {
            radius,
            half_length,
            width,
            thickness,
            ..
        } => ia_bsphere(bounds, radius + half_length + width + thickness),
        SdfNode::Vesica { radius, .. } => ia_bsphere(bounds, *radius),
        SdfNode::InfiniteCone { .. } => Interval::EVERYTHING,
        // Triply periodic minimal surfaces: not distance fields, but Lipschitz
        // with the constants `eval_lipschitz` pins (√3 … 7), so the centre
        // sample ± L·ρ is a sound, finite enclosure. Until 1.10.3 these
        // returned EVERYTHING and interval-based pruning silently did nothing
        // for TPMS scenes (external review 2026-09-15, SDF-R2-5).
        SdfNode::Gyroid { .. }
        | SdfNode::SchwarzP { .. }
        | SdfNode::DiamondSurface { .. }
        | SdfNode::Neovius { .. }
        | SdfNode::Lidinoid { .. }
        | SdfNode::IWP { .. }
        | SdfNode::FRD { .. }
        | SdfNode::FischerKochS { .. }
        | SdfNode::PMY { .. } => ia_lipschitz(node, bounds, eval_lipschitz(node)),
        SdfNode::Heart { size } => ia_bsphere(bounds, *size * 2.0),
        SdfNode::Tube { .. } => ia_lipschitz(node, bounds, 1.0),
        SdfNode::Barrel {
            radius,
            half_height,
            bulge,
        } => ia_bsphere(bounds, (radius + bulge.abs()).hypot(*half_height)),
        SdfNode::Diamond {
            radius,
            half_height,
        } => ia_bsphere(bounds, radius.hypot(*half_height)),
        SdfNode::ChamferedCube { half_extents, .. } => ia_bsphere(bounds, half_extents.length()),
        // Implicit (non-distance) function: interval arithmetic on the formula
        // itself; every step is monotone on non-negative inputs.
        SdfNode::Superellipsoid {
            half_extents,
            e1,
            e2,
        } => {
            let e1 = e1.max(0.02);
            let e2 = e2.max(0.02);
            let q = |iv: Interval, h: f32| {
                let a = (iv * Interval::point(1.0 / h)).abs();
                Interval::new(a.lo.max(1e-10), a.hi.max(1e-10))
            };
            let (qx, qy, qz) = (
                q(bounds.x, half_extents.x),
                q(bounds.y, half_extents.y),
                q(bounds.z, half_extents.z),
            );
            let m1 = 2.0 / e2;
            let m2 = 2.0 / e1;
            let w = powf_nonneg(qx, m1) + powf_nonneg(qz, m1);
            let v = powf_nonneg(w, e2 / e1) + powf_nonneg(qy, m2);
            let f = powf_nonneg(v, e1 * 0.5);
            let min_extent = half_extents.x.min(half_extents.y.min(half_extents.z));
            (f - Interval::point(1.0)) * Interval::point(min_extent * 0.5)
        }
        SdfNode::RoundedX {
            width,
            round_radius,
            half_height,
        } => ia_bsphere(bounds, (width + round_radius).max(*half_height)),
        SdfNode::Pie {
            radius,
            half_height,
            ..
        } => ia_bsphere(bounds, radius.hypot(*half_height)),
        SdfNode::Trapezoid {
            r1,
            r2,
            trap_height,
            half_depth,
        } => ia_bsphere(
            bounds,
            (r1.max(*r2).mul_add(r1.max(*r2), trap_height.powi(2)) + half_depth.powi(2)).sqrt(),
        ),
        SdfNode::Parallelogram {
            width,
            para_height,
            skew,
            half_depth,
        } => ia_bsphere(
            bounds,
            ((width + skew.abs()).mul_add(width + skew.abs(), para_height.powi(2))
                + half_depth.powi(2))
            .sqrt(),
        ),
        // Piecewise (rect / min(rect, dome) above the dome centre): the value is
        // always between the union and the rectangle alone, both 1-Lipschitz.
        SdfNode::Tunnel {
            width,
            height_2d,
            half_depth,
        } => {
            let px = bounds.x.abs();
            let py = bounds.y;
            let dx = px - Interval::point(*width);
            let dy = py.abs() - Interval::point(*height_2d);
            let d_rect = len2(dx.max(Interval::ZERO), dy.max(Interval::ZERO))
                + dx.max(dy).min(Interval::ZERO);
            let d_circle = len2(px, py - Interval::point(*height_2d)) - Interval::point(*width);
            let d_2d = Interval::new(d_rect.min(d_circle).lo, d_rect.hi);
            let d_z = bounds.z.abs() - Interval::point(*half_depth);
            d_2d.max(d_z).min(Interval::ZERO)
                + len2(d_2d.max(Interval::ZERO), d_z.max(Interval::ZERO))
        }
        SdfNode::UnevenCapsule {
            r1,
            r2,
            cap_height,
            half_depth,
        } => ia_bsphere(bounds, r1.max(*r2) + cap_height + half_depth),
        // apex at y = √3·(ra − rb) + ra (IQ sdEgg), base radius ra
        SdfNode::Egg { ra, rb } => ia_bsphere(bounds, 1.732_050_8f32.mul_add(ra - rb, *ra)),
        SdfNode::ArcShape {
            radius,
            thickness,
            half_height,
            ..
        } => ia_bsphere(bounds, (radius + thickness).hypot(*half_height)),
        SdfNode::Moon {
            ra,
            rb,
            half_height,
            ..
        } => ia_bsphere(bounds, ra.max(*rb).hypot(*half_height)),
        SdfNode::CrossShape { .. } => ia_lipschitz(node, bounds, 1.0),
        // Exact SDF (IQ sdBlobbyCross): centre sample ± ρ
        SdfNode::BlobbyCross { .. } => ia_lipschitz(node, bounds, 1.0),
        SdfNode::ParabolaSegment {
            width,
            para_height,
            half_depth,
        } => ia_bsphere(
            bounds,
            (width.powi(2) + para_height.powi(2) + half_depth.powi(2)).sqrt(),
        ),
        SdfNode::RegularPolygon {
            radius,
            half_height,
            ..
        } => ia_bsphere(bounds, radius.hypot(*half_height)),
        SdfNode::StarPolygon {
            radius,
            half_height,
            ..
        } => ia_bsphere(bounds, radius.hypot(*half_height)),
        SdfNode::Stairs {
            step_width,
            step_height,
            n_steps,
            half_depth,
        } => {
            let extent = ((step_width * n_steps)
                .mul_add(step_width * n_steps, (step_height * n_steps).powi(2))
                + half_depth.powi(2))
            .sqrt();
            ia_bsphere(bounds, extent)
        }
        // Nearest-wrap selection makes the tube term discontinuous near the axis,
        // so bound it directly: |dy| ≤ pitch / 2 after wrap selection.
        SdfNode::Helix {
            major_r,
            minor_r,
            pitch,
            half_height,
        } => {
            let d_radial = bounds.length_xz() - Interval::point(*major_r);
            let dy = Interval::new(0.0, pitch.abs() * 0.5);
            let d_tube =
                Interval::new(d_radial.abs().lo - minor_r, len2(d_radial, dy).hi - minor_r);
            let d_cap = bounds.y.abs() - Interval::point(*half_height);
            d_tube.max(d_cap)
        }
        // Generalised distance functions (max of unit-normal planes): 1-Lipschitz
        SdfNode::Tetrahedron { .. }
        | SdfNode::Dodecahedron { .. }
        | SdfNode::Icosahedron { .. }
        | SdfNode::TruncatedOctahedron { .. }
        | SdfNode::TruncatedIcosahedron { .. }
        | SdfNode::BoxFrame { .. } => ia_lipschitz(node, bounds, 1.0),

        // ============ 2D Primitives ============
        // Extruded 2D SDFs: 1-Lipschitz (the old arms mis-modelled the extrusion)
        SdfNode::Circle2D { .. } => ia_lipschitz(node, bounds, 1.0),
        SdfNode::Rect2D { .. } => ia_lipschitz(node, bounds, 1.0),
        SdfNode::Segment2D { .. } | SdfNode::Polygon2D { .. } => Interval::EVERYTHING,
        SdfNode::RoundedRect2D { .. } | SdfNode::Annular2D { .. } => {
            ia_lipschitz(node, bounds, 1.0)
        }
        SdfNode::Terrain { amplitude, .. } => {
            // 地形は y - height なので、amplitude ベースの概算
            Interval::new(bounds.y.lo - amplitude, bounds.y.hi + amplitude)
        }

        // ============ Operations ============
        SdfNode::Union { a, b } => eval_interval(a, bounds).min(eval_interval(b, bounds)),
        SdfNode::Intersection { a, b } => eval_interval(a, bounds).max(eval_interval(b, bounds)),
        SdfNode::Subtraction { a, b } => eval_interval(a, bounds).max(-eval_interval(b, bounds)),
        SdfNode::SmoothUnion { a, b, k } => {
            let sharp = eval_interval(a, bounds).min(eval_interval(b, bounds));
            Interval::new(sharp.lo - k * 0.25, sharp.hi)
        }
        SdfNode::SmoothIntersection { a, b, k } => {
            let sharp = eval_interval(a, bounds).max(eval_interval(b, bounds));
            Interval::new(sharp.lo, sharp.hi + k * 0.25)
        }
        SdfNode::SmoothSubtraction { a, b, k } => {
            let sharp = eval_interval(a, bounds).max(-eval_interval(b, bounds));
            Interval::new(sharp.lo, sharp.hi + k * 0.25)
        }
        // chamfer_min(a, b, r) = min(a, b, (a + b) / √2 - r): exact interval arithmetic
        SdfNode::ChamferUnion { a, b, r } => {
            let ia = eval_interval(a, bounds);
            let ib = eval_interval(b, bounds);
            ia.min(ib).min(
                (ia + ib) * Interval::point(std::f32::consts::FRAC_1_SQRT_2) - Interval::point(*r),
            )
        }
        // chamfer_max(a, b, r) = max(a, b, (a + b) / √2 + r): exact interval arithmetic
        SdfNode::ChamferIntersection { a, b, r } => {
            let ia = eval_interval(a, bounds);
            let ib = eval_interval(b, bounds);
            ia.max(ib).max(
                (ia + ib) * Interval::point(std::f32::consts::FRAC_1_SQRT_2) + Interval::point(*r),
            )
        }
        SdfNode::ChamferSubtraction { a, b, r } => {
            let ia = eval_interval(a, bounds);
            let ib = Interval::ZERO - eval_interval(b, bounds);
            ia.max(ib).max(
                (ia + ib) * Interval::point(std::f32::consts::FRAC_1_SQRT_2) + Interval::point(*r),
            )
        }
        // stairs_*(a, b) is 1-Lipschitz in (a, b) (rotation + mod + min/abs of unit
        // slopes), so evaluate at the (a, b) rectangle centre and widen by its half-diagonal.
        SdfNode::StairsUnion { a, b, r, n } => {
            let ia = eval_interval(a, bounds);
            let ib = eval_interval(b, bounds);
            if !(ia.lo.is_finite() && ia.hi.is_finite() && ib.lo.is_finite() && ib.hi.is_finite()) {
                return Interval::EVERYTHING;
            }
            let (ac, bc) = (ia.midpoint(), ib.midpoint());
            let rho = (ia.width() * 0.5).hypot(ib.width() * 0.5);
            let d = crate::operations::sdf_stairs_union(ac, bc, *r, *n);
            Interval::new(d - rho, d + rho)
        }
        // stairs_*(a, b) is 1-Lipschitz in (a, b) (rotation + mod + min/abs of unit
        // slopes), so evaluate at the (a, b) rectangle centre and widen by its half-diagonal.
        SdfNode::StairsIntersection { a, b, r, n } => {
            let ia = eval_interval(a, bounds);
            let ib = eval_interval(b, bounds);
            if !(ia.lo.is_finite() && ia.hi.is_finite() && ib.lo.is_finite() && ib.hi.is_finite()) {
                return Interval::EVERYTHING;
            }
            let (ac, bc) = (ia.midpoint(), ib.midpoint());
            let rho = (ia.width() * 0.5).hypot(ib.width() * 0.5);
            let d = crate::operations::sdf_stairs_intersection(ac, bc, *r, *n);
            Interval::new(d - rho, d + rho)
        }
        // stairs_*(a, b) is 1-Lipschitz in (a, b) (rotation + mod + min/abs of unit
        // slopes), so evaluate at the (a, b) rectangle centre and widen by its half-diagonal.
        SdfNode::StairsSubtraction { a, b, r, n } => {
            let ia = eval_interval(a, bounds);
            let ib = Interval::ZERO - eval_interval(b, bounds);
            if !(ia.lo.is_finite() && ia.hi.is_finite() && ib.lo.is_finite() && ib.hi.is_finite()) {
                return Interval::EVERYTHING;
            }
            let (ac, bc) = (ia.midpoint(), ib.midpoint());
            let rho = (ia.width() * 0.5).hypot(ib.width() * 0.5);
            let d = crate::operations::sdf_stairs_intersection(ac, bc, *r, *n);
            Interval::new(d - rho, d + rho)
        }
        SdfNode::XOR { a, b } => {
            let ia = eval_interval(a, bounds);
            let ib = eval_interval(b, bounds);
            // XOR = min(a,b).max(-max(a,b)); conservative: union of both bounds
            let min_ab = ia.min(ib);
            let neg_max = -(ia.max(ib));
            min_ab.max(neg_max)
        }
        SdfNode::Morph { a, b, t } => {
            let ia = eval_interval(a, bounds);
            let ib = eval_interval(b, bounds);
            ia * (1.0 - *t) + ib * *t
        }
        SdfNode::ColumnsUnion { a, b, r, .. } => {
            let sharp = eval_interval(a, bounds).min(eval_interval(b, bounds));
            Interval::new(sharp.lo - r, sharp.hi)
        }
        SdfNode::ColumnsIntersection { a, b, r, .. } => {
            let sharp = eval_interval(a, bounds).max(eval_interval(b, bounds));
            Interval::new(sharp.lo, sharp.hi + r)
        }
        SdfNode::ColumnsSubtraction { a, b, r, .. } => {
            let sharp = eval_interval(a, bounds).max(-eval_interval(b, bounds));
            Interval::new(sharp.lo, sharp.hi + r)
        }
        SdfNode::Pipe { a, b, r } => {
            let ia = eval_interval(a, bounds);
            let ib = eval_interval(b, bounds);
            // pipe = sqrt(a²+b²) - r; conservative bound
            let sq = ia * ia + ib * ib;
            Interval::new(sq.lo.max(0.0).sqrt() - r, sq.hi.max(0.0).sqrt() + r.abs())
        }
        SdfNode::Engrave { a, b, r } => {
            let ia = eval_interval(a, bounds);
            let ib = eval_interval(b, bounds);
            // engrave = max(a, (a+r-|b|)*FRAC_1_SQRT_2)
            let abs_b = ib.abs();
            let inner = Interval::new(
                (ia.lo + r - abs_b.hi) * std::f32::consts::FRAC_1_SQRT_2,
                (ia.hi + r - abs_b.lo).max(0.0) * std::f32::consts::FRAC_1_SQRT_2,
            );
            ia.max(inner)
        }
        SdfNode::Groove { a, b, ra, rb } => {
            let ia = eval_interval(a, bounds);
            let ib = eval_interval(b, bounds);
            // groove = max(a, min(a+ra, rb-|b|))
            let abs_b = ib.abs();
            let arm1 = Interval::new(ia.lo + ra, ia.hi + ra);
            let arm2 = Interval::new(rb - abs_b.hi, rb - abs_b.lo);
            ia.max(arm1.min(arm2))
        }
        SdfNode::Tongue { a, b, ra, rb } => {
            let ia = eval_interval(a, bounds);
            let ib = eval_interval(b, bounds);
            // tongue = min(a, max(a-ra, |b|-rb))
            let abs_b = ib.abs();
            let arm1 = Interval::new(ia.lo - ra, ia.hi - ra);
            let arm2 = Interval::new(abs_b.lo - rb, abs_b.hi - rb);
            ia.min(arm1.max(arm2))
        }
        // k·ln(e^{a/k} + e^{b/k}) lies in [max(a, b), max(a, b) + k·ln 2] (and the
        // union form in [min(a, b) - k·ln 2, min(a, b)]).
        SdfNode::ExpSmoothUnion { a, b, k } => {
            let sharp = eval_interval(a, bounds).min(eval_interval(b, bounds));
            Interval::new(sharp.lo - k.max(1e-6) * std::f32::consts::LN_2, sharp.hi)
        }
        SdfNode::ExpSmoothIntersection { a, b, k } => {
            let sharp = eval_interval(a, bounds).max(eval_interval(b, bounds));
            Interval::new(sharp.lo, sharp.hi + k.max(1e-6) * std::f32::consts::LN_2)
        }
        SdfNode::ExpSmoothSubtraction { a, b, k } => {
            let sharp = eval_interval(a, bounds).max(-eval_interval(b, bounds));
            Interval::new(sharp.lo, sharp.hi + k.max(1e-6) * std::f32::consts::LN_2)
        }
        // Affine domain map (`modifier_shear`): interval image is exact
        SdfNode::Shear { child, shear } => {
            let x = bounds.x;
            let y = bounds.y - x * Interval::point(shear.x);
            let z = bounds.z - x * Interval::point(shear.y) - bounds.y * Interval::point(shear.z);
            eval_interval(child, Vec3Interval { x, y, z })
        }
        SdfNode::Animated { child, .. } => eval_interval(child, bounds),

        // ============ Transforms ============
        SdfNode::Translate { child, offset } => eval_interval(child, bounds.translate(*offset)),
        SdfNode::Rotate { child, rotation } => eval_interval(child, bounds.rotate(*rotation)),
        SdfNode::Scale { child, factor } => {
            let inv = 1.0 / *factor;
            let scaled = Vec3Interval {
                x: bounds.x * inv,
                y: bounds.y * inv,
                z: bounds.z * inv,
            };
            eval_interval(child, scaled) * *factor
        }
        SdfNode::ScaleNonUniform { child, factors } => {
            let min_f = factors.x.min(factors.y).min(factors.z);
            let max_f = factors.x.max(factors.y).max(factors.z);
            let scaled = Vec3Interval {
                x: bounds.x * (1.0 / factors.x),
                y: bounds.y * (1.0 / factors.y),
                z: bounds.z * (1.0 / factors.z),
            };
            let d = eval_interval(child, scaled);
            Interval::new(d.lo * min_f, d.hi * max_f)
        }
        SdfNode::ProjectiveTransform {
            child,
            lipschitz_bound,
            ..
        } => {
            // Conservative: scale by Lipschitz bound
            let child_interval = eval_interval(child, bounds);
            Interval::new(
                child_interval.lo / lipschitz_bound,
                child_interval.hi / lipschitz_bound,
            )
        }
        // Outside the lattice bbox the deformation is the identity (1.9.2); inside it
        // the trilinear warp has no usable Jacobian bound, so only the identity
        // region gets a finite interval.
        SdfNode::LatticeDeform {
            child,
            bbox_min,
            bbox_max,
            ..
        } => {
            let outside = bounds.x.hi < bbox_min.x
                || bounds.x.lo > bbox_max.x
                || bounds.y.hi < bbox_min.y
                || bounds.y.lo > bbox_max.y
                || bounds.z.hi < bbox_min.z
                || bounds.z.lo > bbox_max.z;
            if outside {
                eval_interval(child, bounds)
            } else {
                Interval::EVERYTHING
            }
        }
        SdfNode::SdfSkinning { child, .. } => {
            // LBS is approximately distance-preserving
            eval_interval(child, bounds)
        }

        // ============ Modifiers ============
        SdfNode::Twist { child, strength: _ } => {
            // After twist, XZ can reach any rotation → expand to circular
            let max_xz = bounds.length_xz().hi;
            let expanded = Vec3Interval {
                x: Interval::new(-max_xz, max_xz),
                y: bounds.y,
                z: Interval::new(-max_xz, max_xz),
            };
            eval_interval(child, expanded)
        }
        SdfNode::Bend { child, .. } => {
            let max_r = bounds.length().hi;
            eval_interval(
                child,
                Vec3Interval::from_bounds(Vec3::splat(-max_r), Vec3::splat(max_r)),
            )
        }
        SdfNode::RepeatInfinite { child, spacing } => {
            let hs = *spacing * 0.5;
            let cell = Vec3Interval {
                x: if spacing.x > 0.0 {
                    Interval::new(-hs.x, hs.x)
                } else {
                    bounds.x
                },
                y: if spacing.y > 0.0 {
                    Interval::new(-hs.y, hs.y)
                } else {
                    bounds.y
                },
                z: if spacing.z > 0.0 {
                    Interval::new(-hs.z, hs.z)
                } else {
                    bounds.z
                },
            };
            eval_interval(child, cell)
        }
        // q = p - clamp(round(p / s), -limit, limit) · s: the cell index over the box
        // is an integer interval, so the image is [p.lo - cmax·s, p.hi - cmin·s]
        // (this also covers boxes outside the finite repetition range).
        SdfNode::RepeatFinite {
            child,
            count,
            spacing,
        } => {
            let axis = |iv: Interval, s: f32, cnt: u32| {
                if s <= 0.0 {
                    return iv;
                }
                let limit = cnt as f32 * 0.5;
                // `* (1 / s)` like every point evaluator (tie parity)
                let inv = 1.0 / s;
                let cmin = round_half_up(iv.lo * inv).clamp(-limit, limit);
                let cmax = round_half_up(iv.hi * inv).clamp(-limit, limit);
                Interval::new(iv.lo - cmax * s, iv.hi - cmin * s)
            };
            eval_interval(
                child,
                Vec3Interval {
                    x: axis(bounds.x, spacing.x, count[0]),
                    y: axis(bounds.y, spacing.y, count[1]),
                    z: axis(bounds.z, spacing.z, count[2]),
                },
            )
        }
        SdfNode::Noise {
            child, amplitude, ..
        } => eval_interval(child, bounds).expand(*amplitude),
        SdfNode::Round { child, radius } => eval_interval(child, bounds) - Interval::point(*radius),
        SdfNode::Onion { child, thickness } => {
            eval_interval(child, bounds).abs() - Interval::point(*thickness)
        }
        SdfNode::Elongate { child, amount } => {
            let qx = bounds.x - bounds.x.clamp(-amount.x, amount.x);
            let qy = bounds.y - bounds.y.clamp(-amount.y, amount.y);
            let qz = bounds.z - bounds.z.clamp(-amount.z, amount.z);
            eval_interval(
                child,
                Vec3Interval {
                    x: qx,
                    y: qy,
                    z: qz,
                },
            )
        }
        SdfNode::Mirror { child, axes } => eval_interval(child, bounds.mirror(*axes)),
        // |p| sorted descending: the max coordinate lies in [max lo, max hi], the
        // min in [min lo, min hi] and the middle one in the hull.
        SdfNode::OctantMirror { child } => {
            let (ax, ay, az) = (bounds.x.abs(), bounds.y.abs(), bounds.z.abs());
            let lo_min = ax.lo.min(ay.lo).min(az.lo);
            let lo_max = ax.lo.max(ay.lo).max(az.lo);
            let hi_min = ax.hi.min(ay.hi).min(az.hi);
            let hi_max = ax.hi.max(ay.hi).max(az.hi);
            eval_interval(
                child,
                Vec3Interval {
                    x: Interval::new(lo_max, hi_max),
                    y: Interval::new(lo_min, hi_max),
                    z: Interval::new(lo_min, hi_min),
                },
            )
        }
        SdfNode::Revolution { child, offset } => {
            let qx = bounds.length_xz() - Interval::point(*offset);
            eval_interval(
                child,
                Vec3Interval {
                    x: qx,
                    y: bounds.y,
                    z: Interval::ZERO,
                },
            )
        }
        SdfNode::Extrude { child, half_height } => {
            let flat = Vec3Interval {
                x: bounds.x,
                y: bounds.y,
                z: Interval::ZERO,
            };
            let d = eval_interval(child, flat);
            let wy = bounds.z.abs() - Interval::point(*half_height);
            let ox = d.max(Interval::ZERO);
            let oy = wy.max(Interval::ZERO);
            (ox.sqr() + oy.sqr()).sqrt() + d.max(wy).min(Interval::ZERO)
        }
        // q = (x·s, y, z·s) with s = 1 / (1 - factor·y): if the denominator can
        // reach 0 on the box the map is singular (the evaluation itself blows up).
        SdfNode::Taper { child, factor } => {
            let denom = Interval::point(1.0) - bounds.y * Interval::point(*factor);
            if denom.lo <= 0.0 && denom.hi >= 0.0 {
                return Interval::EVERYTHING;
            }
            let s = Interval::new(
                (1.0 / denom.lo).min(1.0 / denom.hi),
                (1.0 / denom.lo).max(1.0 / denom.hi),
            );
            eval_interval(
                child,
                Vec3Interval {
                    x: bounds.x * s,
                    y: bounds.y,
                    z: bounds.z * s,
                },
            )
        }
        SdfNode::Displacement { child, strength } => {
            eval_interval(child, bounds).expand(strength.abs())
        }
        SdfNode::SineDisplacement {
            child, amplitude, ..
        } => eval_interval(child, bounds).expand(amplitude.abs()),
        SdfNode::PolarRepeat { child, count } => {
            let sector = std::f32::consts::TAU / (*count as f32);
            let max_r = bounds.length_xz().hi;
            let ha = sector * 0.5;
            eval_interval(
                child,
                Vec3Interval {
                    x: Interval::new(0.0, max_r),
                    y: bounds.y,
                    z: Interval::new(-max_r * ha.sin(), max_r * ha.sin()),
                },
            )
        }
        SdfNode::SweepBezier { child, p0, p1, p2 } => {
            let bmin_x = p0.x.min(p1.x).min(p2.x);
            let bmax_x = p0.x.max(p1.x).max(p2.x);
            let bmin_z = p0.y.min(p1.y).min(p2.y);
            let bmax_z = p0.y.max(p1.y).max(p2.y);
            // Max perpendicular distance
            let corners = [
                (bounds.x.lo, bounds.z.lo),
                (bounds.x.lo, bounds.z.hi),
                (bounds.x.hi, bounds.z.lo),
                (bounds.x.hi, bounds.z.hi),
            ];
            let curve_corners = [
                (bmin_x, bmin_z),
                (bmin_x, bmax_z),
                (bmax_x, bmin_z),
                (bmax_x, bmax_z),
            ];
            let mut max_d2: f32 = 0.0;
            for &(px, pz) in &corners {
                for &(cx, cz) in &curve_corners {
                    let dx = px - cx;
                    let dz = pz - cz;
                    max_d2 = max_d2.max(dx.mul_add(dx, dz * dz));
                }
            }
            // Min distance: 0 if overlapping, else min box-to-box distance
            let min_perp = {
                let dx = if bounds.x.hi < bmin_x {
                    bmin_x - bounds.x.hi
                } else if bounds.x.lo > bmax_x {
                    bounds.x.lo - bmax_x
                } else {
                    0.0
                };
                let dz = if bounds.z.hi < bmin_z {
                    bmin_z - bounds.z.hi
                } else if bounds.z.lo > bmax_z {
                    bounds.z.lo - bmax_z
                } else {
                    0.0
                };
                dx.hypot(dz)
            };
            eval_interval(
                child,
                Vec3Interval {
                    x: Interval::new(min_perp, max_d2.sqrt()),
                    y: bounds.y,
                    z: Interval::ZERO,
                },
            )
        }

        // The fold is a sequence of reflections (norm-preserving), so the image lies
        // in the cube of half-size |p|.hi; the initial abs keeps the lower corner ≥ 0
        // is not guaranteed after the plane reflections, hence the full cube.
        SdfNode::IcosahedralSymmetry { child } => {
            let l = bounds.length().hi;
            let cube = Interval::new(-l, l);
            eval_interval(
                child,
                Vec3Interval {
                    x: cube,
                    y: cube,
                    z: cube,
                },
            )
        }
        SdfNode::IFS { child, .. } => {
            // Conservative: IFS may contract, so interval is child's interval
            eval_interval(child, bounds)
        }
        SdfNode::HeightmapDisplacement {
            child, amplitude, ..
        } => {
            let child_iv = eval_interval(child, bounds);
            // Displacement shifts distance by at most amplitude
            let amp = *amplitude;
            Interval::new(child_iv.lo - amp.abs(), child_iv.hi + amp.abs())
        }
        // |fbm| ≤ Σ_{i<octaves} 0.5^i = 2 - 2^(1 - octaves) (`modifiers::surface_roughness::fbm_bound`)
        SdfNode::SurfaceRoughness {
            child,
            amplitude,
            octaves,
            ..
        } => {
            let child_iv = eval_interval(child, bounds);
            let amp = amplitude.abs() * crate::modifiers::fbm_bound(*octaves);
            Interval::new(child_iv.lo - amp, child_iv.hi + amp)
        }

        SdfNode::WithMaterial { child, .. } => eval_interval(child, bounds),

        #[allow(unreachable_patterns)]
        _ => Interval::EVERYTHING,
    }
}

/// Box centre and half-diagonal (radius of the box's circumscribed sphere).
#[inline]
fn centre_and_half_diag(bounds: Vec3Interval) -> (Vec3, f32) {
    let c = Vec3::new(
        bounds.x.midpoint(),
        bounds.y.midpoint(),
        bounds.z.midpoint(),
    );
    let h = Vec3::new(bounds.x.width(), bounds.y.width(), bounds.z.width()) * 0.5;
    (c, h.length())
}

/// Lipschitz interval: `eval(centre) ± l · half_diagonal`.
///
/// Sound for every node whose point evaluation is `l`-Lipschitz on the box,
/// independent of the shape's extent or of whether the evaluation is an
/// exact distance or a lower bound (exact SDFs, `max` of unit-normal plane
/// distances, extrusions of 2D SDFs are all 1-Lipschitz). Preferred over
/// [`ia_bsphere`] for primitives whose bounding radius is not trivially
/// derivable from the parameters.
#[inline]
fn ia_lipschitz(node: &SdfNode, bounds: Vec3Interval, l: f32) -> Interval {
    let (c, rho) = centre_and_half_diag(bounds);
    let d = crate::eval::eval(node, c);
    if !d.is_finite() {
        return Interval::EVERYTHING;
    }
    Interval::new(d - l * rho, d + l * rho)
}

/// `x.powf(m)` on a non-negative interval (monotone for `m > 0`).
#[inline]
fn powf_nonneg(iv: Interval, m: f32) -> Interval {
    let lo = iv.lo.max(0.0);
    let hi = iv.hi.max(0.0);
    Interval::new(lo.powf(m), hi.powf(m))
}

/// 2D length of two intervals.
#[inline]
fn len2(a: Interval, b: Interval) -> Interval {
    (a.sqr() + b.sqr()).sqrt()
}

/// Bounding sphere conservative interval: `length(p) - radius`
///
/// Sound only when the shape lies inside the sphere of `radius` **and** the
/// point evaluation never underestimates the Euclidean distance.
#[inline]
fn ia_bsphere(bounds: Vec3Interval, radius: f32) -> Interval {
    let l = bounds.length();
    Interval::new(l.lo - radius, l.hi + radius)
}

// ============================================================
// Lipschitz bound tracking
// ============================================================

/// Compute an upper bound on the Lipschitz constant of the SDF.
///
/// `L` bounds every difference quotient of the field on the closed exterior
/// `{f ≥ 0}`: `|f(p) − f(q)| ≤ L·|p − q|` whenever `f(p) ≥ 0` or `f(q) ≥ 0`.
/// That is exactly what sphere tracing needs — a step of `d / L` from a point
/// with `d = f(p) > 0` cannot cross the surface — and it lets a field be
/// discontinuous strictly inside the solid (the IQ ellipsoid at its centre)
/// without poisoning the bound.
///
/// Every bound here is either analytic (exact SDFs: 1; boolean ops: the
/// hg_sdf forms below; barrel / twist / bend: the larger singular value of
/// the shear Jacobian) or a numerically established supremum pinned with a
/// margin (TPMS implicit functions, Perlin noise). Variants whose law is
/// **not** Lipschitz on the exterior — domain repetition with an arbitrary
/// child, the taper singular plane, laws with a candidate-selection jump —
/// return [`f32::INFINITY`] rather than a guess: `RaymarchConfig::relaxed`
/// then falls back to plain sphere tracing. The claim is checked for every
/// corpus node by `tests/test_evaluator_opcode_parity.rs::lipschitz_claim_bounds_every_difference_quotient`.
///
/// Until 1.10.3 TPMS surfaces and the ellipsoid claimed `L = 1` (measured
/// 1.7–7.0×), so `RaymarchConfig::relaxed` skipped their surface: Neovius
/// and IWP could not be rendered at all (external review, 2026-09-15).
pub fn eval_lipschitz(node: &SdfNode) -> f32 {
    match node {
        // Exact SDFs (|∇f| = 1 almost everywhere, L = 1)
        SdfNode::Sphere { .. }
        | SdfNode::Box3d { .. }
        | SdfNode::Cylinder { .. }
        | SdfNode::Torus { .. }
        | SdfNode::Plane { .. }
        | SdfNode::Capsule { .. }
        | SdfNode::Cone { .. }
        | SdfNode::RoundedCone { .. }
        | SdfNode::Pyramid { .. }
        | SdfNode::Octahedron { .. }
        | SdfNode::HexPrism { .. }
        | SdfNode::Link { .. }
        | SdfNode::Triangle { .. }
        | SdfNode::Bezier { .. }
        | SdfNode::RoundedBox { .. }
        | SdfNode::CappedCone { .. }
        | SdfNode::CappedTorus { .. }
        | SdfNode::RoundedCylinder { .. }
        | SdfNode::TriangularPrism { .. }
        | SdfNode::CutSphere { .. }
        | SdfNode::CutHollowSphere { .. }
        | SdfNode::DeathStar { .. }
        | SdfNode::SolidAngle { .. }
        | SdfNode::Rhombus { .. }
        | SdfNode::Vesica { .. }
        | SdfNode::InfiniteCylinder { .. }
        | SdfNode::InfiniteCone { .. }
        | SdfNode::Heart { .. }
        | SdfNode::Tube { .. }
        | SdfNode::Diamond { .. }
        | SdfNode::ChamferedCube { .. }
        | SdfNode::Superellipsoid { .. }
        | SdfNode::RoundedX { .. }
        | SdfNode::Pie { .. }
        | SdfNode::Trapezoid { .. }
        | SdfNode::Parallelogram { .. }
        | SdfNode::Tunnel { .. }
        | SdfNode::UnevenCapsule { .. }
        | SdfNode::ArcShape { .. }
        | SdfNode::Moon { .. }
        | SdfNode::CrossShape { .. }
        | SdfNode::ParabolaSegment { .. }
        | SdfNode::RegularPolygon { .. }
        | SdfNode::StarPolygon { .. }
        | SdfNode::Tetrahedron { .. }
        | SdfNode::Dodecahedron { .. }
        | SdfNode::Icosahedron { .. }
        | SdfNode::TruncatedOctahedron { .. }
        | SdfNode::TruncatedIcosahedron { .. }
        | SdfNode::BoxFrame { .. }
        | SdfNode::Circle2D { .. }
        | SdfNode::Rect2D { .. }
        | SdfNode::Segment2D { .. }
        | SdfNode::Polygon2D { .. }
        | SdfNode::RoundedRect2D { .. }
        | SdfNode::Annular2D { .. }
        | SdfNode::Egg { .. }
        | SdfNode::Horseshoe { .. }
        | SdfNode::BlobbyCross { .. }
        | SdfNode::Stairs { .. }
        | SdfNode::Helix { .. }
        | SdfNode::Ellipsoid { .. } => 1.0,

        // Triply periodic minimal surfaces: `|F(p·s)| / s − t` with F an
        // implicit trigonometric function, so L = sup|∇F| independent of
        // `scale` (the chain rule cancels it). Suprema located numerically
        // (200k samples + hill climb over ±7, two scales) and pinned with
        // the property test; the closed forms match to 4 digits.
        SdfNode::Gyroid { .. } | SdfNode::SchwarzP { .. } | SdfNode::DiamondSurface { .. } => {
            TPMS_SQRT3
        }
        SdfNode::Neovius { .. } => TPMS_NEOVIUS,
        SdfNode::Lidinoid { .. } => TPMS_LIDINOID,
        SdfNode::IWP { .. } => TPMS_IWP,
        SdfNode::FRD { .. } | SdfNode::FischerKochS { .. } => TPMS_FRD,
        SdfNode::PMY { .. } => TPMS_PMY,

        // Barrel: `r − (radius + bulge·(1 − (y/h)²))` is a shear of the
        // radial coordinate by `k = 2·bulge/h`; the box-corner combination
        // with `|y| − h` has the singular value of [[1, 0], [k, 1]].
        SdfNode::Barrel {
            half_height, bulge, ..
        } => shear_singular_value(2.0 * bulge.abs() / half_height.abs().max(1e-6)),

        // Terrain: `y − amplitude·fbm(scale·xz)`, 3 octaves of 2-D value
        // noise (values in [0, 1], smoothstep weights, gain 0.48, lacunarity
        // 2.1 with a rotation). |∇fbm| ≤ Σ aᵢ·2.1ⁱ · 1.5·√2 = 3.208.
        SdfNode::Terrain { scale, amplitude } => {
            TERRAIN_FBM_GRAD.mul_add((scale * amplitude).abs(), 1.0)
        }

        // min / max and every convex blend (smooth, exp-smooth: the weights
        // on ∇a and ∇b sum to 1) are 1-Lipschitz in (a, b).
        SdfNode::Union { a, b }
        | SdfNode::Intersection { a, b }
        | SdfNode::Subtraction { a, b }
        | SdfNode::SmoothUnion { a, b, .. }
        | SdfNode::SmoothIntersection { a, b, .. }
        | SdfNode::SmoothSubtraction { a, b, .. }
        | SdfNode::XOR { a, b }
        | SdfNode::Groove { a, b, .. }
        | SdfNode::Tongue { a, b, .. }
        | SdfNode::ExpSmoothUnion { a, b, .. }
        | SdfNode::ExpSmoothIntersection { a, b, .. }
        | SdfNode::ExpSmoothSubtraction { a, b, .. } => eval_lipschitz(a).max(eval_lipschitz(b)),
        // Chamfer `(a + b)/√2 − r`, stairs (its 45° branch), engrave
        // `(a + r − |b|)/√2`: |∇a + ∇b| / √2 ≤ √2·max(La, Lb).
        SdfNode::ChamferUnion { a, b, .. }
        | SdfNode::ChamferIntersection { a, b, .. }
        | SdfNode::ChamferSubtraction { a, b, .. }
        | SdfNode::StairsUnion { a, b, .. }
        | SdfNode::StairsIntersection { a, b, .. }
        | SdfNode::StairsSubtraction { a, b, .. }
        | SdfNode::Engrave { a, b, .. } => {
            std::f32::consts::SQRT_2 * eval_lipschitz(a).max(eval_lipschitz(b))
        }
        // Pipe `√(a² + b²) − r`: |(a∇a + b∇b)| / √(a² + b²) ≤ √(La² + Lb²).
        SdfNode::Pipe { a, b, .. } => eval_lipschitz(a).hypot(eval_lipschitz(b)),
        // Morph `a(1 − t) + b·t`.
        SdfNode::Morph { a, b, t } => {
            (1.0 - t).abs() * eval_lipschitz(a) + t.abs() * eval_lipschitz(b)
        }
        // Columns (hg_sdf fOp*Columns): polar modulo of the (a, b) plane
        // — the field jumps between columns, no finite bound.
        SdfNode::ColumnsUnion { .. }
        | SdfNode::ColumnsIntersection { .. }
        | SdfNode::ColumnsSubtraction { .. } => f32::INFINITY,

        // Distance-preserving transforms
        SdfNode::Translate { child, .. } | SdfNode::Rotate { child, .. } => eval_lipschitz(child),
        // Uniform scale: SDF(p/s)*s, gradient = (1/s)*∇child*s = ∇child
        SdfNode::Scale { child, .. } => eval_lipschitz(child),
        // Non-uniform scale: gradient stretches by max/min factor ratio
        SdfNode::ScaleNonUniform { child, factors } => {
            let s = [factors.x.abs(), factors.y.abs(), factors.z.abs()];
            let min_s = s[0].min(s[1]).min(s[2]).max(1e-6);
            let max_s = s[0].max(s[1]).max(s[2]);
            eval_lipschitz(child) * max_s / min_s
        }
        // Projective transform: use user-provided bound
        SdfNode::ProjectiveTransform {
            child,
            lipschitz_bound,
            ..
        } => eval_lipschitz(child) * lipschitz_bound,
        // Lattice deform is the identity outside the lattice box and the
        // deformed field inside: the two disagree on the box faces.
        SdfNode::LatticeDeform { .. } => f32::INFINITY,
        // SDF Skinning: LBS preserves distance (approximately)
        SdfNode::SdfSkinning { child, .. } => eval_lipschitz(child),

        // Modifiers with Jacobian norm ≤ 1
        SdfNode::Round { child, .. }
        | SdfNode::Onion { child, .. }
        | SdfNode::Elongate { child, .. }
        | SdfNode::Mirror { child, .. }
        | SdfNode::OctantMirror { child, .. }
        | SdfNode::Revolution { child, .. }
        | SdfNode::Extrude { child, .. } => eval_lipschitz(child),
        // Sweep along a Bézier: (distance to the curve, y) is 1-Lipschitz in p
        // (closed-form nearest point since 1.11.0).
        SdfNode::SweepBezier { child, .. } => eval_lipschitz(child),

        // Domain repetition evaluates one copy per cell / sector; unless the
        // child is symmetric inside its cell (not knowable here) the field
        // jumps on the cell borders.
        SdfNode::RepeatInfinite { .. }
        | SdfNode::RepeatFinite { .. }
        | SdfNode::PolarRepeat { .. } => f32::INFINITY,

        // Twist about Y: q = R(s·y)·p, ∂q/∂y = s·R'(s·y)·p_xz with |·| = s·r.
        // The Jacobian is a shear by v = s·r whose larger singular value is
        // v/2 + √(1 + v²/4); r is the child's largest XZ radius (its AABB
        // corner), the far-field derivative of a distance field through this
        // map is bounded by the same s·r_max.
        SdfNode::Twist { child, strength } => {
            eval_lipschitz(child)
                * shear_singular_value(strength.abs() * child_radius(child, Axis::Y))
        }
        // Bend about Z: q = R(k·x)·p_xy, ∂q/∂x = R + k·(−q_y, q_x), so
        // ‖J‖ ≤ 1 + k·r with r the child's largest XY radius.
        SdfNode::Bend { child, curvature } => {
            eval_lipschitz(child) * curvature.abs().mul_add(child_radius(child, Axis::Z), 1.0)
        }

        // Displacements add the gradient of the offset field.
        // Perlin: |∇perlin| ≤ PERLIN_GRAD (numerical supremum, 4 seeds).
        SdfNode::Noise {
            child,
            amplitude,
            frequency,
            ..
        } => (amplitude * frequency)
            .abs()
            .mul_add(PERLIN_GRAD, eval_lipschitz(child)),
        // Taper divides XZ by `1 − f·y`, singular on the plane y = 1/f.
        SdfNode::Taper { .. } => f32::INFINITY,
        // `d + s·sin(5x)·sin(5y)·sin(5z)`: |∇(sin·sin·sin)| ≤ 1, times 5.
        SdfNode::Displacement { child, strength } => {
            strength.abs().mul_add(5.0, eval_lipschitz(child))
        }
        // `d + a·sin(fx·x)·sin(fy·y)·sin(fz·z)`: ≤ a·max|f|.
        SdfNode::SineDisplacement {
            child,
            amplitude,
            frequency,
        } => {
            let f = frequency
                .x
                .abs()
                .max(frequency.y.abs())
                .max(frequency.z.abs());
            (amplitude * f).abs() + eval_lipschitz(child)
        }
        SdfNode::Shear { child, shear } => {
            let max_shear = shear.x.abs().max(shear.y.abs()).max(shear.z.abs());
            eval_lipschitz(child) * (1.0 + max_shear)
        }
        SdfNode::Animated { child, .. } => eval_lipschitz(child),
        SdfNode::IcosahedralSymmetry { child } => eval_lipschitz(child),
        SdfNode::IFS { child, .. } => {
            // IFS can have arbitrary scale factors, conservative estimate
            eval_lipschitz(child) * 2.0
        }
        // Heightmap displacement projects on the dominant axis (a jump where
        // the axis changes) and samples texels bilinearly (slope ∝ resolution).
        SdfNode::HeightmapDisplacement { .. } => f32::INFINITY,
        // fbm of value noise: each octave contributes 0.5ⁱ·2ⁱ·|∇hash_noise|
        // with |∇hash_noise_3d| ≤ 2·1.5·√3 (range 2, smoothstep slope 1.5).
        SdfNode::SurfaceRoughness {
            child,
            amplitude,
            frequency,
            octaves,
        } => (amplitude * frequency * *octaves as f32)
            .abs()
            .mul_add(VALUE_NOISE_GRAD, eval_lipschitz(child)),

        SdfNode::WithMaterial { child, .. } => eval_lipschitz(child),
    }
}

/// `√3`: Gyroid, Schwarz P and the diamond surface (`sup|∇F|`, attained at
/// the lattice points where every trig factor is ±1).
const TPMS_SQRT3: f32 = 1.732_051;
/// Neovius `3(cx + cy + cz) + 4·cx·cy·cz`: sup 7.0006 → 7.05.
const TPMS_NEOVIUS: f32 = 7.05;
/// Lidinoid: sup 2.5980 = 3√3/2 → 2.62.
const TPMS_LIDINOID: f32 = 2.62;
/// IWP: sup 5.1966 = 3√3 → 5.25.
const TPMS_IWP: f32 = 5.25;
/// FRD and Fischer–Koch S: sup 2.4439 → 2.47.
const TPMS_FRD: f32 = 2.47;
/// PMY: sup 4.1232 = √17 → 4.17.
const TPMS_PMY: f32 = 4.17;
/// `|∇perlin_noise_3d|` supremum: 3.266 over 4 seeds × 300k samples + hill
/// climb (the analytic bound through the quintic fade, 17.6, is 5× looser).
const PERLIN_GRAD: f32 = 3.5;
/// `|∇hash_noise_3d|` analytic bound: per axis 2 (range) × 1.5 (smoothstep
/// slope), three axes → 3√3. Measured 2.95.
const VALUE_NOISE_GRAD: f32 = 5.196_153;
/// Terrain fbm gradient factor: `(0.5 + 0.5·0.48·2.1 + 0.5·0.48²·2.1²) · 1.5 · √2`.
const TERRAIN_FBM_GRAD: f32 = 3.208;
/// Radius assumed for a twisted / bent child whose extent is unbounded
/// (plane, infinite cylinder, uncompilable tree): the pre-1.10.3 heuristic.
const UNBOUNDED_CHILD_RADIUS: f32 = 10.0;
/// AABB reach beyond which a child counts as unbounded (`AabbPacked::infinite()`
/// uses ±`f32::MAX`; real scenes stay far below 1e6 units).
const UNBOUNDED_CHILD_SENTINEL: f32 = 1e6;

/// Larger singular value of the shear `[[1, 0], [v, 1]]`: `v/2 + √(1 + v²/4)`.
#[inline]
fn shear_singular_value(v: f32) -> f32 {
    (v * 0.25).mul_add(v, 1.0).sqrt() + 0.5 * v
}

/// Axis a deformation rotates about; the radius is measured in the
/// perpendicular plane.
#[derive(Clone, Copy)]
enum Axis {
    Y,
    Z,
}

/// Largest distance from `axis` reached by the child's bounding box
/// (twist / bend Jacobian growth is proportional to it). Falls back to
/// [`UNBOUNDED_CHILD_RADIUS`] when the child has no finite AABB.
fn child_radius(child: &SdfNode, axis: Axis) -> f32 {
    use crate::compiled::{get_scene_aabb, CompiledSdfBvh};
    let Ok(bvh) = CompiledSdfBvh::try_compile(child) else {
        return UNBOUNDED_CHILD_RADIUS;
    };
    let aabb = get_scene_aabb(&bvh);
    let (lo, hi) = (aabb.min(), aabb.max());
    // `AabbPacked::infinite()` is ±f32::MAX (finite): treat anything past the
    // unbounded-shape sentinel as unbounded rather than overflowing to inf.
    if !aabb.is_valid() || !lo.is_finite() || !hi.is_finite() {
        return UNBOUNDED_CHILD_RADIUS;
    }
    let reach = |a: f32, b: f32| a.abs().max(b.abs());
    let r = match axis {
        Axis::Y => reach(lo.x, hi.x).hypot(reach(lo.z, hi.z)),
        Axis::Z => reach(lo.x, hi.x).hypot(reach(lo.y, hi.y)),
    };
    if r > UNBOUNDED_CHILD_SENTINEL {
        UNBOUNDED_CHILD_RADIUS
    } else {
        r
    }
}

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

    /// Verify that scalar eval at sampled points falls within interval bounds
    fn check(node: &SdfNode, min: Vec3, max: Vec3, n: usize) {
        let iv = eval_interval(node, Vec3Interval::from_bounds(min, max));
        let step = (max - min) / ((n - 1) as f32);
        for ix in 0..n {
            for iy in 0..n {
                for iz in 0..n {
                    let p = min + Vec3::new(ix as f32, iy as f32, iz as f32) * step;
                    let d = eval(node, p);
                    assert!(
                        d >= iv.lo - 1e-4 && d <= iv.hi + 1e-4,
                        "d={d} not in [{lo}, {hi}] at {p:?}",
                        lo = iv.lo,
                        hi = iv.hi,
                    );
                }
            }
        }
    }

    #[test]
    fn test_interval_ops() {
        let a = Interval::new(1.0, 3.0);
        let b = Interval::new(2.0, 5.0);
        assert_eq!((a + b).lo, 3.0);
        assert_eq!((a + b).hi, 8.0);
        assert_eq!((a - b).lo, -4.0);
        assert_eq!((a - b).hi, 1.0);
        assert_eq!((-a).lo, -3.0);
        assert_eq!((-a).hi, -1.0);
    }

    #[test]
    fn test_interval_abs() {
        assert_eq!(Interval::new(1.0, 3.0).abs().lo, 1.0);
        assert_eq!(Interval::new(-3.0, -1.0).abs().hi, 3.0);
        assert_eq!(Interval::new(-2.0, 3.0).abs().lo, 0.0);
    }

    #[test]
    fn test_interval_sqr() {
        let s = Interval::new(-2.0, 3.0).sqr();
        assert_eq!(s.lo, 0.0);
        assert_eq!(s.hi, 9.0);
    }

    #[test]
    fn test_sphere() {
        check(
            &SdfNode::sphere(1.0),
            Vec3::splat(-2.0),
            Vec3::splat(2.0),
            5,
        );
    }

    #[test]
    fn test_box3d() {
        check(
            &SdfNode::box3d(1.0, 0.5, 0.75),
            Vec3::splat(-2.0),
            Vec3::splat(2.0),
            5,
        );
    }

    #[test]
    fn test_cylinder() {
        check(
            &SdfNode::cylinder(0.5, 1.0),
            Vec3::splat(-2.0),
            Vec3::splat(2.0),
            5,
        );
    }

    #[test]
    fn test_torus() {
        check(
            &SdfNode::torus(1.0, 0.25),
            Vec3::splat(-2.0),
            Vec3::splat(2.0),
            5,
        );
    }

    #[test]
    fn test_plane() {
        check(
            &SdfNode::plane(Vec3::Y, 0.0),
            Vec3::splat(-2.0),
            Vec3::splat(2.0),
            5,
        );
    }

    #[test]
    fn test_rounded_box() {
        check(
            &SdfNode::rounded_box(1.0, 0.5, 0.75, 0.1),
            Vec3::splat(-2.0),
            Vec3::splat(2.0),
            4,
        );
    }

    #[test]
    fn test_union() {
        let node = SdfNode::sphere(1.0)
            .translate(1.0, 0.0, 0.0)
            .union(SdfNode::sphere(1.0).translate(-1.0, 0.0, 0.0));
        check(&node, Vec3::splat(-3.0), Vec3::splat(3.0), 4);
    }

    #[test]
    fn test_smooth_union() {
        let node = SdfNode::sphere(1.0)
            .translate(0.5, 0.0, 0.0)
            .smooth_union(SdfNode::sphere(1.0).translate(-0.5, 0.0, 0.0), 0.5);
        check(&node, Vec3::splat(-3.0), Vec3::splat(3.0), 4);
    }

    #[test]
    fn test_subtraction() {
        let node = SdfNode::sphere(1.5).subtract(SdfNode::sphere(1.0).translate(0.5, 0.0, 0.0));
        check(&node, Vec3::splat(-3.0), Vec3::splat(3.0), 4);
    }

    #[test]
    fn test_translate() {
        check(
            &SdfNode::sphere(1.0).translate(2.0, 0.0, 0.0),
            Vec3::new(0.0, -2.0, -2.0),
            Vec3::new(4.0, 2.0, 2.0),
            4,
        );
    }

    #[test]
    fn test_rotate() {
        let node = SdfNode::box3d(1.0, 0.5, 0.3).rotate(Quat::from_rotation_y(0.7));
        check(&node, Vec3::splat(-2.0), Vec3::splat(2.0), 4);
    }

    #[test]
    fn test_scale() {
        check(
            &SdfNode::sphere(1.0).scale(2.0),
            Vec3::splat(-4.0),
            Vec3::splat(4.0),
            4,
        );
    }

    #[test]
    fn test_round() {
        let node = SdfNode::box3d(1.0, 1.0, 1.0).round(0.2);
        check(&node, Vec3::splat(-2.0), Vec3::splat(2.0), 4);
    }

    #[test]
    fn test_onion() {
        let node = SdfNode::sphere(1.0).onion(0.1);
        check(&node, Vec3::splat(-2.0), Vec3::splat(2.0), 4);
    }

    #[test]
    fn test_revolution() {
        let node = SdfNode::Revolution {
            child: std::sync::Arc::new(SdfNode::sphere(0.3)),
            offset: 1.0,
        };
        check(&node, Vec3::splat(-2.0), Vec3::splat(2.0), 4);
    }

    #[test]
    fn test_extrude() {
        let node = SdfNode::Extrude {
            child: std::sync::Arc::new(SdfNode::sphere(0.5)),
            half_height: 1.0,
        };
        check(&node, Vec3::splat(-2.0), Vec3::splat(2.0), 4);
    }

    #[test]
    fn test_mirror() {
        let node = SdfNode::box3d(1.0, 0.5, 0.5).mirror(true, false, false);
        check(&node, Vec3::splat(-2.0), Vec3::splat(2.0), 4);
    }

    #[test]
    fn test_pruning() {
        let sphere = SdfNode::sphere(1.0);
        // Far region: entirely positive
        let far = Vec3Interval::from_bounds(Vec3::splat(5.0), Vec3::splat(6.0));
        assert!(eval_interval(&sphere, far).is_positive());
        // Inside region: entirely negative
        let inside = Vec3Interval::from_bounds(Vec3::splat(-0.1), Vec3::splat(0.1));
        assert!(eval_interval(&sphere, inside).is_negative());
    }

    #[test]
    fn test_chamfer_blend() {
        let node = SdfNode::sphere(1.0).chamfer_union(SdfNode::box3d(0.8, 0.8, 0.8), 0.3);
        check(&node, Vec3::splat(-2.0), Vec3::splat(2.0), 4);
    }

    #[test]
    fn test_stairs_blend() {
        let node = SdfNode::sphere(1.0).stairs_union(SdfNode::box3d(0.8, 0.8, 0.8), 0.3, 4.0);
        check(&node, Vec3::splat(-2.0), Vec3::splat(2.0), 4);
    }

    #[test]
    fn test_conservative_primitives() {
        // These use bounding sphere — should still contain all scalar values
        check(
            &SdfNode::cone(0.5, 1.0),
            Vec3::splat(-3.0),
            Vec3::splat(3.0),
            4,
        );
        check(
            &SdfNode::octahedron(1.0),
            Vec3::splat(-2.0),
            Vec3::splat(2.0),
            4,
        );
        check(
            &SdfNode::pyramid(1.0),
            Vec3::splat(-3.0),
            Vec3::splat(3.0),
            4,
        );
    }

    // ============ Lipschitz tests ============

    /// Check that finite-difference gradient ≤ L at sampled points
    fn check_lipschitz(node: &SdfNode, min: Vec3, max: Vec3, n: usize) {
        let lip = eval_lipschitz(node);
        assert!(lip >= 1.0, "Lipschitz should be >= 1.0, got {lip}");
        let eps = 1e-3;
        let step = (max - min) / ((n - 1) as f32);
        let mut max_grad: f32 = 0.0;
        for ix in 0..n {
            for iy in 0..n {
                for iz in 0..n {
                    let p = min + step * Vec3::new(ix as f32, iy as f32, iz as f32);
                    let d = eval(node, p);
                    // Finite-difference gradient (central differences)
                    let gx = (eval(node, p + Vec3::X * eps) - eval(node, p - Vec3::X * eps))
                        / (2.0 * eps);
                    let gy = (eval(node, p + Vec3::Y * eps) - eval(node, p - Vec3::Y * eps))
                        / (2.0 * eps);
                    let gz = (eval(node, p + Vec3::Z * eps) - eval(node, p - Vec3::Z * eps))
                        / (2.0 * eps);
                    let grad_mag = (gx * gx + gy * gy + gz * gz).sqrt();
                    max_grad = max_grad.max(grad_mag);
                    // Allow 5% tolerance for finite-difference error
                    assert!(
                        grad_mag <= lip.mul_add(1.05, 0.01),
                        "|∇SDF|={grad_mag:.4} > L={lip:.4} at {p:?} (d={d:.4})"
                    );
                }
            }
        }
    }

    #[test]
    fn test_lipschitz_sphere() {
        check_lipschitz(
            &SdfNode::sphere(1.0),
            Vec3::splat(-3.0),
            Vec3::splat(3.0),
            6,
        );
    }

    #[test]
    fn test_lipschitz_box() {
        check_lipschitz(
            &SdfNode::box3d(1.0, 0.5, 0.75),
            Vec3::splat(-3.0),
            Vec3::splat(3.0),
            6,
        );
    }

    #[test]
    fn test_lipschitz_union() {
        let node = SdfNode::sphere(1.0)
            .translate(1.0, 0.0, 0.0)
            .union(SdfNode::sphere(1.0).translate(-1.0, 0.0, 0.0));
        assert_eq!(eval_lipschitz(&node), 1.0);
        check_lipschitz(&node, Vec3::splat(-3.0), Vec3::splat(3.0), 5);
    }

    #[test]
    fn test_lipschitz_smooth_union() {
        let node = SdfNode::sphere(1.0)
            .translate(0.5, 0.0, 0.0)
            .smooth_union(SdfNode::sphere(1.0).translate(-0.5, 0.0, 0.0), 0.5);
        assert_eq!(eval_lipschitz(&node), 1.0);
        check_lipschitz(&node, Vec3::splat(-3.0), Vec3::splat(3.0), 5);
    }

    #[test]
    fn test_lipschitz_scale() {
        let node = SdfNode::sphere(1.0).scale(2.0);
        assert_eq!(eval_lipschitz(&node), 1.0);
    }

    #[test]
    fn test_lipschitz_non_uniform_scale() {
        let node = SdfNode::ScaleNonUniform {
            child: std::sync::Arc::new(SdfNode::sphere(1.0)),
            factors: Vec3::new(2.0, 1.0, 1.0),
        };
        assert_eq!(eval_lipschitz(&node), 2.0);
    }

    #[test]
    fn test_lipschitz_noise() {
        let node = SdfNode::Noise {
            child: std::sync::Arc::new(SdfNode::sphere(1.0)),
            amplitude: 0.2,
            frequency: 3.0,
            seed: 0,
        };
        // L = 1.0 + 0.2 * 3.0 * PERLIN_GRAD (|∇perlin| ≤ 3.5)
        assert!((eval_lipschitz(&node) - (1.0 + 0.6 * PERLIN_GRAD)).abs() < 1e-5);
    }

    #[test]
    fn test_lipschitz_twist() {
        let node = SdfNode::Twist {
            child: std::sync::Arc::new(SdfNode::sphere(1.0)),
            strength: 0.5,
        };
        // Unit sphere: XZ reach of its AABB corner is √2, v = 0.5·√2, the
        // shear singular value is v/2 + √(1 + v²/4).
        let v = 0.5 * 2.0_f32.sqrt();
        let expected = 0.5 * v + (1.0 + 0.25 * v * v).sqrt();
        assert!((eval_lipschitz(&node) - expected).abs() < 1e-3);
    }

    #[test]
    fn test_lipschitz_round() {
        let node = SdfNode::box3d(1.0, 1.0, 1.0).round(0.2);
        assert_eq!(eval_lipschitz(&node), 1.0);
        check_lipschitz(&node, Vec3::splat(-3.0), Vec3::splat(3.0), 5);
    }

    #[test]
    fn interval_width() {
        let iv = Interval::new(1.0, 3.0);
        assert!((iv.width() - 2.0).abs() < 1e-6);
    }

    #[test]
    fn interval_midpoint() {
        let iv = Interval::new(1.0, 5.0);
        assert!((iv.midpoint() - 3.0).abs() < 1e-6);
    }

    #[test]
    fn interval_contains_value() {
        let iv = Interval::new(-1.0, 2.0);
        assert!(iv.contains(0.0));
        assert!(iv.contains(-1.0));
        assert!(iv.contains(2.0));
        assert!(!iv.contains(3.0));
        assert!(!iv.contains(-2.0));
    }

    #[test]
    fn interval_overlaps() {
        let a = Interval::new(0.0, 2.0);
        let b = Interval::new(1.0, 3.0);
        let c = Interval::new(3.0, 5.0);
        assert!(a.overlaps(b));
        assert!(!a.overlaps(c));
    }

    #[test]
    fn interval_intersection() {
        let a = Interval::new(0.0, 3.0);
        let b = Interval::new(1.0, 5.0);
        let c = a.intersect(b);
        assert!((c.lo - 1.0).abs() < 1e-6);
        assert!((c.hi - 3.0).abs() < 1e-6);
    }

    #[test]
    fn interval_hull() {
        let a = Interval::new(1.0, 2.0);
        let b = Interval::new(4.0, 5.0);
        let h = a.hull(b);
        assert!((h.lo - 1.0).abs() < 1e-6);
        assert!((h.hi - 5.0).abs() < 1e-6);
    }

    #[test]
    fn test_lipschitz_subtraction() {
        let node = SdfNode::sphere(1.0).subtract(SdfNode::box3d(0.5, 0.5, 0.5));
        assert_eq!(eval_lipschitz(&node), 1.0);
    }

    #[test]
    fn test_lipschitz_intersection() {
        let node = SdfNode::sphere(1.0).intersection(SdfNode::box3d(0.5, 0.5, 0.5));
        assert_eq!(eval_lipschitz(&node), 1.0);
    }

    #[test]
    fn test_lipschitz_onion() {
        let node = SdfNode::Onion {
            child: std::sync::Arc::new(SdfNode::sphere(1.0)),
            thickness: 0.1,
        };
        assert_eq!(eval_lipschitz(&node), 1.0);
    }

    #[test]
    fn test_interval_morph() {
        let node = SdfNode::Morph {
            a: std::sync::Arc::new(SdfNode::sphere(1.0)),
            b: std::sync::Arc::new(SdfNode::box3d(1.0, 1.0, 1.0)),
            t: 0.5,
        };
        check(&node, Vec3::splat(-3.0), Vec3::splat(3.0), 4);
    }
}