alice-sdf 5.0.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,
}

/// The next representable `f32` below `x` (toward −∞).
///
/// Rust has no stable rounding-mode control, so the outward rounding below is
/// done by nudging the result of a round-to-nearest operation by one ulp. That
/// is sound because a single IEEE-754 operation is within half an ulp of the
/// exact result. `next_down`/`next_up` in std are newer than this crate's MSRV,
/// hence the bit arithmetic.
#[inline(always)]
#[must_use]
pub fn next_down(x: f32) -> f32 {
    if x.is_nan() || x == f32::NEG_INFINITY {
        return x;
    }
    if x == 0.0 {
        // both zeros step to the smallest negative subnormal
        return -f32::from_bits(1);
    }
    let bits = x.to_bits();
    f32::from_bits(if x > 0.0 { bits - 1 } else { bits + 1 })
}

/// The next representable `f32` above `x` (toward +∞).
#[inline(always)]
#[must_use]
pub fn next_up(x: f32) -> f32 {
    if x.is_nan() || x == f32::INFINITY {
        return x;
    }
    if x == 0.0 {
        return f32::from_bits(1);
    }
    let bits = x.to_bits();
    f32::from_bits(if x > 0.0 { bits + 1 } else { bits - 1 })
}

impl Interval {
    /// The whole real line `(-inf, +inf)`: "nothing is known about the value".
    ///
    /// The only sound enclosure of a field that is NaN or has no finite bound (a scale of 0 / NaN /
    /// infinity). It is built directly (no outward rounding, no `lo <= hi` assertion needed).
    #[inline(always)]
    #[must_use]
    pub const fn entire() -> Self {
        // ALLOW-RAW-INTERVAL: the infinities are already the widest possible bounds
        Self {
            lo: f32::NEG_INFINITY,
            hi: f32::INFINITY,
        }
    }

    /// Create a new interval from computed bounds.
    ///
    /// The bounds are rounded **outward** by one ulp ([`Interval::outward`]).
    /// This is the safe default: a call site that was not audited still widens
    /// rather than narrows, so an unaudited site costs an ulp of tightness
    /// instead of soundness. Operations that only *select* existing values
    /// (`min` / `max` / `hull` / `clamp` / `abs` / `neg`) build `Self { .. }`
    /// directly and stay exact.
    #[inline(always)]
    pub fn new(lo: f32, hi: f32) -> Self {
        debug_assert!(lo <= hi + 1e-6, "lo ({}) > hi ({})", lo, hi);
        Self::outward(lo, hi)
    }

    /// Create an interval from *computed* bounds, rounded **outward** by one
    /// ulp on each side.
    ///
    /// Every operation that computes a new value must go through this. Without
    /// it the enclosure is only as good as round-to-nearest, and the error
    /// amplifies: the inclusion drift measured over the corpus before this
    /// landed was 2.086·10⁻⁴ relative (`stairs_intersection`, where a child
    /// box's interval escaped its parent's) — roughly a thousand ulps, because
    /// the stairs blend divides by small numbers and takes `floor`.
    ///
    /// Selection-only operations (`min`, `max`, `hull`, `clamp`, `abs`, `neg`)
    /// return existing values unchanged and must **not** be widened.
    ///
    /// What this does *not* cover: the rounding of the point evaluator itself.
    /// Arms built on a centre sample (see `ia_lipschitz`) add an explicit
    /// margin for that instead.
    #[inline(always)]
    #[must_use]
    pub fn outward(lo: f32, hi: f32) -> Self {
        // ALLOW-RAW-INTERVAL: this *is* the widening; the ulp step is applied here
        Self {
            lo: next_down(lo),
            hi: next_up(hi),
        }
    }

    /// Create a point interval [v, v]
    #[inline(always)]
    pub const fn point(v: f32) -> Self {
        // ALLOW-RAW-INTERVAL: [v, v] is exact — v is the value, not a rounding of one
        Self { lo: v, hi: v }
    }

    /// The entire real line
    // ALLOW-RAW-INTERVAL: the infinities are already the widest possible bounds
    pub const EVERYTHING: Self = Self {
        lo: f32::NEG_INFINITY,
        hi: f32::INFINITY,
    };

    /// Zero interval
    // ALLOW-RAW-INTERVAL: 0 is exact in f32; widening it would lose the sign test
    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 const fn midpoint(self) -> f32 {
        f32::midpoint(self.lo, self.hi)
    }

