Merge hands/cultmath-interval-ops: CultMath intervals r1-r4, the tile-amortized march proved on fields with empty space

iv_frustum_ball widened by its derived float32 bound, envelope-then-noise
fields, the analytic constant-density body, one shipped void (Rh 198, S 50),
LOD agreement against a converged march, and the wells contract as ruled
(cultmath-tapes:ruling:wells-overhead-plus-floor). Verdict
cut-interval-ops.s3 holds; one Low finding (LodStep footprint branch unpinned)
is recorded for the next interval cut.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

packages/cultmath/README.md

@@ -33,8 +33,20 @@ shader semantics.
 - `shaders/CultMath.hlsl`, a canonical HLSL mirror include for shader-side
   parity. HLSL already owns `float2`, `float3`, and `float4`; the include
   exposes `cultmath_*` functions for shared semantics HLSL does not provide.
-  It includes `shaders/CultMath.Phacelle.hlsl` (MPL-2.0) from the same
-  directory, so ship both files.
+  It includes `shaders/CultMath.Phacelle.hlsl` (MPL-2.0) and
+  `shaders/CultMath.Interval.hlsl` from the same directory, so ship all three.
+- interval arithmetic for culling empty space with a proof: `iv_*` functions
+  over `float2(lo, hi)` (`iv_add`, `iv_mul`, `iv_exp`, `iv_smoothstep`, ...)
+  and `iv_snoise_ball`/`iv_fbm_ball`, which bound `snoise` and `fbm_grad` over a
+  ball, and `iv_frustum_ball`, the ball around a screen tile's rays over a depth
+  segment, so one probe serves the whole tile (its radius carries its own
+  float32 rounding bound). Every one returns an interval
+  containing the function's value at every point of its input, so a raymarch
+  can skip a segment whose `hi` is below its cutoff; `cultmath_iv_*` in HLSL.
+  The bounds compose with a consumer's own analytic envelope (a height fog, a
+  carved sphere): its interval over the segment proves empty space in one probe,
+  and the noise ball is consulted only near the surface. See `docs/design.md`,
+  "Intervals".
 - `Voronoi.SampleTones`, a C# batch surface that calls the Rust
   `cultmath-core` native kernel when `cultmath_core` is available and falls back
   to the managed parity path otherwise.

packages/cultmath/docs/design.md

@@ -341,6 +341,230 @@ Known C# friction: under `using static CultMath.math;` the constructor
 functions (`float3(...)`, `float3x3(...)`) hide the type names in member access,
 so static members need a qualified type (`CultMath.float3x3.Euler(...)`).
 