    /// 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 const fn intersect(self, other: Self) -> Self {
        // ALLOW-RAW-INTERVAL: selection only — both bounds are existing endpoints
        Self {
            lo: self.lo.max(other.lo),
            hi: self.hi.min(other.hi),
        }
    }

    /// Hull (union) of two intervals (widest encompassing)
    #[inline(always)]
    pub const fn hull(self, other: Self) -> Self {
        // ALLOW-RAW-INTERVAL: selection only — both bounds are existing endpoints
        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 {
            // ALLOW-RAW-INTERVAL: f32 negation is exact, so the swapped endpoints are too
            Self {
                lo: -self.hi,
                hi: -self.lo,
            }
        } else {
            // ALLOW-RAW-INTERVAL: 0 is an exact lower bound for |x| over a straddling
            // interval, and the upper bound selects an exactly negated endpoint
            Self {
                lo: 0.0,
                hi: self.hi.max(-self.lo),
            }
        }
    }

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

    /// Square of an interval
    #[inline(always)]
    pub fn sqr(self) -> Self {
        if self.lo >= 0.0 {
            Self::outward(self.lo * self.lo, self.hi * self.hi)
        } else if self.hi <= 0.0 {
            Self::outward(self.hi * self.hi, self.lo * self.lo)
        } else {
            // ALLOW-RAW-INTERVAL: 0 is an exact lower bound (a square is non-negative),
            // so only the computed upper bound is widened
            Self {
                lo: 0.0,
                hi: next_up((self.lo * self.lo).max(self.hi * self.hi)),
            }
        }
    }

    /// Minimum of two intervals
    #[inline(always)]
    pub const fn min(self, other: Self) -> Self {
        // ALLOW-RAW-INTERVAL: selection only — both bounds are existing endpoints
        Self {
            lo: self.lo.min(other.lo),
            hi: self.hi.min(other.hi),
        }
    }

    /// Maximum of two intervals
    #[inline(always)]
    pub const fn max(self, other: Self) -> Self {
        // ALLOW-RAW-INTERVAL: selection only — both bounds are existing endpoints
        Self {
            lo: self.lo.max(other.lo),
            hi: self.hi.max(other.hi),
        }
    }

    /// Clamp to scalar range
    #[inline(always)]
    pub const fn clamp(self, lo: f32, hi: f32) -> Self {
        // ALLOW-RAW-INTERVAL: selection only — each bound is either an existing
        // endpoint or one of the caller's exact limits
        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::outward(self.lo - amount, self.hi + amount)
    }
}

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

impl Sub for Interval {
    type Output = Self;
    #[inline(always)]
    fn sub(self, rhs: Self) -> Self {
        Self::outward(self.lo - rhs.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::outward(a.min(b).min(c).min(d), a.max(b).max(c).max(d))
    }
}

impl Neg for Interval {
    type Output = Self;
    #[inline(always)]
    fn neg(self) -> Self {
        // ALLOW-RAW-INTERVAL: f32 negation is exact, so the swapped endpoints are too
        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::outward(self.lo * rhs, self.hi * rhs)
        } else {
            Self::outward(self.hi * rhs, self.lo * rhs)
        }
    }
}

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

    /// Image of the box under `q = A p + b` (`a` is the 3×3 in glam column
    /// order: a[0..3] = column x, …): bounding box of the 8 transformed corners.
    pub fn affine(self, a: &[f32; 9], b: Vec3) -> Self {
        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 map = |c: Vec3| {
            Vec3::new(
                a[6].mul_add(c.z, a[3].mul_add(c.y, a[0] * c.x)) + b.x,
                a[7].mul_add(c.z, a[4].mul_add(c.y, a[1] * c.x)) + b.y,
                a[8].mul_add(c.z, a[5].mul_add(c.y, a[2] * c.x)) + b.z,
            )
        };
        let first = map(corners[0]);
        let mut min = first;
        let mut max = first;
        for &c in &corners[1..] {
            let r = map(c);
            min = min.min(r);
            max = max.max(r);
        }
        Self::from_bounds(min, max)
    }

    /// Smallest box containing both.
    pub const fn hull(&self, other: &Self) -> Self {
        Self {
            x: self.x.hull(other.x),
            y: self.y.hull(other.y),
            z: self.z.hull(other.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::MetricBall {
            radius,
            w_l1,
            w_l2,
            w_linf,
        } => {
            // a non-negative combination of the three norm intervals; each
            // basis norm is monotone in |x|, |y|, |z|
            let ax = bounds.x.abs();
            let ay = bounds.y.abs();
            let az = bounds.z.abs();
            let l1 = ax + ay + az;
            let l2 = bounds.length();
            let linf = ax.max(ay).max(az);
            l1 * Interval::point(*w_l1)
                + l2 * Interval::point(*w_l2)
                + linf * Interval::point(*w_linf)
                - Interval::point(*radius)
        }

        SdfNode::MetricBlend { inner, outer, .. } => {
            // the result is a convex combination of the two fields, so it
            // lies between their pointwise extremes
            let a = eval_interval(inner, bounds);
            let b = eval_interval(outer, bounds);
            Interval::new(a.lo.min(b.lo), a.hi.max(b.hi))
        }

        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),
        // `round_radius` inflates the 2D (x, z) distance, so it grows the
        // bounding sphere too — dropping it made the enclosure too tight by
        // exactly `round_radius` (0.05 in the corpus).
        SdfNode::Rhombus {
            la,
            lb,
            half_height,
            round_radius,
        } => ia_bsphere(bounds, (la.max(*lb) + round_radius).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))
        }
        // The law doubles `cap_height` (`h = cap_height * 2`) and anchors the
        // lower cap at y = 0, so the shape reaches y = h + r2 — the old
        // `max(r1, r2) + cap_height + half_depth` under-bounded that by
        // `cap_height` (0.057 too tight in the corpus).
        SdfNode::UnevenCapsule {
            r1,
            r2,
            cap_height,
            half_depth,
        } => ia_bsphere(
            bounds,
            r1.max(2.0f32.mul_add(*cap_height, *r2)).hypot(*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),
            )
        }
        // The stairs laws are **not** Lipschitz in `(a, b)`: `glsl_mod` jumps by a
        // full tooth at each period boundary, so the old midpoint + half-diagonal
        // form could claim a narrower range than the law actually takes (measured:
        // a child box's interval escaping its parent's by 2.086·10⁻⁴ relative).
        // Evaluate the law op by op on intervals instead ([`ia_stairs_min`]).
        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;
            }
            ia_stairs_min(ia, ib, *r, *n)
        }
        // `stairs_max(a, b) = −stairs_min(−a, −b)`; negation on intervals is exact
        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;
            }
            -ia_stairs_min(-ia, -ib, *r, *n)
        }
        // `stairs_subtraction(a, b) = −stairs_min(−a, b)` (the point law negates the
        // first argument only; the arm below keeps the same shape on intervals)
        SdfNode::StairsSubtraction { 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;
            }
            -ia_stairs_min(-ia, ib, *r, *n)
        }
        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(
                k.max(1e-6).mul_add(-std::f32::consts::LN_2, sharp.lo),
                sharp.hi,
            )
        }
        SdfNode::ExpSmoothIntersection { a, b, k } => {
            let sharp = eval_interval(a, bounds).max(eval_interval(b, bounds));
            Interval::new(
                sharp.lo,
                k.max(1e-6).mul_add(std::f32::consts::LN_2, sharp.hi),
            )
        }
        SdfNode::ExpSmoothSubtraction { a, b, k } => {
            let sharp = eval_interval(a, bounds).max(-eval_interval(b, bounds));
            Interval::new(
                sharp.lo,
                k.max(1e-6).mul_add(std::f32::consts::LN_2, sharp.hi),
            )
        }
        // 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 } => {
            // A zero / NaN / infinite factor (or one whose reciprocal overflows) makes the field NaN
            // or constant, so no finite enclosure exists: dividing by it built `[-1e-45, NaN]` or an
            // inverted interval. "Nothing is known" is the only sound answer.
            let inv = 1.0 / *factor;
            if !(factor.is_finite() && inv.is_finite()) {
                return Interval::entire();
            }
            let scaled = Vec3Interval {
                x: bounds.x * inv,
                y: bounds.y * inv,
                z: bounds.z * inv,
            };
            eval_interval(child, scaled) * *factor
        }
        SdfNode::ScaleNonUniform { child, factors } => {
            // same reasoning as `Scale`, per component
            if !(factors.x.is_finite()
                && factors.y.is_finite()
                && factors.z.is_finite()
                && (1.0 / factors.x).is_finite()
                && (1.0 / factors.y).is_finite()
                && (1.0 / factors.z).is_finite())
            {
                return Interval::entire();
            }
            let min_f = factors.x.min(factors.y).min(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);
            // `eval` multiplies by exactly `min_f`
            // ([`crate::transforms::scale::transform_scale_nonuniform`]), so the
            // enclosure is the child's interval times that one factor. Scaling
            // `hi` by `max_f` instead was looser and, on a wholly negative child
            // interval, **inverted** the bounds (`lo > hi`, e.g. `[−0.5, −0.25]`
            // with factors 0.8 / 2.5 → `[−0.4, −0.625]`): the debug assert in
            // `Interval::new` fired and release builds carried an interval whose
            // `hi` sat below the values it was supposed to enclose.
            let (lo, hi) = if min_f >= 0.0 {
                (d.lo * min_f, d.hi * min_f)
            } else {
                (d.hi * min_f, d.lo * min_f)
            };
            Interval::new(lo, hi)
        }
        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
            }
        }
        // LBS with fixed bones is one affine map q = A p + b (Σ w_i · bind_i ·
        // current_i / Σ w); push the box through it exactly (`affine_map`).
        SdfNode::SdfSkinning { child, bones } => match skinning_affine(bones) {
            Some((a, b)) => eval_interval(child, bounds.affine(&a, b)),
            None => 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).
        // The child interval is then pushed through `real::taper_bound`
        // (Jacobian ball + cone ∩ slab) in interval arithmetic.
        SdfNode::Taper {
            child,
            factor,
            reach,
        } => {
            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),
            );
            let d = eval_interval(
                child,
                Vec3Interval {
                    x: bounds.x * s,
                    y: bounds.y,
                    z: bounds.z * s,
                },
            );
            taper_bound_interval(d, bounds, *factor, *reach)
        }
        SdfNode::Displacement { child, strength } => {
            eval_interval(child, bounds).expand(strength.abs())
        }
        SdfNode::SineDisplacement {
            child, amplitude, ..
        } => eval_interval(child, bounds).expand(amplitude.abs()),
        // The fold sends `p` to the wedge `|atan2(z, x)| ≤ ha` at the same
        // radius, so `x = r·cos a`, `z = r·sin a` with `|a| ≤ ha`. Below
        // `ha = π/2` that is `x ∈ [0, max_r]`, `|z| ≤ max_r·sin ha`; from
        // `ha ≥ π/2` (count ≤ 2) the wedge reaches behind the axis and past
        // the ±z extremes, so both only bound by `max_r`. Using `sin ha`
        // unconditionally collapsed the count = 1 wedge to `z = 0`
        // (`sin π ≈ 0`): the enclosure then claimed `+0.05` on a cell where
        // the field is `−0.71`, and `alice-lol` proved "no overlap" there.
        SdfNode::PolarRepeat { child, count } => {
            let sector = std::f32::consts::TAU / (*count as f32);
            let ha = sector * 0.5;
            let max_r = bounds.length_xz().hi;
            let (x_lo, z_abs) = if ha >= std::f32::consts::FRAC_PI_2 {
                (-max_r, max_r)
            } else {
                (0.0, max_r * ha.sin())
            };
            eval_interval(
                child,
                Vec3Interval {
                    x: Interval::new(x_lo, max_r),
                    y: bounds.y,
                    z: Interval::new(-z_abs, z_abs),
                },
            )
        }
        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,
                },
            )
        }
        // Each iteration keeps the transform whose image is nearest the origin
        // (or the point itself): the reachable set after `iterations` steps is
        // inside the hull of the box and its images, iterated. The child is
        // evaluated over that hull and divided by the accumulated X-axis
        // scale, which lies between the smallest and largest of {1, |x_i|}^n.
        SdfNode::IFS {
            child,
            transforms,
            iterations,
        } => {
            let mut hull = bounds;
            let (mut s_lo, mut s_hi) = (1.0f32, 1.0f32);
            for t in transforms {
                let sx = t[2].mul_add(t[2], t[1].mul_add(t[1], t[0] * t[0])).sqrt();
                s_lo = s_lo.min(sx);
                s_hi = s_hi.max(sx);
            }
            for _ in 0..*iterations {
                let mut next = hull;
                for t in transforms {
                    let a = [t[0], t[1], t[2], t[4], t[5], t[6], t[8], t[9], t[10]];
                    let img = hull.affine(&a, Vec3::new(t[12], t[13], t[14]));
                    next = next.hull(&img);
                }
                hull = next;
            }
            let d = eval_interval(child, hull);
            let n = *iterations as i32;
            let scale = Interval::new(s_lo.powi(n).max(1e-6), s_hi.powi(n).max(1e-6));
            let recip = Interval::new(1.0 / scale.hi, 1.0 / scale.lo);
            d * recip
        }
        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(l.mul_add(-rho, d), l.mul_add(rho, d))
}