+## Intervals
+
+`math.Interval.cs` and `shaders/CultMath.Interval.hlsl` (included from
+`CultMath.hlsl` after the Phacelle include) carry interval arithmetic for
+culling empty space with a proof. Their invariant, **intervals-enclose**: every
+interval function returns a `[lo, hi]` that contains `f(x)` for every `x` in
+its input region, where `f(x)` is the same function evaluated pointwise in
+float32. A consumer that skips a region because `hi < cutoff` has skipped
+nothing that was there.
+
+Representation: an interval is `float2(lo, hi)` with `lo <= hi`, both finite.
+There is no empty interval; a caller that needs "provably empty" tests
+`hi < cutoff`. `iv_point(x)` is `float2(x, x)`.
+
+Ops: `iv_point`, `iv_add`, `iv_sub`, `iv_neg`, `iv_mul` (the least and greatest
+of the four corner products), `iv_scale(iv, float)`, `iv_abs` and `iv_sqr`
+(an interval straddling zero starts at 0), `iv_min`, `iv_max`, `iv_sqrt`
+(`lo` clamped at 0), `iv_exp`, `iv_clamp(iv, float, float)`, `iv_saturate`,
+`iv_smoothstep(float, float, iv)` (its ordered endpoints) and
+`iv_lerp(iv, iv, float)`. No others until a consumer names one. Each has an
+HLSL mirror named `cultmath_iv_*`, compared bit for bit by
+`HlslSourceCompatibilityTests`.
+
+The ulp-widening rule. Each bound is computed with the pointwise function's
+own float32 operations. IEEE round-to-nearest is monotone, so an op that is
+monotone on its region encloses its float32 evaluation exactly, with no
+widening; that covers every op above except two, which widen by exactly one
+ulp, in the implementation and never in the test:
+
+- `iv_exp`: `exp` is not required to be correctly rounded, so a libm need not
+  be monotone. .NET's `MathF.Exp` on Linux steps down 0 times over every float
+  in [-87, 88.7] (`IntervalTests.MeasureMonotonicity`), so on that CPU the
+  widening is a margin; it is there for other libms and GPUs.
+- `iv_smoothstep`: `t * t * (3 - 2t)` in float32 steps down 337,095 times over
+  the floats in [0, 1], each time by at most 1 ulp (`MeasureMonotonicity`).
+  Since `t = saturate((x - minimum) / (maximum - minimum))` is itself monotone
+  in `x`, that covers every `smoothstep`.
+
+`iv_lerp` is the natural interval extension of `x + (y - x) * amount`, which is
+how `lerp` evaluates; every step is monotone, so it encloses with no widening
+for any `amount`, and it is tight only when `a` is a point, because `x` appears
+twice.
+
+`IntervalTests` runs every op over 10,000 seeded intervals (centres in
+[-20, 20], widths log-uniform in [1e-4, 10], every eighth a point) and 16
+points inside each, endpoints and zero included, with no tolerance; and for
+every op it raises `lo` and lowers `hi` by one ulp and requires the harness to
+catch each. A new op is one entry in that table.
+
+Noise bounds. `iv_snoise_ball(c, r)` is `[n - L r, n + L r]` intersected with
+`[-1, 1]`, where `n = snoise(c)` and `L = SNOISE_LIPSCHITZ`: one `snoise`, value
+only. No interval function reads `snoise_grad`, and none has a value-only or
+gradient twin. `iv_fbm_ball(c, r, octaves, lacunarity, gain)` is the octave sum
+of `iv_snoise_ball(c * f_i, r * |f_i|)` scaled by `a_i`, compounding frequency
+and amplitude exactly as `fbm_grad` does, octaves clamped to [0, 16]. A domain
+warp that moves a point by at most `D` is enclosed by the radius `r + D`
+(`NoiseBoundTests.WarpedPointsStayEnclosed`); this is what the gamecult.org
+ground's march rests on.
+
+`iv_frustum_ball(m_c, z0, z1, footprintPerDepth, warp)` is the ball one probe
+serves a whole screen tile with. In a pinhole camera frame looking down `+z`,
+where the ray of slope `m` is the points `(m z, z)`, it encloses every point of
+every ray whose slope lies within `footprintPerDepth` of `m_c` (`N / (sqrt(2) f)`
+for an `N x N` tile at focal length `f` pixels: the full pixel footprint, so
+sub-pixel jitter is covered), over depths `[z0, z1]`, each point moved by a warp
+of length at most `warp`. The centre is `(m_c z_m, z_m)`, `z_m` the mid-depth;
+the radius is `((z1 - z0) / 2) |(m_c, 1)| + z1 footprintPerDepth + warp`. A point
+of slope `m_c + e` at depth `z` is the centre plus `(z - z_m)(m_c, 1)` plus
+`z (e, 0)`. The returned ball is float32, and on a degenerate (`z0 == z1`),
+thin or far segment the rounding of the centre alone exceeds the gap the
+triangle inequality leaves, so the radius is widened by
+`(|c|_1 + radius) 2^-20`. With `u = 2^-24`, each centre component carries at most
+two roundings (`z_m`, then the product), so the float centre is within
+`2u |c|_1` of the exact one; the radius carries at most seven (`z1 - z0`, the
+axis length at `3u`, the product, `z1 footprintPerDepth`, the sum, `+ warp`), so
+the exact radius is at most `(1 + 7u)` times the float one; the widening rounds
+once more. `2^-20 = 16u` covers `8u (|c|_1 + radius)` twice over, and fma
+contraction only removes roundings, so HLSL and GLSL compilers are covered by
+the same bound. The caller rotates and translates the centre into world space and scales
+centre and radius by its noise frequency; the rounding of that transform is the
+caller's. `NoiseBoundTests.TileBallEnclosesEveryRaySegment` checks 2,000 seeded
+tiles x 64 points of the r2 domain and 2,000 tiles x 8 corner points of each of
+seven extreme families (`z0 == z1`, `z1 - z0 = 1e-6 z1`, `|m_c|` to 1000, `z` to
+`1e7`, `z0 = 0` with `z1` to `1e-3`, `f` down to 1, warp to `1e6`), corner flows
+full-length and outward, distances in double, with no tolerance. It pins the ball
+to the derivation both ways (`r <= ball.w <= r + 2^-19 (r + |c|_1)`), so a ball
+that forgets the widening fails and one grown past it fails too. A moving camera needs nothing
+more: each frame's probes use that frame's camera, and no ball, mask or probe
+result is carried from one frame to the next.
+
+The noise ball alone proves little: at `L` near 10 it bounds `snoise` to a width
+of `2 L r`, which covers all of `[-1, 1]` once `L r` reaches 1. Empty space
+becomes provable through an analytic envelope. The consumer's density is a
+pre-distortion SDF (a height fog, a well map, a carved sphere) that noise
+displaces only near its surface, and the consumer composes the bound from `iv_*`
+ops. The envelope is taken over the slice's exact box (`iv_mul` of the depth
+range by each slope range, then the SDF's own monotone pieces), together with its
+fade. The noise is `iv_snoise_ball` over `iv_frustum_ball`, scaled by the noise
+frequency. The warp moves only the noise argument, so it enlarges the ball and
+the envelope's box takes none of it. Where the fade's upper end is 0 the noise
+term is identically 0 and no `snoise` is evaluated. Where the envelope already
+proves the slice empty with the noise at its full range, the ball is not
+evaluated either. A probe is worth attempting only when the best result any
+centre value could give would prove the slice, and a ball with `L r >= 2` is
+`[-1, 1]` whatever `snoise(c)` is, so it costs nothing. CultMath owns none of the
+envelope: the interval files have no envelope function, and each consumer composes
+its own. `NoiseBoundTests.EnvelopeBoundEnclosesDensity` checks that composition,
+the one the site and Aetheria use, over 2,000 seeded slices x 64 warped points
+with no tolerance.
+
+`L` has provenance, not a proof. `NoiseBoundTests.MeasureLipschitz` (slow,
+explicit) takes the largest `|snoise_grad|` over 1e6 seeded points in
+[-256, 256]^3 (7.1933103), refines the 1e4 largest by gradient ascent
+(9.181147), and multiplies by 1.10: `L = 10.099261`. It also reaches
+`|snoise| = 0.9718335` by the same ascent on the value, which is what lets the
+bound intersect with `[-1, 1]`. `LipschitzConstantPinsSampledGradients`
+checks 1e5 more seeded points against `L` and climbs again from the start
+`MeasureLipschitz` prints as its witness, failing if `L` is not 1.10 times that
+maximum to within half a percent, so a 1% move of `L` fails. The enclosure tests
+over 2,000 balls x 64 points are the defence. A change to the `snoise` kernel
+must re-run `MeasureLipschitz` and re-pin `L` and the witness.
+
+What the ball bounds guarantee: enclosure of `snoise` and `fbm_grad(...).w` over
+the ball, in float32. What they do not: anything about tapes, pruning, choice
+tracking or affine forms; those are later steps of
+`docs/cultmath-tape-target.md` (CultLib root). How loose they are is a number:
+`NoiseBoundTests.TightnessReport` prints the mean interval width over the
+sampled range: 5.08 for `iv_snoise_ball` at `r = 0.05`, 1.96 at 0.25 and 1.28
+at 1.0; 6.64 and 2.82 for 4-octave `iv_fbm_ball` at 0.05 and 0.25.
+
+The saving is measured where there is empty space to find.
+`NoiseBoundTests.IntervalSkipHalvesEvaluations` draws screen tiles, each with
+its own camera, and runs one pre-pass per tile over the dense grid's cell ranges:
+it starts with the whole grid, skips and doubles a range it proves empty, halves
+one it cannot, and runs to the grid's end. It then draws the tile's `N^2` ray
+slopes, one per pixel, jittered inside it. Each ray is marched three times with
+the same integrator and the same early-out at transmittance 0.02: over every
+cell (dense), over the tile's unproven cells (the tile march), and over the
+oracle mask, the cells whose mid-depth density is nonzero for at least one of
+the tile's rays. The tile march and the oracle march integrate the dense march's
+nonzero cells, so their transmittance agrees with it exactly in every
+configuration. The oracle is a ceiling the pre-pass cannot move: dense cost over
+oracle cost. Against it the test prints the efficiency (oracle cost over the
+tile march's cost with its probes) and the cull fraction (the oracle-empty cells
+the pre-pass proves empty). The probe overhead is the tile march's cost over that
+cost plus the probes'. Cost is counted in `snoise` evaluations, with an envelope
+evaluation or probe weighted 0.2. That weight is an estimate, and both counts are
+printed. The flow is not counted. At `N = 8`, in Aetheria's units on its
+256-cell quadratic grid out to the Main Camera's far plane, 2048:
+
+- (a) Height fog, camera above the safety band, rays level and up: the whole
+  grid is one envelope probe per tile, against 256 envelope evaluations per ray,
+  16384x. Every cell is oracle-empty and the pre-pass proves them all.
+- (b) Inside the fog: 1.15x, and 0.99x with the warp. Most cells are fog, and
+  the pre-pass probes each one.
+- (c) The zone bowl and four wells, the camera gazing across. With no warp:
+  2.74x, oracle ceiling 7.26x, efficiency 0.378, cull fraction 0.952, probe
+  overhead 0.961. With the warp `D = 60`: 1.39x, oracle ceiling 7.17x,
+  efficiency 0.194, cull fraction 0.842, probe overhead 0.993. With the warp the
+  pre-pass cuts envelope evaluations from 177.46 to 33.49 per ray, but the
+  oracle march costs 0.194 of the tile march, so about four fifths of the tile
+  march's cost goes to cells no ray samples nonzero: the warp bound swamps the
+  noise ball, so the ball cannot prove those cells empty. That gap belongs to
+  the tape target's affine forms, not to a better pre-pass.
+- (d) r1's uniform slab, which has no envelope: 1.23x, and 1.48x with no warp.
+- (e) The void at the shipped point (ruling `operator-shipped-void`): `Rh = 198`,
+  `S = 50`, `Rc = Rh + S = 248`, `e = -ln 1.5 / ln(1 - (Rh / Rc)^2)`, the
+  shipped camera, 200 tiles. The fixed-step reference takes 36.34 steps per
+  pixel. The same cells masked by the pre-pass take 20.13, with identical
+  transmittance. The footprint-aware march takes 13.54: 2.68x fewer. `snoise`
+  per pixel goes from 32.98 to 26.69, and combined cost from 40.25 to 29.52 with
+  the probes, 1.36x. Every one of the 200 tiles proves a body. Before the body
+  both marches are measured against a converged one (midpoint quadrature at
+  `h = 0.5` over every cell, trusting no pre-pass): the relative optical-depth
+  error at the body start is at most 0.0608 for the fixed-step march and 0.08639
+  for the footprint march, 1.42x. Over the grid of hollow radius, camera offset
+  and ramp width the steps ratio runs from 1.57x (`Rh = 100`, low, `S = 150`) to
+  27.45x (`Rh = 1000`, centred, `S = 10`), the cost ratio from 1.26x to 3.89x,
+  and the depth-error ratio from 0.85 to 4.72 (`Rh = 200`, centred, `S = 150`,
+  where the errors are 0.01504 and 0.003188). LOD-far takes 14.19x fewer steps;
+  its depth errors, 0.2849 and 0.205, are not comparable, because its
+  fixed-step reference keeps the octaves the footprint march truncates.
+
+The analytic body (ruling `operator-analytic-wall-body`). When a single dense
+cell's own bound is a positive point, the pre-pass spends one more probe over
+the rest of the grid. If that bound is a positive point too, every point of
+every ray of the tile from there on has exactly that density, so the footprint
+march takes one closed-form step, `exp(-BodyDepth)`, from the body start to the
+grid's end and stops. `BodyDepth` is density x extinction x length x `|dir|`.
+`NoiseBoundTests.AnalyticBodyMatchesFineMarch` pins it: the bound, the mask and
+the pointwise densities agree exactly at the body, and `BodyDepth` equals the
+summed fixed-step depth of the body cells within 1e-4 relative. The comparison
+is in the depth domain, so LOD-far's deep bodies do not underflow.
+
+The footprint-truncated void weights its octaves by `w_i(z)`, decreasing in `z`.
+Its bound takes each weight over `[w_i(z1), 1]` rather than
+`[w_i(z1), w_i(z0)]`, because one bound serves both the full field (`w = 1`) and
+the truncated one. `[w_i(z1), 1]` contains both, and `iv_mul` is monotone in its
+interval argument, so the composed interval encloses both fields.
+`NoiseBoundTests.WeightBoundEnclosesOneSignedNoise` checks it on LOD-far over
+2,000 segments: each weight lies in its bound, and `w(z) n` lies in
+`iv_mul(bound, n)` for a one-signed `n`. A lower end of `w(z0)` fails it.
+
+The contracts, at `N = 8`, the lower of warp 0 and `D`:
+
+- (a) at least 2x combined cost;
+- (c) (ruling `wells-overhead-plus-floor`): probe overhead, the tile march over
+  the tile march plus its probes, at least 0.90; and, hard, the pre-pass proves
+  at least half of the oracle-empty cells empty, so a pre-pass that skips
+  nothing fails. The oracle ceiling, the efficiency against it and the cull
+  fraction are printed, not asserted;
+- (e) at the shipped void: the masked fixed-step march agrees with the reference
+  exactly; the footprint march's optical depth at the body start is off the
+  converged march's by at most twice the fixed-step march's error; and it takes
+  at most half the reference's steps per pixel.
+
+This run meets every contract: (a) 16384x, (c) probe overhead 0.961 and cull
+fraction 0.842, (e) 2.68x fewer steps and a depth error 1.42x the fixed-step
+march's. The fields, the grids and the 0.2 weight are not tuned toward the
+contracts.
+
+Consumers: the gamecult.org ground shader (through the GLSL lowering) and
+Aetheria's nebula raymarch at its CultMath pin bump.
+
 ## Rules
 
 - Keep type names HLSL-shaped: `float2`, `float3`, `float4`, `float3x3`, and `math`.