/// `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()
}

/// `glsl_mod(x, step)` over an interval: `x − step·floor(x / step)`.
///
/// The sawtooth is exact while the interval stays inside one period (a shift by
/// a multiple of `step`), and spans the whole `[0, step]` the moment it crosses
/// a period boundary. That crossing is a **jump**, which is why the stairs arms
/// can no longer assume a Lipschitz constant: the midpoint form
/// `d(centre) ± ρ` is invalid across it, and the corpus measured the resulting
/// inclusion violation (2.086·10⁻⁴ relative on `stairs_intersection`, child
/// interval escaping its parent, 2026-09-27).
#[inline]
fn mod_interval(x: Interval, step: f32) -> Interval {
    if step <= 0.0 || !step.is_finite() || !(x.lo.is_finite() && x.hi.is_finite()) {
        return Interval::new(0.0, step.max(0.0));
    }
    let k_lo = (x.lo / step).floor();
    let k_hi = (x.hi / step).floor();
    if k_lo == k_hi {
        // same period: a shift, so the width is preserved
        Interval::outward(x.lo - step * k_lo, x.hi - step * k_lo)
    } else {
        // crosses at least one jump: the range is the full tooth
        Interval::new(0.0, step)
    }
}

/// Stairs (stepped) minimum over intervals — the interval form of
/// [`crate::operations::stairs_min_r`], op by op.
///
/// Interval arithmetic loses the correlation between `px` and `py` (both are
/// rotations of the same `(a, b)`), so the result is wider than the point law's
/// range. It is *sound* though, which the midpoint + Lipschitz form was not.
fn ia_stairs_min(ia: Interval, ib: Interval, r: f32, n: f32) -> Interval {
    let n = n.max(1.0);
    let s = std::f32::consts::FRAC_1_SQRT_2;
    let sqrt2 = std::f32::consts::SQRT_2;
    let d = ia.min(ib);
    let py = (ia + ib) * s;
    let px = (ib - ia) * s;
    let rn = r / n;
    let off = (r - rn) * 0.5 * sqrt2;
    let px = px - (off - 0.5 * sqrt2 * rn);
    let py = py - off;
    let step = r * sqrt2 / n;
    let hs = step * 0.5;
    let px = mod_interval(px + Interval::point(hs), step) - hs;
    let d = d.min(py);
    let npx = (px + py) * s;
    let npy = (py - px) * s;
    let edge = 0.5 * rn;
    d.min((npx - edge).max(npy - edge))
}

/// 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
        }
        // The weighted norm's Lipschitz constant over the Euclidean unit
        // sphere is exact (alice-det-math derives it and checks it against a
        // brute-force sweep): `sqrt((w1+winf)^2 + 2*w1^2) + w2`. Weights that
        // are not a metric cannot reach here — `SdfNode::metric_ball` takes a
        // validated `MetricWeights` — but a hand-built node or a decoded
        // bytecode could, and then no bound is claimed.
        SdfNode::MetricBall {
            w_l1, w_l2, w_linf, ..
        } => alice_det_math::metric::MetricWeights::new(*w_l1, *w_l2, *w_linf)
            .map_or(f32::INFINITY, |w| w.lipschitz()),

        // Blending two unrelated fields across a skin adds
        // `|f_in - f_out| * |grad w|` to the gradient, and the first factor is
        // a property of the two fields over that skin, not of this node. No
        // finite bound follows from the law, so none is claimed; measure it
        // with `crate::measure::measure_tension`.
        SdfNode::MetricBlend { .. } => f32::INFINITY,

        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 `hash_noise_3d`
        // sampled on the xz plane (values in [-1, 1], smoothstep weights,
        // gain 0.48, lacunarity 2.1 with a rotation). Per-axis bound is
        // range(2)·slope(1.5) = 3, two active axes (x, z; y held at 0
        // contributes no gradient) combine as `3·√2`.
        // |∇fbm| ≤ Σ aᵢ·2.1ⁱ · 3·√2 = 6.415 → 6.42.
        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 } => t
            .abs()
            .mul_add(eval_lipschitz(b), (1.0 - t).abs() * eval_lipschitz(a)),
        // 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)
        // One affine map: L(child) · ‖A‖ (Frobenius ≥ spectral norm).
        SdfNode::SdfSkinning { child, bones } => match skinning_affine(bones) {
            Some((a, _)) => eval_lipschitz(child) * a.iter().map(|v| v * v).sum::<f32>().sqrt(),
            None => 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. The
        // value is a distance *bound* (`real::taper_bound`) safe for plain
        // sphere tracing, but it has no finite global Lipschitz constant.
        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),
        // The nearest-image choice switches between transforms: the folded
        // point jumps across the switching surfaces, so there is no finite
        // global constant (like domain repetition). `2 · L(child)` until 2.2.0
        // was unsound by 280× on a scale / rotate IFS.
        SdfNode::IFS { .. } => f32::INFINITY,
        // 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²) · 3 · √2`