packages/cultmath/scripts/build-unity-package.ps1

@@ -82,7 +82,9 @@ $requiredTemplateFiles = @(
   "Shaders\CultMath.hlsl",
   "Shaders\CultMath.hlsl.meta",
   "Shaders\CultMath.Phacelle.hlsl",
-  "Shaders\CultMath.Phacelle.hlsl.meta"
+  "Shaders\CultMath.Phacelle.hlsl.meta",
+  "Shaders\CultMath.Interval.hlsl",
+  "Shaders\CultMath.Interval.hlsl.meta"
 )
 foreach ($relativePath in $requiredTemplateFiles) {
   if (-not (Test-Path -LiteralPath (Join-Path $templateRoot $relativePath) -PathType Leaf)) {

packages/cultmath/shaders/CultMath.Interval.hlsl

@@ -0,0 +1,122 @@
+#ifndef CULTMATH_INTERVAL_HLSL
+#define CULTMATH_INTERVAL_HLSL
+
+// Interval arithmetic over float2(lo, hi), lo <= hi, both finite. Every cultmath_iv_* function returns
+// an interval containing f(x) for every x in its input intervals, mirroring math.Interval.cs bit for
+// bit; see that file and docs/design.md ("Intervals") for each op's enclosure and its ulp widening.
+// Included by CultMath.hlsl after its own definitions of cultmath_smoothstep and cultmath_snoise.
+
+// Lipschitz constant of cultmath_snoise(float3); provenance in math.Interval.cs (SNOISE_LIPSCHITZ).
+static const float CULTMATH_SNOISE_LIPSCHITZ = 10.099261;
+
+float2 cultmath_iv_point(float value) { return float2(value, value); }
+float2 cultmath_iv_add(float2 a, float2 b) { return float2(a.x + b.x, a.y + b.y); }
+float2 cultmath_iv_sub(float2 a, float2 b) { return float2(a.x - b.y, a.y - b.x); }
+float2 cultmath_iv_neg(float2 a) { return float2(-a.y, -a.x); }
+
+float2 cultmath_iv_mul(float2 a, float2 b)
+{
+    float p0 = a.x * b.x;
+    float p1 = a.x * b.y;
+    float p2 = a.y * b.x;
+    float p3 = a.y * b.y;
+    return float2(min(min(p0, p1), min(p2, p3)), max(max(p0, p1), max(p2, p3)));
+}
+
+float2 cultmath_iv_scale(float2 a, float s) { return s >= 0.0 ? float2(a.x * s, a.y * s) : float2(a.y * s, a.x * s); }
+
+float2 cultmath_iv_abs(float2 a)
+{
+    if (a.x >= 0.0)
+        return a;
+    if (a.y <= 0.0)
+        return float2(-a.y, -a.x);
+    return float2(0.0, max(-a.x, a.y));
+}
+
+float2 cultmath_iv_min(float2 a, float2 b) { return float2(min(a.x, b.x), min(a.y, b.y)); }
+float2 cultmath_iv_max(float2 a, float2 b) { return float2(max(a.x, b.x), max(a.y, b.y)); }
+
+float2 cultmath_iv_sqr(float2 a)
+{
+    float l = a.x * a.x;
+    float h = a.y * a.y;
+    if (a.x >= 0.0)
+        return float2(l, h);
+    if (a.y <= 0.0)
+        return float2(h, l);
+    return float2(0.0, max(l, h));
+}
+
+float2 cultmath_iv_sqrt(float2 a) { return float2(sqrt(max(a.x, 0.0)), sqrt(max(a.y, 0.0))); }
+
+// Widened by one ulp on each side: exp is not required to be correctly rounded.
+float2 cultmath_iv_exp(float2 a)
+{
+    float lo = exp(a.x);
+    float hi = exp(a.y);
+    return float2(lo > 0.0 ? asfloat(asuint(lo) - 1u) : 0.0, asfloat(asuint(hi) + 1u));
+}
+
+float2 cultmath_iv_clamp(float2 a, float minimum, float maximum)
+{
+    return float2(cultmath_clamp(a.x, minimum, maximum), cultmath_clamp(a.y, minimum, maximum));
+}
+
+float2 cultmath_iv_saturate(float2 a) { return float2(cultmath_saturate(a.x), cultmath_saturate(a.y)); }
+
+// The endpoints, ordered, widened by one ulp inside [0, 1]: the float32 evaluation of smoothstep is
+// monotone only to within one ulp.
+float2 cultmath_iv_smoothstep(float minimum, float maximum, float2 a)
+{
+    float s0 = cultmath_smoothstep(minimum, maximum, a.x);
+    float s1 = cultmath_smoothstep(minimum, maximum, a.y);
+    float lo = min(s0, s1);
+    float hi = max(s0, s1);
+    return float2(lo > 0.0 ? asfloat(asuint(lo) - 1u) : 0.0, hi < 1.0 ? asfloat(asuint(hi) + 1u) : 1.0);
+}
+
+// The natural interval extension of x + (y - x) * amount, which is how cultmath_lerp evaluates.
+float2 cultmath_iv_lerp(float2 a, float2 b, float amount) { return cultmath_iv_add(a, cultmath_iv_scale(cultmath_iv_sub(b, a), amount)); }
+
+// [n - L r, n + L r] intersected with [-1, 1], n = cultmath_snoise(centre). One snoise, value only.
+float2 cultmath_iv_snoise_ball(float3 centre, float radius)
+{
+    float n = cultmath_snoise(centre);
+    float e = CULTMATH_SNOISE_LIPSCHITZ * radius;
+    return float2(max(n - e, -1.0), min(n + e, 1.0));
+}
+
+// Octave sum of cultmath_iv_snoise_ball, compounding frequency and amplitude as the value-and-gradient fBm does.
+float2 cultmath_iv_fbm_ball(float3 centre, float radius, int octaves, float lacunarity, float gain)
+{
+    octaves = clamp(octaves, 0, 16);
+    float amplitude = 1.0;
+    float frequency = 1.0;
+    float2 sum = float2(0.0, 0.0);
+    for (int i = 0; i < octaves; i++)
+    {
+        sum = cultmath_iv_add(sum, cultmath_iv_scale(cultmath_iv_snoise_ball(centre * frequency, radius * abs(frequency)), amplitude));
+        frequency *= lacunarity;
+        amplitude *= gain;
+    }
+
+    return sum;
+}
+
+// A ball (xyz centre, w radius) enclosing every point of every ray of a screen tile over the depth segment
+// [z0, z1], warped by at most warp. Camera frame: pinhole at the origin looking down +z, ray (m z, z).
+// footprintPerDepth = N / (sqrt(2) f) for an N x N tile at focal length f. The radius is widened by
+// (|c|_1 + radius) 2^-20, which covers the ball's own float32 rounding. Derivation in math.Interval.cs.
+float4 cultmath_iv_frustum_ball(float2 centreSlope, float z0, float z1, float footprintPerDepth, float warp)
+{
+    float zm = (z0 + z1) * 0.5;
+    float axis = sqrt(centreSlope.x * centreSlope.x + centreSlope.y * centreSlope.y + 1.0);
+    float cx = centreSlope.x * zm;
+    float cy = centreSlope.y * zm;
+    float radius = (z1 - z0) * 0.5 * axis + z1 * footprintPerDepth + warp;
+    radius += (abs(cx) + abs(cy) + abs(zm) + radius) * 9.5367431640625e-7;
+    return float4(cx, cy, zm, radius);
+}
+
+#endif

packages/cultmath/shaders/CultMath.Interval.hlsl.meta

@@ -0,0 +1,7 @@
+fileFormatVersion: 2
+guid: 4d8a899561654ac990cab0bf48658cc1
+ShaderIncludeImporter:
+  externalObjects: {}
+  userData: 
+  assetBundleName: 
+  assetBundleVariant: 

packages/cultmath/shaders/CultMath.hlsl

@@ -504,6 +504,7 @@ CultCellular cultmath_cellular(float3 p)
 }
 
 #include "CultMath.Phacelle.hlsl"
+#include "CultMath.Interval.hlsl"
 
 float cultmath_value_noise(float2 position)
 {

packages/cultmath/src/CultMath/CultMath.csproj

@@ -24,6 +24,7 @@
     <None Include="..\..\README.md" Pack="true" PackagePath="\" />
     <None Include="..\..\shaders\CultMath.hlsl" Pack="true" PackagePath="contentFiles\any\any\shaders\CultMath.hlsl" CopyToOutputDirectory="PreserveNewest" />
     <None Include="..\..\shaders\CultMath.Phacelle.hlsl" Pack="true" PackagePath="contentFiles\any\any\shaders\CultMath.Phacelle.hlsl" CopyToOutputDirectory="PreserveNewest" />
+    <None Include="..\..\shaders\CultMath.Interval.hlsl" Pack="true" PackagePath="contentFiles\any\any\shaders\CultMath.Interval.hlsl" CopyToOutputDirectory="PreserveNewest" />
     <None Include="..\..\THIRD-PARTY-NOTICES.md" Pack="true" PackagePath="\" />
   </ItemGroup>
 

packages/cultmath/src/CultMath/math.Interval.cs

@@ -0,0 +1,199 @@
+namespace CultMath;
+
+public static partial class math
+{
+    // Interval arithmetic over float2(lo, hi), lo <= hi, both finite (design.md, "Intervals"). Every
+    // iv_* function returns an interval that contains f(x) for every x in its input intervals, where
+    // f(x) is the same function evaluated pointwise in float32, mirrored by cultmath_iv_* in
+    // shaders/CultMath.Interval.hlsl. Each bound is computed with the pointwise function's own float32
+    // operations; because IEEE round-to-nearest is monotone, an op that is monotone on its region
+    // encloses its float32 evaluation with no widening. Where float32 evaluation is not monotone, or
+    // the libm is not required to be, the bound is widened by exactly one ulp, and the doc comment
+    // says so. There is no empty interval: a caller that needs "provably empty" tests hi < cutoff.
+
+    // Lipschitz constant of snoise(float3) in its input: an upper bound on the length of its analytic
+    // gradient. Provenance: NoiseBoundTests.MeasureLipschitz (slow, explicit) takes the largest gradient
+    // length over 1e6 seeded points, refines the 1e4 largest by gradient ascent, and multiplies the
+    // refined maximum by 1.10; NoiseBoundTests.LipschitzConstantPinsSampledGradients pins this value
+    // against a fast re-measure and the ascent from MeasureLipschitz's witness.
+    // Empirical with a margin, not a proof; the enclosure tests are the defence. A change to the
+    // snoise kernel must re-run MeasureLipschitz and re-pin this constant.
+    public const float SNOISE_LIPSCHITZ = 10.099261f;
+
+    /// <summary>The point interval [x, x].</summary>
+    public static float2 iv_point(float value) => new(value, value);
+
+    /// <summary>Encloses x + y for x in a, y in b. Monotone: no widening.</summary>
+    public static float2 iv_add(float2 a, float2 b) => new(a.x + b.x, a.y + b.y);
+
+    /// <summary>Encloses x - y for x in a, y in b. Monotone: no widening.</summary>
+    public static float2 iv_sub(float2 a, float2 b) => new(a.x - b.y, a.y - b.x);
+
+    /// <summary>Encloses -x for x in a. Exact.</summary>
+    public static float2 iv_neg(float2 a) => new(-a.y, -a.x);
+
+    /// <summary>
+    /// Encloses x * y for x in a, y in b: the least and greatest of the four corner products. The
+    /// product is monotone in each factor for a fixed sign of the other, so no widening.
+    /// </summary>
+    public static float2 iv_mul(float2 a, float2 b)
+    {
+        var p0 = a.x * b.x;
+        var p1 = a.x * b.y;
+        var p2 = a.y * b.x;
+        var p3 = a.y * b.y;
+        return new(min(min(p0, p1), min(p2, p3)), max(max(p0, p1), max(p2, p3)));
+    }
+
+    /// <summary>Encloses x * s for x in a. Monotone: no widening.</summary>
+    public static float2 iv_scale(float2 a, float s) => s >= 0.0f ? new(a.x * s, a.y * s) : new(a.y * s, a.x * s);
+
+    /// <summary>Encloses abs(x) for x in a; an interval straddling zero gives [0, max(-lo, hi)]. Exact.</summary>
+    public static float2 iv_abs(float2 a)
+    {
+        if (a.x >= 0.0f)
+            return a;
+        if (a.y <= 0.0f)
+            return new(-a.y, -a.x);
+        return new(0.0f, max(-a.x, a.y));
+    }
+
+    /// <summary>Encloses min(x, y) for x in a, y in b. Monotone: no widening.</summary>
+    public static float2 iv_min(float2 a, float2 b) => new(min(a.x, b.x), min(a.y, b.y));
+
+    /// <summary>Encloses max(x, y) for x in a, y in b. Monotone: no widening.</summary>
+    public static float2 iv_max(float2 a, float2 b) => new(max(a.x, b.x), max(a.y, b.y));
+
+    /// <summary>
+    /// Encloses x * x for x in a; an interval straddling zero gives [0, max(lo * lo, hi * hi)]. Tighter
+    /// than iv_mul(a, a), which cannot know both factors are the same x. Monotone on each side of zero:
+    /// no widening.
+    /// </summary>
+    public static float2 iv_sqr(float2 a)
+    {
+        var l = a.x * a.x;
+        var h = a.y * a.y;
+        if (a.x >= 0.0f)
+            return new(l, h);
+        if (a.y <= 0.0f)
+            return new(h, l);
+        return new(0.0f, max(l, h));
+    }
+
+    /// <summary>
+    /// Encloses sqrt(x) for x in a with x >= 0 (sqrt's domain): lo is clamped at 0, so [-4, 9] gives [0, 3]
+    /// and an interval wholly below zero gives [0, 0]. sqrt is correctly rounded, hence monotone: no
+    /// widening.
+    /// </summary>
+    public static float2 iv_sqrt(float2 a) => new(sqrt(max(a.x, 0.0f)), sqrt(max(a.y, 0.0f)));
+
+    /// <summary>
+    /// Encloses exp(x) for x in a. Widened by one ulp on each side, because exp is not required to be
+    /// correctly rounded, so a libm need not be monotone. Where exp(hi) overflows to infinity, hi is NaN,
+    /// which a cull test hi &lt; cutoff reads as "not provably below", the safe side.
+    /// </summary>
+    public static float2 iv_exp(float2 a)
+    {
+        var lo = exp(a.x);
+        var hi = exp(a.y);
+        return new(lo > 0.0f ? asfloat(asuint(lo) - 1u) : 0.0f, asfloat(asuint(hi) + 1u));
+    }
+
+    /// <summary>Encloses clamp(x, minimum, maximum) for x in a. Monotone: no widening.</summary>
+    public static float2 iv_clamp(float2 a, float minimum, float maximum) =>
+        new(clamp(a.x, minimum, maximum), clamp(a.y, minimum, maximum));
+
+    /// <summary>Encloses saturate(x) for x in a. Monotone: no widening.</summary>
+    public static float2 iv_saturate(float2 a) => new(saturate(a.x), saturate(a.y));
+
+    /// <summary>
+    /// Encloses smoothstep(minimum, maximum, x) for x in a. smoothstep is monotone in x (increasing when
+    /// minimum &lt; maximum, decreasing when minimum &gt; maximum), so the bound is its value at the two
+    /// endpoints, ordered. Its float32 evaluation t * t * (3 - 2t) is not monotone, by at most one ulp
+    /// (IntervalTests.MeasureMonotonicity), so the bound is widened by one ulp inside [0, 1].
+    /// </summary>
+    public static float2 iv_smoothstep(float minimum, float maximum, float2 a)
+    {
+        var s0 = smoothstep(minimum, maximum, a.x);
+        var s1 = smoothstep(minimum, maximum, a.y);
+        var lo = min(s0, s1);
+        var hi = max(s0, s1);
+        return new(lo > 0.0f ? asfloat(asuint(lo) - 1u) : 0.0f, hi < 1.0f ? asfloat(asuint(hi) + 1u) : 1.0f);
+    }
+
+    /// <summary>
+    /// Encloses lerp(x, y, amount) = x + (y - x) * amount for x in a, y in b, by evaluating that expression
+    /// with iv_sub, iv_scale and iv_add. Each step is monotone, so the result encloses the float32 lerp
+    /// with no widening, for any amount. It is the natural interval extension: x appears twice, so for
+    /// non-point a the width is (width(a) + width(b)) * |amount| + width(a), wider than the true range.
+    /// </summary>
+    public static float2 iv_lerp(float2 a, float2 b, float amount) => iv_add(a, iv_scale(iv_sub(b, a), amount));
+
+    /// <summary>
+    /// Encloses snoise(x) for every x within distance radius of centre: [n - L r, n + L r] intersected
+    /// with [-1, 1], where n = snoise(centre) and L = SNOISE_LIPSCHITZ. radius must be at least 0. One
+    /// snoise evaluation, value only. A domain warp that moves x by at most D is enclosed by passing
+    /// radius + D (NoiseBoundTests.WarpedPointsStayEnclosed).
+    /// </summary>
+    public static float2 iv_snoise_ball(float3 centre, float radius)
+    {
+        var n = snoise(centre);
+        var e = SNOISE_LIPSCHITZ * radius;
+        return new(max(n - e, -1.0f), min(n + e, 1.0f));
+    }
+
+    /// <summary>
+    /// Encloses the fBm value, sum of a_i * snoise(x * f_i) over the octaves, for every x within distance
+    /// radius of centre: the octave sum of iv_snoise_ball(centre * f_i, radius * |f_i|) scaled by a_i,
+    /// with frequency f_i and amplitude a_i compounded exactly as the value-and-gradient fBm compounds
+    /// them (its .w is what NoiseBoundTests.FbmBallEnclosesPoints encloses). octaves is clamped to [0, 16]
+    /// as there.
+    /// </summary>
+    public static float2 iv_fbm_ball(float3 centre, float radius, int octaves, float lacunarity, float gain)
+    {
+        octaves = clamp(octaves, 0, 16);
+        var amplitude = 1.0f;
+        var frequency = 1.0f;
+        var sum = new float2(0.0f, 0.0f);
+        for (var i = 0; i < octaves; i++)
+        {
+            sum = iv_add(sum, iv_scale(iv_snoise_ball(centre * frequency, radius * abs(frequency)), amplitude));
+            frequency *= lacunarity;
+            amplitude *= gain;
+        }
+
+        return sum;
+    }
+
+    /// <summary>
+    /// A ball (xyz centre, w radius) enclosing every point of every ray of a screen tile over the depth
+    /// segment [z0, z1], each point then moved by a warp of length at most warp: the ball to hand
+    /// iv_snoise_ball (scaled by the noise frequency) when one probe serves the whole tile. Camera frame:
+    /// a pinhole at the origin looking down +z, where the ray of slope m is the points (m z, z).
+    /// centreSlope is the tile's central slope m_c; footprintPerDepth is the largest distance of any
+    /// ray's slope from m_c, N / (sqrt(2) f) for an N x N pixel tile at focal length f pixels (the full
+    /// pixel footprint, so sub-pixel jitter is covered). The centre is (m_c z_m, z_m), z_m = (z0 + z1) / 2;
+    /// the radius is ((z1 - z0) / 2) |(m_c, 1)| + z1 footprintPerDepth + warp. Derivation: a point of slope
+    /// m_c + e (|e| &lt;= footprintPerDepth) at depth z is c + (z - z_m)(m_c, 1) + z (e, 0); the first term is
+    /// at most the first radius term, the second at most the second, and the warp adds at most warp.
+    /// The ball also carries its own float32 rounding: the radius is widened by (|c|_1 + radius) 2^-20.
+    /// With u = 2^-24, each centre component carries at most two roundings (z_m, then the product), so the
+    /// float centre is within 2u |c|_1 of the exact one; the radius carries at most seven, so the exact
+    /// radius is at most (1 + 7u) times the float one; the widening rounds once more. 2^-20 = 16u covers
+    /// 8u (|c|_1 + radius) twice over, and fma contraction only removes roundings. So every exact point of
+    /// every ray, at any depth and slope, lies in the returned float ball
+    /// (NoiseBoundTests.TileBallEnclosesEveryRaySegment, degenerate and far segments included). The caller
+    /// rotates and translates the centre into world space (the radius is unchanged) and scales both by its
+    /// noise frequency; the rounding of that transform is the caller's.
+    /// </summary>
+    public static float4 iv_frustum_ball(float2 centreSlope, float z0, float z1, float footprintPerDepth, float warp)
+    {
+        var zm = (z0 + z1) * 0.5f;
+        var axis = sqrt(centreSlope.x * centreSlope.x + centreSlope.y * centreSlope.y + 1.0f);
+        var cx = centreSlope.x * zm;
+        var cy = centreSlope.y * zm;
+        var radius = (z1 - z0) * 0.5f * axis + z1 * footprintPerDepth + warp;
+        radius += (abs(cx) + abs(cy) + abs(zm) + radius) * 9.5367431640625e-7f;
+        return new float4(cx, cy, zm, radius);
+    }
+}

packages/cultmath/src/CultMath/math.Interval.cs.meta

@@ -0,0 +1,2 @@
+fileFormatVersion: 2
+guid: 3eb757caae684b388921e28e3c87ba3a

packages/cultmath/src/CultMath/math.cs

@@ -388,6 +388,7 @@ public static partial class math
     public static float hash(float3 value) => hash(dot(value, new float3(127.1f, 311.7f, 74.7f)));
 
     public static uint asuint(float value) => (uint)BitConverter.SingleToInt32Bits(value);
+    public static float asfloat(uint value) => BitConverter.Int32BitsToSingle((int)value);
 
     // PCG hashes from Jarzynski and Olano, "Hash Functions for GPU Rendering" (JCGT 9(3), 2020). Integer-only, so
     // dxc and C# agree bit for bit. pcg3d/pcg4d carry uint bit patterns in int3/int4; float overloads hash IEEE bits.
@@ -599,6 +600,8 @@ public static partial class math
     // src/noise3D.glsl at 22434e04d7 (kernel radius^2 0.5, scale 105, permute (34x+10)x); see
     // THIRD-PARTY-NOTICES.md. Continuous across simplex-cell boundaries
     // (NoiseGradTests.SnoiseAndItsGradientAreContinuousAcrossSimplexCellBoundaries).
+    // SNOISE_LIPSCHITZ (math.Interval.cs) is measured from this kernel: a change here must re-run
+    // NoiseBoundTests.MeasureLipschitz and re-pin it.
     public static float snoise(float3 value)
     {
         const float cx = 1.0f / 6.0f;

packages/cultmath/tests/CultMath.Tests/HlslMirrorTests.cs

@@ -65,6 +65,7 @@ public sealed class HlslMirrorTests
     [Theory]
     [InlineData("CultMath.hlsl")]
     [InlineData("CultMath.Phacelle.hlsl")]
+    [InlineData("CultMath.Interval.hlsl")]
     public void UnityPackageShaderIsIdenticalToTheCanonicalShader(string name)
     {
         var root = GetCultMathRoot();

packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs

@@ -31,6 +31,8 @@ public sealed class HlslSourceCompatibilityTests
         Assert.NotNull(ShaderFunction(shader, "cultmath_fbm_grad", typeof(float3), typeof(int), typeof(float), typeof(float)));
         Assert.NotNull(ShaderFunction(shader, "cultmath_ridged_grad", typeof(float3), typeof(int), typeof(float), typeof(float)));
         Assert.NotNull(ShaderFunction(shader, "cultmath_phacelle", typeof(float3), typeof(float3), typeof(float), typeof(float)));
+        Assert.NotNull(ShaderFunction(shader, "cultmath_iv_snoise_ball", typeof(float3), typeof(float)));
+        Assert.NotNull(ShaderFunction(shader, "cultmath_iv_fbm_ball", typeof(float3), typeof(float), typeof(int), typeof(float), typeof(float)));
     }
 
     // HLSL functions carry no access modifier, so they compile as private instance methods.
@@ -131,7 +133,7 @@ public sealed class HlslSourceCompatibilityTests
         }
 
         Assert.True(mismatches.Count == 0, $"{mismatches.Count} mismatches:{Environment.NewLine}{string.Join(Environment.NewLine, mismatches.Take(12))}");
-        Assert.Equal(34, compared.Count);
+        Assert.Equal(53, compared.Count);
     }
 
     /// <summary>

packages/cultmath/tests/CultMath.Tests/IntervalTests.cs

@@ -0,0 +1,312 @@
+using CultMath;
+using Xunit;
+using static CultMath.math;
+
+namespace CultMath.Tests;
+
+/// <summary>
+/// intervals-enclose for every iv_* op (design.md, "Intervals"): for each op, intervals with seeded random
+/// centres and widths and the points inside them; f(x), evaluated pointwise in float32 by the same C# math
+/// the op mirrors, must lie in op(I) with no tolerance.
+///
+/// Stryker mutates operators, not float bounds, so the mutants that matter most here are written by hand
+/// and run by <see cref="EveryOneUlpShrinkIsCaught"/>: for every op, the lo bound raised by one ulp and the
+/// hi bound lowered by one ulp must each fail the harness. Two pairs are equivalent mutants on this CPU and
+/// are not in that list:
+///   - iv_exp's one-ulp widening, both sides: MathF.Exp is monotone over every float on the Linux runtime
+///     (MeasureMonotonicity prints 0 violations), so removing the widening encloses here. It is kept because
+///     exp is not required to be correctly rounded, and another libm or a GPU need not be monotone.
+///   - iv_smoothstep's one-ulp widening: the generic harness never meets one of smoothstep's 1-ulp
+///     non-monotone steps; <see cref="SmoothstepWideningCoversItsFloatNonMonotonicity"/> constructs one on
+///     each side and kills both mutants.
+/// </summary>
+public sealed class IntervalTests
+{
+    private readonly ITestOutputHelper output;
+
+    public IntervalTests(ITestOutputHelper output) => this.output = output;
+
+    private const int IntervalCount = 10000;
+    private const int PointsPerInterval = 16;
+
+    // One row of the shared harness: how many interval arguments the op takes, how to draw its scalar
+    // arguments, the op itself and the pointwise function it encloses. A new op is one entry here.
+    internal sealed record Op(
+        string Name,
+        int Intervals,
+        Func<System.Random, float[]> Scalars,
+        Func<float2[], float[], float2> Bound,
+        Func<float[], float[], float> Pointwise,
+        bool NonNegativeDomain = false);
+
+    private static float[] None(System.Random random) => Array.Empty<float>();
+
+    private static float Uniform(System.Random random, float lo, float hi) => lo + (hi - lo) * random.NextSingle();
+
+    internal static readonly Op[] Ops =
+    {
+        new("iv_point", 0, r => new[] { Uniform(r, -20.0f, 20.0f) }, (i, s) => iv_point(s[0]), (x, s) => s[0]),
+        new("iv_add", 2, None, (i, s) => iv_add(i[0], i[1]), (x, s) => x[0] + x[1]),
+        new("iv_sub", 2, None, (i, s) => iv_sub(i[0], i[1]), (x, s) => x[0] - x[1]),
+        new("iv_neg", 1, None, (i, s) => iv_neg(i[0]), (x, s) => -x[0]),
+        new("iv_mul", 2, None, (i, s) => iv_mul(i[0], i[1]), (x, s) => x[0] * x[1]),
+        new("iv_scale", 1, r => new[] { Uniform(r, -4.0f, 4.0f) }, (i, s) => iv_scale(i[0], s[0]), (x, s) => x[0] * s[0]),
+        new("iv_abs", 1, None, (i, s) => iv_abs(i[0]), (x, s) => abs(x[0])),
+        new("iv_min", 2, None, (i, s) => iv_min(i[0], i[1]), (x, s) => min(x[0], x[1])),
+        new("iv_max", 2, None, (i, s) => iv_max(i[0], i[1]), (x, s) => max(x[0], x[1])),
+        new("iv_sqr", 1, None, (i, s) => iv_sqr(i[0]), (x, s) => x[0] * x[0]),
+        new("iv_sqrt", 1, None, (i, s) => iv_sqrt(i[0]), (x, s) => sqrt(x[0]), NonNegativeDomain: true),
+        new("iv_exp", 1, None, (i, s) => iv_exp(i[0]), (x, s) => exp(x[0])),
+        new("iv_clamp", 1, r => new[] { Uniform(r, -25.0f, 25.0f), Uniform(r, -25.0f, 25.0f) },
+            (i, s) => iv_clamp(i[0], s[0], s[1]), (x, s) => clamp(x[0], s[0], s[1])),
+        new("iv_saturate", 1, None, (i, s) => iv_saturate(i[0]), (x, s) => saturate(x[0])),
+        new("iv_smoothstep", 1, r => new[] { Uniform(r, -25.0f, 25.0f), Uniform(r, -25.0f, 25.0f) },
+            (i, s) => iv_smoothstep(s[0], s[1], i[0]), (x, s) => smoothstep(s[0], s[1], x[0])),
+        new("iv_lerp", 2, r => new[] { Uniform(r, -0.5f, 1.5f) }, (i, s) => iv_lerp(i[0], i[1], s[0]), (x, s) => lerp(x[0], x[1], s[0])),
+    };
+
+    // Centres in [-20, 20], widths log-uniform in [1e-4, 10]. Every eighth interval is a point
+    // interval: a bound that is only reached when an argument is a point (iv_lerp's, where x appears
+    // twice) is still reached, so a one-ulp shrink of it is still caught.
+    private static float2 RandomInterval(System.Random random, int index)
+    {
+        var centre = Uniform(random, -20.0f, 20.0f);
+        if (index % 8 == 7)
+            return iv_point(centre);
+        var width = MathF.Pow(10.0f, Uniform(random, -4.0f, 1.0f));
+        return new float2(centre - width * 0.5f, centre + width * 0.5f);
+    }
+
+    // The points inside one argument: both endpoints, zero when the interval holds it, and the rest
+    // uniform. Endpoints and zero are where a tight bound is reached.
+    private static float[] PointsInside(System.Random random, float2 interval, bool nonNegative)
+    {
+        var lo = nonNegative ? max(interval.x, 0.0f) : interval.x;
+        var hi = nonNegative ? max(interval.y, 0.0f) : interval.y;
+        var points = new float[PointsPerInterval];
+        points[0] = lo;
+        points[1] = hi;
+        var next = 2;
+        if (lo < 0.0f && hi > 0.0f)
+            points[next++] = 0.0f;
+        for (; next < points.Length; next++)
+            points[next] = min(max(lo + (hi - lo) * random.NextSingle(), lo), hi);
+        return points;
+    }
+
+    /// <summary>Runs one op over the seeded set; returns the violations, at most a few, described.</summary>
+    internal static List<string> Violations(Op op, Func<float2, float2>? mutate = null, int seed = 0x1A7E)
+    {
+        var random = new System.Random(seed);
+        var violations = new List<string>();
+        var count = 0;
+        for (var k = 0; k < IntervalCount; k++)
+        {
+            var intervals = new float2[op.Intervals];
+            for (var a = 0; a < op.Intervals; a++)
+                intervals[a] = RandomInterval(random, k + a);
+            var scalars = op.Scalars(random);
+            var bound = op.Bound(intervals, scalars);
+            if (mutate is not null)
+                bound = mutate(bound);
+
+            var points = intervals.Select(i => PointsInside(random, i, op.NonNegativeDomain)).ToArray();
+            for (var p = 0; p < PointsPerInterval; p++)
+            {
+                // Binary ops visit the four endpoint corners first (p = 0..3), then pair points by index.
+                var x = new float[op.Intervals];
+                for (var a = 0; a < op.Intervals; a++)
+                    x[a] = op.Intervals == 2 && p < 4 ? points[a][a == 0 ? p & 1 : p >> 1] : points[a][p];
+                var value = op.Pointwise(x, scalars);
+                if (!(value >= bound.x && value <= bound.y))
+                {
+                    if (count++ < 4)
+                        violations.Add($"{op.Name}({string.Join(", ", intervals)}; {string.Join(", ", scalars)}) = {bound}, f({string.Join(", ", x)}) = {value:R}");
+                }
+            }
+        }
+
+        if (count > violations.Count)
+            violations.Add($"... {count} violations in all");
+        return violations;
+    }
+
+    [Fact]
+    public void EveryOpEnclosesItsPoints()
+    {
+        var failures = Ops.SelectMany(op => Violations(op)).ToList();
+        Assert.True(failures.Count == 0, string.Join(Environment.NewLine, failures));
+    }
+
+    internal static float NextUp(float v) => v == 0.0f ? asfloat(1u) : asfloat(v > 0.0f ? asuint(v) + 1u : asuint(v) - 1u);
+    private static float NextDown(float v) => -NextUp(-v);
+
+    /// <summary>
+    /// The hand-mutation list: every op's lo raised by one ulp, and its hi lowered by one ulp, must each
+    /// fail the harness. The two equivalent pairs are named in the class comment.
+    /// </summary>
+    [Fact]
+    public void EveryOneUlpShrinkIsCaught()
+    {
+        var survivors = new List<string>();
+        foreach (var op in Ops.Where(o => o.Name is not "iv_exp" and not "iv_smoothstep"))
+        {
+            if (Violations(op, b => new float2(NextUp(b.x), b.y)).Count == 0)
+                survivors.Add(op.Name + " lo + 1 ulp");
+            if (Violations(op, b => new float2(b.x, NextDown(b.y))).Count == 0)
+                survivors.Add(op.Name + " hi - 1 ulp");
+        }
+
+        Assert.True(survivors.Count == 0, "surviving one-ulp shrinks: " + string.Join(", ", survivors));
+    }
+
+    [Fact]
+    public void MulCoversEverySignCase()
+    {
+        // Both positive, both negative, mixed both ways, and a factor straddling zero on either side.
+        Assert.Equal(new float2(2.0f, 12.0f), iv_mul(new float2(1.0f, 3.0f), new float2(2.0f, 4.0f)));
+        Assert.Equal(new float2(2.0f, 12.0f), iv_mul(new float2(-3.0f, -1.0f), new float2(-4.0f, -2.0f)));
+        Assert.Equal(new float2(-12.0f, -2.0f), iv_mul(new float2(1.0f, 3.0f), new float2(-4.0f, -2.0f)));
+        Assert.Equal(new float2(-12.0f, -2.0f), iv_mul(new float2(-3.0f, -1.0f), new float2(2.0f, 4.0f)));
+        Assert.Equal(new float2(-8.0f, 12.0f), iv_mul(new float2(-2.0f, 3.0f), new float2(2.0f, 4.0f)));
+        Assert.Equal(new float2(-12.0f, 8.0f), iv_mul(new float2(-2.0f, 3.0f), new float2(-4.0f, -2.0f)));
+        Assert.Equal(new float2(-15.0f, 12.0f), iv_mul(new float2(-2.0f, 3.0f), new float2(-5.0f, 4.0f)));
+    }
+
+    [Fact]
+    public void SqrtClampsANegativeLoAtZero()
+    {
+        Assert.Equal(new float2(0.0f, 3.0f), iv_sqrt(new float2(-4.0f, 9.0f)));
+        Assert.Equal(new float2(0.0f, 0.0f), iv_sqrt(new float2(-4.0f, -1.0f)));
+        Assert.Equal(new float2(2.0f, 3.0f), iv_sqrt(new float2(4.0f, 9.0f)));
+    }
+
+    [Fact]
+    public void AbsAndSqrStraddlingZeroStartAtZero()
+    {
+        Assert.Equal(new float2(0.0f, 3.0f), iv_abs(new float2(-3.0f, 2.0f)));
+        Assert.Equal(new float2(0.0f, 5.0f), iv_abs(new float2(-2.0f, 5.0f)));
+        Assert.Equal(new float2(2.0f, 5.0f), iv_abs(new float2(-5.0f, -2.0f)));
+        Assert.Equal(new float2(0.0f, 9.0f), iv_sqr(new float2(-3.0f, 2.0f)));
+        Assert.Equal(new float2(4.0f, 25.0f), iv_sqr(new float2(-5.0f, -2.0f)));
+    }
+
+    [Fact]
+    public void MinAndMaxAreThePointwiseMinAndMaxOfTheEnds()
+    {
+        var a = new float2(-1.0f, 4.0f);
+        var b = new float2(2.0f, 3.0f);
+        Assert.Equal(new float2(-1.0f, 3.0f), iv_min(a, b));
+        Assert.Equal(new float2(2.0f, 4.0f), iv_max(a, b));
+    }
+
+    [Fact]
+    public void SmoothstepBoundsAreItsOrderedEndpointsWidenedByOneUlp()
+    {
+        var rising = iv_smoothstep(0.0f, 1.0f, new float2(0.25f, 0.5f));
+        Assert.Equal(NextDown(smoothstep(0.0f, 1.0f, 0.25f)), rising.x);
+        Assert.Equal(NextUp(smoothstep(0.0f, 1.0f, 0.5f)), rising.y);
+
+        // minimum > maximum: smoothstep falls, so the endpoints swap.
+        var falling = iv_smoothstep(1.0f, 0.0f, new float2(0.25f, 0.5f));
+        Assert.Equal(NextDown(smoothstep(1.0f, 0.0f, 0.5f)), falling.x);
+        Assert.Equal(NextUp(smoothstep(1.0f, 0.0f, 0.25f)), falling.y);
+
+        // The widening never leaves [0, 1].
+        Assert.Equal(new float2(0.0f, 1.0f), iv_smoothstep(0.0f, 1.0f, new float2(-2.0f, 3.0f)));
+    }
+
+    /// <summary>
+    /// smoothstep's float32 evaluation steps down by one ulp in places (MeasureMonotonicity). Three
+    /// consecutive floats t0 &lt; t1 &lt; t2 with f(t1) below both neighbours put f(t1) outside the ordered
+    /// endpoints of [t0, t2]; three with f(t1) above both put it outside on the other side. The widened
+    /// bound encloses both; without the widening it does not.
+    /// </summary>
+    [Fact]
+    public void SmoothstepWideningCoversItsFloatNonMonotonicity()
+    {
+        float F(float t) => smoothstep(0.0f, 1.0f, t);
+        float2? dip = null, peak = null;
+        for (var t = 1.0e-4f; t < 0.01f && (dip is null || peak is null); t = NextUp(t))
+        {
+            var t1 = NextUp(t);
+            var t2 = NextUp(t1);
+            if (dip is null && F(t1) < F(t) && F(t1) < F(t2))
+                dip = new float2(t, t2);
+            if (peak is null && F(t1) > F(t) && F(t1) > F(t2))
+                peak = new float2(t, t2);
+        }
+
+        Assert.NotNull(dip);
+        Assert.NotNull(peak);
+        foreach (var (interval, low) in new[] { (dip!.Value, true), (peak!.Value, false) })
+        {
+            var mid = F(NextUp(interval.x));
+            var unwidened = new float2(min(F(interval.x), F(interval.y)), max(F(interval.x), F(interval.y)));
+            Assert.False(mid >= unwidened.x && mid <= unwidened.y, $"the {(low ? "dip" : "peak")} at {interval} is inside the endpoints");
+            var bound = iv_smoothstep(0.0f, 1.0f, interval);
+            Assert.True(mid >= bound.x && mid <= bound.y, $"iv_smoothstep{interval} = {bound} misses {mid:R}");
+        }
+    }
+
+    [Fact]
+    public void LerpIsTheNaturalExtensionAndTightForPointArguments()
+    {
+        Assert.Equal(new float2(lerp(1.0f, 3.0f, 0.25f), lerp(1.0f, 7.0f, 0.25f)), iv_lerp(iv_point(1.0f), new float2(3.0f, 7.0f), 0.25f));
+        Assert.Equal(new float2(-1.0f + (2.0f - 1.0f) * 0.5f, 1.0f + (4.0f + 1.0f) * 0.5f), iv_lerp(new float2(-1.0f, 1.0f), new float2(2.0f, 4.0f), 0.5f));
+    }
+
+    [Fact]
+    public void ExpIsWidenedByOneUlpOnEachSide()
+    {
+        var bound = iv_exp(new float2(-1.0f, 2.0f));
+        Assert.Equal(NextDown(exp(-1.0f)), bound.x);
+        Assert.Equal(NextUp(exp(2.0f)), bound.y);
+        Assert.Equal(0.0f, iv_exp(new float2(-200.0f, 0.0f)).x);
+    }
+
+    /// <summary>
+    /// The monotonicity the widening rules rest on, measured over every float: exp over [-87, 88.7], the
+    /// range where it neither underflows to a subnormal nor overflows, and smoothstep(0, 1, t) over [0, 1],
+    /// which covers every smoothstep because t = saturate((x - minimum) / (maximum - minimum)) is itself
+    /// monotone in x. Prints the number of steps down and the largest, in ulps. Slow: about a minute.
+    /// </summary>
+    [Fact(Explicit = true)]
+    [Trait("Category", "Slow")]
+    public void MeasureMonotonicity()
+    {
+        long expDown = 0, expWorst = 0;
+        var previous = exp(-87.0f);
+        for (var x = NextUp(-87.0f); x <= 88.7f; x = NextUp(x))
+        {
+            var e = exp(x);
+            if (e < previous)
+            {
+                expDown++;
+                expWorst = Math.Max(expWorst, (long)asuint(previous) - asuint(e));
+            }
+
+            previous = e;
+        }
+
+        long smoothDown = 0, smoothWorst = 0;
+        var runningMax = 0.0f;
+        for (var bits = 0u; bits <= asuint(1.0f); bits++)
+        {
+            var s = smoothstep(0.0f, 1.0f, asfloat(bits));
+            if (s < runningMax)
+            {
+                smoothDown++;
+                smoothWorst = Math.Max(smoothWorst, (long)asuint(runningMax) - asuint(s));
+            }
+            else
+            {
+                runningMax = s;
+            }
+        }
+
+        output.WriteLine($"IV-REPORT monotonicity: exp steps down {expDown} times, worst {expWorst} ulp; smoothstep steps down {smoothDown} times, worst {smoothWorst} ulp");
+        Assert.True(expWorst <= 1, $"exp steps down by {expWorst} ulp; one ulp of widening is not enough");
+        Assert.True(smoothWorst <= 1, $"smoothstep steps down by {smoothWorst} ulp; one ulp of widening is not enough");
+    }
+}

packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl

@@ -0,0 +1,122 @@
+#ifndef CULTMATH_INTERVAL_HLSL
+#define CULTMATH_INTERVAL_HLSL
+
+// Interval arithmetic over float2(lo, hi), lo <= hi, both finite. Every cultmath_iv_* function returns
+// an interval containing f(x) for every x in its input intervals, mirroring math.Interval.cs bit for
+// bit; see that file and docs/design.md ("Intervals") for each op's enclosure and its ulp widening.
+// Included by CultMath.hlsl after its own definitions of cultmath_smoothstep and cultmath_snoise.
+
+// Lipschitz constant of cultmath_snoise(float3); provenance in math.Interval.cs (SNOISE_LIPSCHITZ).
+static const float CULTMATH_SNOISE_LIPSCHITZ = 10.099261;
+
+float2 cultmath_iv_point(float value) { return float2(value, value); }
+float2 cultmath_iv_add(float2 a, float2 b) { return float2(a.x + b.x, a.y + b.y); }
+float2 cultmath_iv_sub(float2 a, float2 b) { return float2(a.x - b.y, a.y - b.x); }
+float2 cultmath_iv_neg(float2 a) { return float2(-a.y, -a.x); }
+
+float2 cultmath_iv_mul(float2 a, float2 b)
+{
+    float p0 = a.x * b.x;
+    float p1 = a.x * b.y;
+    float p2 = a.y * b.x;
+    float p3 = a.y * b.y;
+    return float2(min(min(p0, p1), min(p2, p3)), max(max(p0, p1), max(p2, p3)));
+}
+
+float2 cultmath_iv_scale(float2 a, float s) { return s >= 0.0 ? float2(a.x * s, a.y * s) : float2(a.y * s, a.x * s); }
+
+float2 cultmath_iv_abs(float2 a)
+{
+    if (a.x >= 0.0)
+        return a;
+    if (a.y <= 0.0)
+        return float2(-a.y, -a.x);
+    return float2(0.0, max(-a.x, a.y));
+}
+
+float2 cultmath_iv_min(float2 a, float2 b) { return float2(min(a.x, b.x), min(a.y, b.y)); }
+float2 cultmath_iv_max(float2 a, float2 b) { return float2(max(a.x, b.x), max(a.y, b.y)); }
+
+float2 cultmath_iv_sqr(float2 a)
+{
+    float l = a.x * a.x;
+    float h = a.y * a.y;
+    if (a.x >= 0.0)
+        return float2(l, h);
+    if (a.y <= 0.0)
+        return float2(h, l);
+    return float2(0.0, max(l, h));
+}
+
+float2 cultmath_iv_sqrt(float2 a) { return float2(sqrt(max(a.x, 0.0)), sqrt(max(a.y, 0.0))); }
+
+// Widened by one ulp on each side: exp is not required to be correctly rounded.
+float2 cultmath_iv_exp(float2 a)
+{
+    float lo = exp(a.x);
+    float hi = exp(a.y);
+    return float2(lo > 0.0 ? asfloat(asuint(lo) - 1u) : 0.0, asfloat(asuint(hi) + 1u));
+}
+
+float2 cultmath_iv_clamp(float2 a, float minimum, float maximum)
+{
+    return float2(cultmath_clamp(a.x, minimum, maximum), cultmath_clamp(a.y, minimum, maximum));
+}
+
+float2 cultmath_iv_saturate(float2 a) { return float2(cultmath_saturate(a.x), cultmath_saturate(a.y)); }
+
+// The endpoints, ordered, widened by one ulp inside [0, 1]: the float32 evaluation of smoothstep is
+// monotone only to within one ulp.
+float2 cultmath_iv_smoothstep(float minimum, float maximum, float2 a)
+{
+    float s0 = cultmath_smoothstep(minimum, maximum, a.x);
+    float s1 = cultmath_smoothstep(minimum, maximum, a.y);
+    float lo = min(s0, s1);
+    float hi = max(s0, s1);
+    return float2(lo > 0.0 ? asfloat(asuint(lo) - 1u) : 0.0, hi < 1.0 ? asfloat(asuint(hi) + 1u) : 1.0);
+}
+
+// The natural interval extension of x + (y - x) * amount, which is how cultmath_lerp evaluates.
+float2 cultmath_iv_lerp(float2 a, float2 b, float amount) { return cultmath_iv_add(a, cultmath_iv_scale(cultmath_iv_sub(b, a), amount)); }
+
+// [n - L r, n + L r] intersected with [-1, 1], n = cultmath_snoise(centre). One snoise, value only.
+float2 cultmath_iv_snoise_ball(float3 centre, float radius)
+{
+    float n = cultmath_snoise(centre);
+    float e = CULTMATH_SNOISE_LIPSCHITZ * radius;
+    return float2(max(n - e, -1.0), min(n + e, 1.0));
+}
+
+// Octave sum of cultmath_iv_snoise_ball, compounding frequency and amplitude as the value-and-gradient fBm does.
+float2 cultmath_iv_fbm_ball(float3 centre, float radius, int octaves, float lacunarity, float gain)
+{
+    octaves = clamp(octaves, 0, 16);
+    float amplitude = 1.0;
+    float frequency = 1.0;
+    float2 sum = float2(0.0, 0.0);
+    for (int i = 0; i < octaves; i++)
+    {
+        sum = cultmath_iv_add(sum, cultmath_iv_scale(cultmath_iv_snoise_ball(centre * frequency, radius * abs(frequency)), amplitude));
+        frequency *= lacunarity;
+        amplitude *= gain;
+    }
+
+    return sum;
+}
+
+// A ball (xyz centre, w radius) enclosing every point of every ray of a screen tile over the depth segment
+// [z0, z1], warped by at most warp. Camera frame: pinhole at the origin looking down +z, ray (m z, z).
+// footprintPerDepth = N / (sqrt(2) f) for an N x N tile at focal length f. The radius is widened by
+// (|c|_1 + radius) 2^-20, which covers the ball's own float32 rounding. Derivation in math.Interval.cs.
+float4 cultmath_iv_frustum_ball(float2 centreSlope, float z0, float z1, float footprintPerDepth, float warp)
+{
+    float zm = (z0 + z1) * 0.5;
+    float axis = sqrt(centreSlope.x * centreSlope.x + centreSlope.y * centreSlope.y + 1.0);
+    float cx = centreSlope.x * zm;
+    float cy = centreSlope.y * zm;
+    float radius = (z1 - z0) * 0.5 * axis + z1 * footprintPerDepth + warp;
+    radius += (abs(cx) + abs(cy) + abs(zm) + radius) * 9.5367431640625e-7;
+    return float4(cx, cy, zm, radius);
+}
+
+#endif

packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl.meta

@@ -0,0 +1,7 @@
+fileFormatVersion: 2
+guid: 18b031192f344246ac523cea014b8f24
+ShaderIncludeImporter:
+  externalObjects: {}
+  userData: 
+  assetBundleName: 
+  assetBundleVariant: 

packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl

@@ -504,6 +504,7 @@ CultCellular cultmath_cellular(float3 p)
 }
 
 #include "CultMath.Phacelle.hlsl"
+#include "CultMath.Interval.hlsl"
 
 float cultmath_value_noise(float2 position)
 {