/// (range-2 `hash_noise_3d`, not the old range-1 sin hash — see the
/// `SdfNode::Terrain` arm above).
const TERRAIN_FBM_GRAD: f32 = 6.42;
/// 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 {
    0.5f32.mul_add(v, (v * 0.25).mul_add(v, 1.0).sqrt())
}

/// 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
    }
}

/// Interval form of `real::taper_bound` (same steps, each in interval
/// arithmetic; `d` is the child interval over the tapered box).
fn taper_bound_interval(d: Interval, p: Vec3Interval, factor: f32, reach: [f32; 2]) -> Interval {
    let f_abs = factor.abs();
    if f_abs == 0.0 || d.lo == f32::NEG_INFINITY || d.hi == f32::INFINITY {
        return d;
    }
    let one = Interval::point(1.0);
    let eps = Interval::point(1e-6);
    // Reciprocal / quotient for positive intervals (all denominators here are
    // ≥ 1e-6 or ≥ 1).
    let recip = |a: Interval| Interval::new(1.0 / a.hi, 1.0 / a.lo);
    let div = |a: Interval, b: Interval| a * recip(b);

    let den_abs = (one - p.y * Interval::point(factor)).abs().max(eps);
    let s = recip(den_abs);
    let rho = (p.x.sqr() + p.z.sqr()).sqrt();
    let d_abs = d.abs();

    let j0 = s.max(one) + s.sqr() * rho * f_abs;
    let r = div(d_abs, j0).min(den_abs * (0.5 / f_abs));
    let den1 = (den_abs - r * f_abs).max(eps);
    let s1 = recip(den1);
    let j1 = s1.max(one) + s1.sqr() * (rho + r) * f_abs;
    let mag = div(d_abs, j1).min(r);
    let d_j = if d.hi < 0.0 {
        -mag
    } else if d.lo >= 0.0 {
        mag
    } else {
        (-mag).hull(mag)
    };

    if reach[0] >= 1e30 || reach[1] >= 1e30 {
        return d_j;
    }
    let k = reach[0] * f_abs;
    let inv_n = 1.0 / k.mul_add(k, 1.0).sqrt();
    let big_y = p.y - 1.0 / factor;
    let d_cone = (rho - big_y.abs() * k) * inv_n;
    let d_slab = p.y.abs() - reach[1];
    let d_region = d_cone.max(d_slab);

    if d.hi < 0.0 {
        d_j
    } else if d.lo >= 0.0 {
        d_j.max(d_region)
    } else {
        d_j.hull(d_j.max(d_region))
    }
}

/// The single affine map `q = A p + b` a fixed-pose skinning node applies
/// (`sdf_skinning`: Σ w_i · bind_i(current_i(p)) / Σ w over bones with
/// w ≥ 1e-6); `None` when no bone carries weight (identity).
fn skinning_affine(
    bones: &[crate::transforms::skinning::BoneTransform],
) -> Option<([f32; 9], Vec3)> {
    use glam::Mat4;
    let active: Vec<_> = bones.iter().filter(|b| b.weight >= 1e-6).collect();
    let total: f32 = active.iter().map(|b| b.weight).sum();
    if active.is_empty() || total <= 1e-6 {
        return None;
    }
    let mut m = Mat4::ZERO;
    for b in &active {
        let combined =
            Mat4::from_cols_array(&b.inv_bind_pose) * Mat4::from_cols_array(&b.current_pose);
        m += combined * (b.weight / total);
    }
    let c = m.to_cols_array();
    Some((
        [c[0], c[1], c[2], c[4], c[5], c[6], c[8], c[9], c[10]],
        Vec3::new(c[12], c[13], c[14]),
    ))
}

/// Reach of a taper child in its own (q) space, `[r_xz, r_y]`: `|q_xz| ≤ r_xz`
/// and `|q_y| ≤ r_y` over the child's AABB. `[INFINITY; 2]` when the child has
/// no finite AABB (a plane, an infinite repeat, …), which disables the cone
/// bound in `real::taper_bound`.
pub(crate) fn taper_reach(child: &SdfNode) -> [f32; 2] {
    use crate::compiled::{get_scene_aabb, CompiledSdfBvh};
    let Ok(bvh) = CompiledSdfBvh::try_compile(child) else {
        return [f32::INFINITY; 2];
    };
    let aabb = get_scene_aabb(&bvh);
    let (lo, hi) = (aabb.min(), aabb.max());
    if !aabb.is_valid() || !lo.is_finite() || !hi.is_finite() {
        return [f32::INFINITY; 2];
    }
    let reach = |a: f32, b: f32| a.abs().max(b.abs());
    let r_xz = reach(lo.x, hi.x).hypot(reach(lo.z, hi.z));
    let r_y = reach(lo.y, hi.y);
    if r_xz > UNBOUNDED_CHILD_SENTINEL || r_y > UNBOUNDED_CHILD_SENTINEL {
        [f32::INFINITY; 2]
    } else {
        [r_xz, r_y]
    }
}

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

    /// 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,
                    );
                }
            }
        }
    }

    /// Computed bounds must **contain** the exact value and sit within a few
    /// ulps of it — the contract of [`Interval::outward`] (3.1.1). Asserting
    /// exact equality would pin round-to-nearest, which is what made the
    /// enclosure non-rigorous in the first place.
    ///
    /// The slack is 8 ulps, not 1: the inputs are built with `Interval::new`,
    /// which already rounds outward, so a two-operation chain accumulates a
    /// handful (measured: `[1, 3] − [2, 5]` gives an upper bound of 1.0000005,
    /// i.e. 4 ulps above 1).
    #[track_caller]
    fn encloses(bound: f32, exact: f32, side: Ordering) {
        let ulps = 8.0 * f32::EPSILON * exact.abs().max(1.0);
        match side {
            Ordering::Less => assert!(
                bound <= exact && bound >= exact - ulps,
                "lower bound {bound} must be at or just below {exact}"
            ),
            _ => assert!(
                bound >= exact && bound <= exact + ulps,
                "upper bound {bound} must be at or just above {exact}"
            ),
        }
    }

    #[test]
    fn test_interval_ops() {
        let a = Interval::new(1.0, 3.0);
        let b = Interval::new(2.0, 5.0);
        encloses((a + b).lo, 3.0, Ordering::Less);
        encloses((a + b).hi, 8.0, Ordering::Greater);
        encloses((a - b).lo, -4.0, Ordering::Less);
        encloses((a - b).hi, 1.0, Ordering::Greater);
        // negation only swaps and flips signs, so it stays exact
        assert_eq!((-a).lo, -a.hi);
        assert_eq!((-a).hi, -a.lo);
    }

    #[test]
    fn test_interval_abs() {
        // `abs` selects existing bounds (or 0), so it is exact
        let a = Interval::new(1.0, 3.0);
        assert_eq!(a.abs().lo, a.lo);
        let b = Interval::new(-3.0, -1.0);
        assert_eq!(b.abs().hi, -b.lo);
        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();
        // 0 is an exact lower bound for a square; the upper bound is computed
        assert_eq!(s.lo, 0.0);
        encloses(s.hi, 9.0, Ordering::Greater);
    }

    #[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) - 0.6f32.mul_add(PERLIN_GRAD, 1.0)).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.5f32.mul_add(v, (0.25 * v).mul_add(v, 1.0).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);
    }
}