diff --git a/packages/cultmath/README.md b/packages/cultmath/README.md
index f70df445..81d15859 100644
--- a/packages/cultmath/README.md
+++ b/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.
diff --git a/packages/cultmath/docs/design.md b/packages/cultmath/docs/design.md
index 0de8eb14..d102353b 100644
--- a/packages/cultmath/docs/design.md
+++ b/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`.
diff --git a/packages/cultmath/scripts/build-unity-package.ps1 b/packages/cultmath/scripts/build-unity-package.ps1
index 4af5862e..7fb6f7a6 100644
--- a/packages/cultmath/scripts/build-unity-package.ps1
+++ b/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)) {
diff --git a/packages/cultmath/shaders/CultMath.Interval.hlsl b/packages/cultmath/shaders/CultMath.Interval.hlsl
new file mode 100644
index 00000000..b6b16665
--- /dev/null
+++ b/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
diff --git a/packages/cultmath/shaders/CultMath.Interval.hlsl.meta b/packages/cultmath/shaders/CultMath.Interval.hlsl.meta
new file mode 100644
index 00000000..a22ec627
--- /dev/null
+++ b/packages/cultmath/shaders/CultMath.Interval.hlsl.meta
@@ -0,0 +1,7 @@
+fileFormatVersion: 2
+guid: 4d8a899561654ac990cab0bf48658cc1
+ShaderIncludeImporter:
+ externalObjects: {}
+ userData:
+ assetBundleName:
+ assetBundleVariant:
diff --git a/packages/cultmath/shaders/CultMath.hlsl b/packages/cultmath/shaders/CultMath.hlsl
index af2d3741..da3362a6 100644
--- a/packages/cultmath/shaders/CultMath.hlsl
+++ b/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)
{
diff --git a/packages/cultmath/src/CultMath/CultMath.csproj b/packages/cultmath/src/CultMath/CultMath.csproj
index 0b1ee09f..2a82f8db 100644
--- a/packages/cultmath/src/CultMath/CultMath.csproj
+++ b/packages/cultmath/src/CultMath/CultMath.csproj
@@ -24,6 +24,7 @@
+
diff --git a/packages/cultmath/src/CultMath/math.Interval.cs b/packages/cultmath/src/CultMath/math.Interval.cs
new file mode 100644
index 00000000..91e99f87
--- /dev/null
+++ b/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;
+
+ /// The point interval [x, x].
+ public static float2 iv_point(float value) => new(value, value);
+
+ /// Encloses x + y for x in a, y in b. Monotone: no widening.
+ public static float2 iv_add(float2 a, float2 b) => new(a.x + b.x, a.y + b.y);
+
+ /// Encloses x - y for x in a, y in b. Monotone: no widening.
+ public static float2 iv_sub(float2 a, float2 b) => new(a.x - b.y, a.y - b.x);
+
+ /// Encloses -x for x in a. Exact.
+ public static float2 iv_neg(float2 a) => new(-a.y, -a.x);
+
+ ///
+ /// 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.
+ ///
+ 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)));
+ }
+
+ /// Encloses x * s for x in a. Monotone: no widening.
+ 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);
+
+ /// Encloses abs(x) for x in a; an interval straddling zero gives [0, max(-lo, hi)]. Exact.
+ 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));
+ }
+
+ /// Encloses min(x, y) for x in a, y in b. Monotone: no widening.
+ public static float2 iv_min(float2 a, float2 b) => new(min(a.x, b.x), min(a.y, b.y));
+
+ /// Encloses max(x, y) for x in a, y in b. Monotone: no widening.
+ public static float2 iv_max(float2 a, float2 b) => new(max(a.x, b.x), max(a.y, b.y));
+
+ ///
+ /// 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.
+ ///
+ 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));
+ }
+
+ ///
+ /// 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.
+ ///
+ public static float2 iv_sqrt(float2 a) => new(sqrt(max(a.x, 0.0f)), sqrt(max(a.y, 0.0f)));
+
+ ///
+ /// 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 < cutoff reads as "not provably below", the safe side.
+ ///
+ 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));
+ }
+
+ /// Encloses clamp(x, minimum, maximum) for x in a. Monotone: no widening.
+ public static float2 iv_clamp(float2 a, float minimum, float maximum) =>
+ new(clamp(a.x, minimum, maximum), clamp(a.y, minimum, maximum));
+
+ /// Encloses saturate(x) for x in a. Monotone: no widening.
+ public static float2 iv_saturate(float2 a) => new(saturate(a.x), saturate(a.y));
+
+ ///
+ /// Encloses smoothstep(minimum, maximum, x) for x in a. smoothstep is monotone in x (increasing when
+ /// minimum < maximum, decreasing when minimum > 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].
+ ///
+ 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);
+ }
+
+ ///
+ /// 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.
+ ///
+ public static float2 iv_lerp(float2 a, float2 b, float amount) => iv_add(a, iv_scale(iv_sub(b, a), amount));
+
+ ///
+ /// 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).
+ ///
+ 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));
+ }
+
+ ///
+ /// 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.
+ ///
+ 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;
+ }
+
+ ///
+ /// 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| <= 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.
+ ///
+ 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);
+ }
+}
diff --git a/packages/cultmath/src/CultMath/math.Interval.cs.meta b/packages/cultmath/src/CultMath/math.Interval.cs.meta
new file mode 100644
index 00000000..eb3b76d1
--- /dev/null
+++ b/packages/cultmath/src/CultMath/math.Interval.cs.meta
@@ -0,0 +1,2 @@
+fileFormatVersion: 2
+guid: 3eb757caae684b388921e28e3c87ba3a
diff --git a/packages/cultmath/src/CultMath/math.cs b/packages/cultmath/src/CultMath/math.cs
index 23916b6b..d6c85847 100644
--- a/packages/cultmath/src/CultMath/math.cs
+++ b/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;
diff --git a/packages/cultmath/tests/CultMath.Tests/HlslMirrorTests.cs b/packages/cultmath/tests/CultMath.Tests/HlslMirrorTests.cs
index 6b69a908..9e2dc9af 100644
--- a/packages/cultmath/tests/CultMath.Tests/HlslMirrorTests.cs
+++ b/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();
diff --git a/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs b/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs
index 2cb8ae49..a5965a3d 100644
--- a/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs
+++ b/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);
}
///
diff --git a/packages/cultmath/tests/CultMath.Tests/IntervalTests.cs b/packages/cultmath/tests/CultMath.Tests/IntervalTests.cs
new file mode 100644
index 00000000..4f165eab
--- /dev/null
+++ b/packages/cultmath/tests/CultMath.Tests/IntervalTests.cs
@@ -0,0 +1,312 @@
+using CultMath;
+using Xunit;
+using static CultMath.math;
+
+namespace CultMath.Tests;
+
+///
+/// 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 : 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; constructs one on
+/// each side and kills both mutants.
+///
+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 Scalars,
+ Func Bound,
+ Func Pointwise,
+ bool NonNegativeDomain = false);
+
+ private static float[] None(System.Random random) => Array.Empty();
+
+ 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;
+ }
+
+ /// Runs one op over the seeded set; returns the violations, at most a few, described.
+ internal static List Violations(Op op, Func? mutate = null, int seed = 0x1A7E)
+ {
+ var random = new System.Random(seed);
+ var violations = new List();
+ 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);
+
+ ///
+ /// 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.
+ ///
+ [Fact]
+ public void EveryOneUlpShrinkIsCaught()
+ {
+ var survivors = new List();
+ 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)));
+ }
+
+ ///
+ /// smoothstep's float32 evaluation steps down by one ulp in places (MeasureMonotonicity). Three
+ /// consecutive floats t0 < t1 < 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.
+ ///
+ [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);
+ }
+
+ ///
+ /// 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.
+ ///
+ [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");
+ }
+}
diff --git a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs
new file mode 100644
index 00000000..2b00ca88
--- /dev/null
+++ b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs
@@ -0,0 +1,1596 @@
+using CultMath;
+using Xunit;
+using static CultMath.math;
+
+namespace CultMath.Tests;
+
+///
+/// The noise bounds (design.md, "Intervals"): iv_snoise_ball and iv_fbm_ball enclose their functions over
+/// a ball, a domain warp bounded by D is enclosed by enlarging the radius by D, SNOISE_LIPSCHITZ has
+/// measured provenance, iv_frustum_ball encloses a screen tile's rays over a depth segment, the bound
+/// composed from an analytic envelope and the noise ball encloses the density, and a tile march over
+/// fields with authored empty space is measured against a dense one for the same transmittance. Lines
+/// starting IV-REPORT carry the numbers the cut report quotes.
+///
+/// Hand mutations of iv_frustum_ball that TileBallEnclosesEveryRaySegment must kill: z0 for z1 in the
+/// lateral term (z0 * footprintPerDepth); |m_c| for |(m_c, 1)| in the depth term; the warp dropped; the
+/// rounding widening dropped (which the degenerate family's enclosure check alone must catch). Hand
+/// mutations of the test's own fields and marches: w(z0) for w(z1) in VoidField.WeightBound must fail
+/// WeightBoundEnclosesOneSignedNoise; CarveRadius for HollowRadius in the lo test of VoidField.Compose
+/// and BodyDepth without |dir| must fail AnalyticBodyMatchesFineMarch; the footprint march's pre-body
+/// integration dropped, or LodStep doubled, must fail IntervalSkipHalvesEvaluations's depth agreement.
+///
+public sealed class NoiseBoundTests
+{
+ private readonly ITestOutputHelper output;
+
+ public NoiseBoundTests(ITestOutputHelper output) => this.output = output;
+
+ private static float Uniform(System.Random random, float lo, float hi) => lo + (hi - lo) * random.NextSingle();
+
+ private static float LogUniform(System.Random random, float lo, float hi) =>
+ MathF.Exp(Uniform(random, MathF.Log(lo), MathF.Log(hi)));
+
+ private static float3 UnitVector(System.Random random)
+ {
+ float3 v;
+ do
+ {
+ v = new float3(Uniform(random, -1.0f, 1.0f), Uniform(random, -1.0f, 1.0f), Uniform(random, -1.0f, 1.0f));
+ }
+ while (dot(v, v) > 1.0f || dot(v, v) < 1.0e-4f);
+
+ return normalize(v);
+ }
+
+ // The first eight points of a ball are on its surface, where the bound is closest to binding; the
+ // rest are uniform in its volume.
+ private static float3 PointInBall(System.Random random, float3 centre, float radius, int index) =>
+ centre + UnitVector(random) * (index < 8 ? radius : radius * MathF.Cbrt(random.NextSingle()));
+
+ private static float3 RandomCentre(System.Random random) =>
+ new(Uniform(random, -20.0f, 20.0f), Uniform(random, -20.0f, 20.0f), Uniform(random, -20.0f, 20.0f));
+
+ private static bool Inside(float value, float2 bound) => value >= bound.x && value <= bound.y;
+
+ [Fact]
+ public void SnoiseBallEnclosesPoints()
+ {
+ var random = new System.Random(0xBA11);
+ for (var b = 0; b < 2000; b++)
+ {
+ var centre = RandomCentre(random);
+ var radius = LogUniform(random, 1.0e-3f, 1.0e3f);
+ var bound = iv_snoise_ball(centre, radius);
+ for (var i = 0; i < 64; i++)
+ {
+ var x = PointInBall(random, centre, radius, i);
+ Assert.True(Inside(snoise(x), bound), $"iv_snoise_ball({centre}, {radius:R}) = {bound} misses snoise({x}) = {snoise(x):R}");
+ }
+ }
+ }
+
+ [Fact]
+ public void FbmBallEnclosesPoints()
+ {
+ var random = new System.Random(0xFB11);
+ for (var b = 0; b < 2000; b++)
+ {
+ var centre = RandomCentre(random);
+ var radius = LogUniform(random, 1.0e-3f, 1.0e3f);
+ var octaves = random.Next(0, 7);
+ var lacunarity = Uniform(random, 1.5f, 2.5f);
+ var gain = Uniform(random, 0.3f, 0.7f);
+ var bound = iv_fbm_ball(centre, radius, octaves, lacunarity, gain);
+ for (var i = 0; i < 64; i++)
+ {
+ var x = PointInBall(random, centre, radius, i);
+ var value = fbm_grad(x, octaves, lacunarity, gain).w;
+ Assert.True(Inside(value, bound), $"iv_fbm_ball({centre}, {radius:R}, {octaves}, {lacunarity:R}, {gain:R}) = {bound} misses {value:R} at {x}");
+ }
+ }
+ }
+
+ ///
+ /// The site's march rests on this: a point x in the ball (c, r), displaced by a flow of length at most
+ /// D, lies in the ball (c, r + D), so snoise of the warped point lies in iv_snoise_ball(c, r + D). The
+ /// geometric check allows one part in 1e5 for the rounding of x + flow.
+ ///
+ [Fact]
+ public void WarpedPointsStayEnclosed()
+ {
+ var random = new System.Random(0x3A4B);
+ for (var b = 0; b < 2000; b++)
+ {
+ var centre = RandomCentre(random);
+ var radius = LogUniform(random, 1.0e-3f, 2.0f);
+ var warp = Uniform(random, 0.0f, 1.0f);
+ var bound = iv_snoise_ball(centre, radius + warp);
+ for (var i = 0; i < 64; i++)
+ {
+ var x = PointInBall(random, centre, radius, i);
+ var flow = UnitVector(random) * (warp * MathF.Sqrt(random.NextSingle()));
+ var warped = x + flow;
+ Assert.True(length(warped - centre) <= (radius + warp) * (1.0f + 1.0e-5f), $"warped point {warped} leaves the ball ({centre}, {radius + warp:R})");
+ Assert.True(Inside(snoise(warped), bound), $"iv_snoise_ball({centre}, {radius + warp:R}) = {bound} misses snoise({warped}) = {snoise(warped):R}");
+ }
+ }
+ }
+
+ // ---- The tile ball and the composed bound ----
+
+ // One tile case of TileBallEnclosesEveryRaySegment: the r2 domain, or one of the extreme families
+ // where float32 rounding of the centre is largest against the radius (Soul's verdict s1).
+ private static (int N, float Focal, float2 Slope, float Z0, float Z1, float Warp) DrawBallCase(System.Random random, string family)
+ {
+ int[] sizes = { 1, 4, 8, 16 };
+ float[] depths = { 12.0f, 1000.0f, 2400.0f };
+ var n = sizes[random.Next(sizes.Length)];
+ var focal = Uniform(random, 500.0f, 3000.0f);
+ var slope = new float2(Uniform(random, -0.77f, 0.77f), Uniform(random, -0.77f, 0.77f));
+ var depth = depths[random.Next(depths.Length)];
+ var z0 = random.Next(4) == 0 ? 0.0f : Uniform(random, 0.0f, depth * 0.9f);
+ var z1 = z0 + LogUniform(random, 1.0e-3f * depth, depth - z0);
+ var warp = random.Next(4) == 0 ? 0.0f : LogUniform(random, 1.0e-3f, 60.0f);
+ float Signed(float magnitude) => random.Next(2) == 0 ? -magnitude : magnitude;
+ switch (family)
+ {
+ case "degenerate":
+ z0 = z1 = LogUniform(random, 1.0e-3f * depth, depth);
+ break;
+ case "thin":
+ z0 = z1 - 1.0e-6f * z1;
+ break;
+ case "wide slope":
+ slope = new float2(Signed(LogUniform(random, 1.0e-3f, 707.0f)), Signed(LogUniform(random, 1.0e-3f, 707.0f)));
+ break;
+ case "far":
+ z1 = LogUniform(random, 1.0e4f, 1.0e7f);
+ z0 = z1 - z1 * LogUniform(random, 1.0e-6f, 1.0f);
+ break;
+ case "from the camera":
+ z0 = 0.0f;
+ z1 = LogUniform(random, 1.0e-3f, 1.0f);
+ break;
+ case "tiny focal":
+ focal = LogUniform(random, 1.0f, 500.0f);
+ break;
+ case "huge warp":
+ warp = LogUniform(random, 60.0f, 1.0e6f);
+ break;
+ }
+
+ return (n, focal, slope, z0, z1, warp);
+ }
+
+ ///
+ /// iv_frustum_ball encloses every exact point of every ray of an N x N tile over a depth segment,
+ /// sub-pixel jitter included, each point then moved by a flow of length at most warp, with no
+ /// tolerance: the distance is taken in double from the float32 ball to the point (m z, z) formed in
+ /// double from the float32 slope and depth. 2,000 seeded tiles of the r2 domain (N in {1, 4, 8, 16},
+ /// f in [500, 3000], |m_c| up to 1.1, segments over the slab's, the fog's and the void's depth ranges,
+ /// a quarter starting at the camera) x 64 points, and 2,000 tiles of each extreme family (z0 == z1;
+ /// z1 - z0 = 1e-6 z1; |m_c| up to 1000; z up to 1e7; z0 = 0, z1 down to 1e-3; f down to 1; warp up to
+ /// 1e6) x the 8 points at the footprint's corners at both depth ends. Every corner point's flow is
+ /// full-length and points outward from the ball's centre. It also pins the ball to the derivation in
+ /// double, both ways: the centre within 2^-22 |c|_1, and r <= ball.w <= r + 2^-19 (r + |c|_1), so a
+ /// ball that forgets its rounding widening fails and one that grows past it fails. Misses and pin
+ /// failures are counted per family and reported together. Each tile's camera is drawn on its own, as a
+ /// moving camera's frames are: no ball, mask or probe result is carried from one to the next.
+ ///
+ [Fact]
+ public void TileBallEnclosesEveryRaySegment()
+ {
+ var random = new System.Random(0x711E);
+ string[] families = { "r2 domain", "degenerate", "thin", "wide slope", "far", "from the camera", "tiny focal", "huge warp" };
+ var misses = new Dictionary();
+ var pins = new Dictionary();
+ var witness = new List();
+ foreach (var family in families)
+ {
+ misses[family] = 0;
+ pins[family] = 0;
+ for (var t = 0; t < 2000; t++)
+ {
+ var (n, focal, slope, z0, z1, warp) = DrawBallCase(random, family);
+ var half = n / (2.0f * focal);
+ var footprint = n / (MathF.Sqrt(2.0f) * focal);
+ var ball = iv_frustum_ball(slope, z0, z1, footprint, warp);
+ var (cx, cy, cz, w) = ((double)ball.x, (double)ball.y, (double)ball.z, (double)ball.w);
+
+ var zm = (z0 + (double)z1) / 2.0;
+ var (ex, ey) = (slope.x * zm, slope.y * zm);
+ var c1 = Math.Abs(ex) + Math.Abs(ey) + Math.Abs(zm);
+ var radius = (z1 - (double)z0) / 2.0 * Math.Sqrt((double)slope.x * slope.x + (double)slope.y * slope.y + 1.0) + z1 * (double)footprint + warp;
+ var centreOff = Math.Abs(cx - ex) + Math.Abs(cy - ey) + Math.Abs(cz - zm);
+ if (!(centreOff <= Math.ScaleB(c1, -22) && radius <= w && w <= radius + Math.ScaleB(radius + c1, -19)))
+ {
+ pins[family]++;
+ if (witness.Count < 8)
+ witness.Add($"{family} pin: iv_frustum_ball({slope}, {z0:R}, {z1:R}, {footprint:R}, {warp:R}) = {ball}; derivation centre ({ex:R}, {ey:R}, {zm:R}), radius {radius:R}");
+ }
+
+ for (var i = 0; i < (family == "r2 domain" ? 64 : 8); i++)
+ {
+ var offset = i < 8
+ ? new float2((i & 1) == 0 ? -half : half, (i & 2) == 0 ? -half : half)
+ : new float2(Uniform(random, -half, half), Uniform(random, -half, half));
+ var z = i < 8 ? ((i & 4) == 0 ? z0 : z1) : Uniform(random, z0, z1);
+ var m = slope + offset;
+ var (px, py, pz) = (m.x * (double)z, m.y * (double)z, (double)z);
+ double fx, fy, fz;
+ if (i < 8)
+ {
+ // Outward: along the corner's direction from the centre, full length.
+ var (dx, dy, dz) = (px - cx, py - cy, pz - cz);
+ var away = Math.Sqrt(dx * dx + dy * dy + dz * dz);
+ (fx, fy, fz) = away > 0.0 ? (dx / away * warp, dy / away * warp, dz / away * warp) : (warp, 0.0, 0.0);
+ }
+ else
+ {
+ var flow = UnitVector(random) * (warp * MathF.Sqrt(random.NextSingle()));
+ (fx, fy, fz) = (flow.x, flow.y, flow.z);
+ }
+
+ var (qx, qy, qz) = (px + fx - cx, py + fy - cy, pz + fz - cz);
+ var distance = Math.Sqrt(qx * qx + qy * qy + qz * qz);
+ if (!(distance <= w))
+ {
+ misses[family]++;
+ if (witness.Count < 8)
+ witness.Add($"{family}: iv_frustum_ball({slope}, {z0:R}, {z1:R}, {footprint:R}, {warp:R}) = {ball} misses slope {m}, depth {z:R} by {(distance - w) / w:R} r");
+ }
+ }
+ }
+ }
+
+ Assert.True(misses.Values.Sum() + pins.Values.Sum() == 0,
+ "misses " + string.Join(", ", families.Select(f => $"{f} {misses[f]}")) + "; pin failures " + string.Join(", ", families.Select(f => $"{f} {pins[f]}"))
+ + Environment.NewLine + string.Join(Environment.NewLine, witness));
+ }
+
+ ///
+ /// intervals-enclose for the composition the site and Aetheria use: over 2,000 seeded tile slices
+ /// (Aetheria's wells, inside the fog, the void at every grid point, the shipped void and LOD-far;
+ /// N in {1, 4, 8, 16}; a camera drawn per tile) the bound composed from the envelope over the slice's
+ /// box, its fade, and the noise ball enlarged by the warp encloses the density at 64 points of the
+ /// slice, each sampled through the field's own warp, with no tolerance. The void's footprint-truncated density (octave weights) is checked
+ /// against the same bound; a quarter of the slices are LOD-far's, where the weights are strictly
+ /// between 0 and 1 and at 0, and the test asserts it sampled both.
+ ///
+ [Fact]
+ public void EnvelopeBoundEnclosesDensity()
+ {
+ var random = new System.Random(0xE7B0);
+ int[] sizes = { 1, 4, 8, 16 };
+ var fields = new List { FogField.Wells(FogField.AetheriaWarp), FogField.InsideFog(FogField.AetheriaWarp) };
+ foreach (var rh in VoidField.HollowRadii)
+ foreach (var offset in VoidField.CameraOffsets)
+ foreach (var ramp in VoidField.Ramps)
+ fields.Add(new VoidField(rh, ramp, offset.Fraction, offset.Name));
+ fields.Add(VoidField.Shipped());
+ fields.Add(VoidField.LodFar());
+ var counts = new Counts();
+ var (partial, zero) = (0, 0);
+ for (var t = 0; t < 2000; t++)
+ {
+ var field = t % 2 == 0 ? fields[(t / 2) % 2] : t % 4 == 1 ? fields[^1] : fields[2 + random.Next(fields.Count - 2)];
+ var tile = field.DrawTile(random, sizes[random.Next(sizes.Length)]);
+ var near = field.Grid[0];
+ var far = field.Grid[^1];
+ var z0 = Uniform(random, near, far * 0.9f);
+ var z1 = z0 + LogUniform(random, 1.0e-3f * far, far - z0);
+ var bound = field.Bound(tile, z0, z1, gated: false, counts);
+ for (var i = 0; i < 64; i++)
+ {
+ var m = tile.PixelSlope(random.Next(tile.N), random.Next(tile.N), random.NextSingle(), random.NextSingle());
+ var z = i < 2 ? (i == 0 ? z0 : z1) : Uniform(random, z0, z1);
+ var density = field.Density(tile, m, z, lod: false, counts);
+ Assert.True(Inside(density, bound), $"{field.Name}: bound {bound} over [{z0:R}, {z1:R}] misses density {density:R} at slope {m}, depth {z:R}");
+ if (field is VoidField voidField)
+ {
+ var (wc, wd) = voidField.Weights(z);
+ partial += (wc > 0.0f && wc < 1.0f) || (wd > 0.0f && wd < 1.0f) ? 1 : 0;
+ zero += wd == 0.0f ? 1 : 0;
+ var truncated = field.Density(tile, m, z, lod: true, counts);
+ Assert.True(Inside(truncated, bound), $"{field.Name}: bound {bound} misses the truncated density {truncated:R} at slope {m}, depth {z:R}");
+ }
+ }
+ }
+
+ // The truncated check is not vacuous: some samples carry an octave weight strictly inside (0, 1)
+ // and some a fine weight of 0 (LOD-far reaches both).
+ Assert.True(partial > 0 && zero > 0, $"octave weights in (0, 1) at {partial} samples and fine weight 0 at {zero}");
+ }
+
+ ///
+ /// VoidField's octave weight bound over [z0, z1] encloses each octave's weight at every depth of the
+ /// segment, and, for a one-signed noise interval n = [0.25, 1], the weighted noise w(z) n_s lies in
+ /// iv_mul(bound, n). 2,000 LOD-far segments with z0 ~ U(1460, 0.95 grid end), z1 ~ U(z0, grid end),
+ /// 16 depths each including both ends. A bound whose lower end is w(z0) misses w(z1) 0.25 wherever the
+ /// weight falls across the segment.
+ ///
+ [Fact]
+ public void WeightBoundEnclosesOneSignedNoise()
+ {
+ var field = VoidField.LodFar();
+ var random = new System.Random(0x3B0D);
+ var end = field.Grid[^1];
+ var n = new float2(0.25f, 1.0f);
+ for (var i = 0; i < 2000; i++)
+ {
+ var z0 = Uniform(random, 1460.0f, 0.95f * end);
+ var z1 = Uniform(random, z0, end);
+ var (coarseBound, fineBound) = field.WeightBounds(z0, z1);
+ for (var s = 0; s < 16; s++)
+ {
+ var z = s == 0 ? z0 : s == 15 ? z1 : Uniform(random, z0, z1);
+ var (coarse, fine) = field.Weights(z);
+ foreach (var (octave, w, bound) in new[] { ("coarse", coarse, coarseBound), ("fine", fine, fineBound) })
+ {
+ Assert.True(Inside(w, bound), $"{octave} weight {w:R} at depth {z:R} lies outside its bound {bound} over [{z0:R}, {z1:R}]");
+ var weighted = iv_mul(bound, n);
+ foreach (var ns in new[] { 0.25f, Uniform(random, 0.25f, 1.0f), 1.0f })
+ Assert.True(Inside(w * ns, weighted), $"{octave} weight {w:R} x noise {ns:R} at depth {z:R} lies outside iv_mul({bound}, {n}) = {weighted} over [{z0:R}, {z1:R}]");
+ }
+ }
+ }
+ }
+
+ // ---- The Lipschitz constant ----
+
+ private static float GradientNorm(float3 p)
+ {
+ var g = snoise_grad(p);
+ return MathF.Sqrt(g.x * g.x + g.y * g.y + g.z * g.z);
+ }
+
+ // Hill-climbs |snoise_grad| from p: the ascent direction is the central difference of |grad|, the
+ // step grows on success and halves on failure. Returns the largest |grad| it reached.
+ private static float AscendGradientNorm(float3 p)
+ {
+ const float h = 1.0e-3f;
+ var step = 1.0e-2f;
+ var best = GradientNorm(p);
+ for (var k = 0; k < 200 && step > 1.0e-6f; k++)
+ {
+ var d = new float3(
+ GradientNorm(p + new float3(h, 0.0f, 0.0f)) - GradientNorm(p - new float3(h, 0.0f, 0.0f)),
+ GradientNorm(p + new float3(0.0f, h, 0.0f)) - GradientNorm(p - new float3(0.0f, h, 0.0f)),
+ GradientNorm(p + new float3(0.0f, 0.0f, h)) - GradientNorm(p - new float3(0.0f, 0.0f, h)));
+ if (dot(d, d) == 0.0f)
+ break;
+ var q = p + normalize(d) * step;
+ var value = GradientNorm(q);
+ if (value > best)
+ {
+ p = q;
+ best = value;
+ step *= 1.5f;
+ }
+ else
+ {
+ step *= 0.5f;
+ }
+ }
+
+ return best;
+ }
+
+ // Same climb on |snoise|, along the value's own analytic gradient.
+ private static float AscendAbsValue(float3 p)
+ {
+ var step = 1.0e-2f;
+ var best = MathF.Abs(snoise(p));
+ for (var k = 0; k < 200 && step > 1.0e-6f; k++)
+ {
+ var g = snoise_grad(p);
+ var d = new float3(g.x, g.y, g.z) * (g.w < 0.0f ? -1.0f : 1.0f);
+ if (dot(d, d) == 0.0f)
+ break;
+ var q = p + normalize(d) * step;
+ var value = MathF.Abs(snoise(q));
+ if (value > best)
+ {
+ p = q;
+ best = value;
+ step *= 1.5f;
+ }
+ else
+ {
+ step *= 0.5f;
+ }
+ }
+
+ return best;
+ }
+
+ // The largest |snoise_grad| over `samples` seeded points in [-256, 256]^3, the largest after refining
+ // the `refine` largest of them by ascent, and the start whose ascent reached it (the witness).
+ private static (float Sampled, float Refined, float3 Witness) MeasureGradient(int seed, int samples, int refine)
+ {
+ var random = new System.Random(seed);
+ var points = new float3[samples];
+ var norms = new float[samples];
+ for (var i = 0; i < samples; i++)
+ {
+ points[i] = new float3(Uniform(random, -256.0f, 256.0f), Uniform(random, -256.0f, 256.0f), Uniform(random, -256.0f, 256.0f));
+ norms[i] = GradientNorm(points[i]);
+ }
+
+ var order = Enumerable.Range(0, samples).OrderByDescending(i => norms[i]).Take(Math.Max(refine, 1)).ToArray();
+ var refined = norms[order[0]];
+ var witness = points[order[0]];
+ foreach (var i in order.Take(refine))
+ {
+ var value = AscendGradientNorm(points[i]);
+ if (value > refined)
+ {
+ refined = value;
+ witness = points[i];
+ }
+ }
+
+ return (norms[order[0]], refined, witness);
+ }
+
+ ///
+ /// The provenance of SNOISE_LIPSCHITZ: the largest |snoise_grad| over 1e6 seeded points, the 1e4
+ /// largest refined by gradient ascent, times 1.10. Prints the start whose ascent reached the maximum,
+ /// which LipschitzConstantPinsSampledGradients climbs from again. Also prints the largest |snoise|
+ /// reached the same way, which is what lets iv_snoise_ball intersect with [-1, 1]. Slow; run
+ /// explicitly after any change to the snoise kernel and re-pin the constant and the witness.
+ ///
+ [Fact(Explicit = true)]
+ [Trait("Category", "Slow")]
+ public void MeasureLipschitz()
+ {
+ var (sampled, refined, witness) = MeasureGradient(0x11B5, 1_000_000, 10_000);
+ var random = new System.Random(0xAB5);
+ var valueStarts = Enumerable.Range(0, 10_000)
+ .Select(_ => new float3(Uniform(random, -256.0f, 256.0f), Uniform(random, -256.0f, 256.0f), Uniform(random, -256.0f, 256.0f)));
+ var maxValue = valueStarts.Max(AscendAbsValue);
+ output.WriteLine($"IV-REPORT lipschitz: max |snoise_grad| sampled {sampled:R}, refined {refined:R}, x 1.10 = {refined * 1.10f:R}, witness ({witness.x:R}, {witness.y:R}, {witness.z:R}); max |snoise| refined {maxValue:R}");
+ Assert.True(refined <= SNOISE_LIPSCHITZ, $"measured {refined:R} exceeds SNOISE_LIPSCHITZ {SNOISE_LIPSCHITZ:R}");
+ Assert.True(maxValue < 1.0f, $"|snoise| reaches {maxValue:R}; iv_snoise_ball may not intersect with [-1, 1]");
+ }
+
+ // The start MeasureLipschitz printed as its witness.
+ private static readonly float3 LipschitzWitness = new(-241.98447f, -249.98927f, -60.058212f);
+
+ ///
+ /// Pins SNOISE_LIPSCHITZ: no |snoise_grad| over 1e5 seeded points exceeds it, and the ascent from
+ /// MeasureLipschitz's witness reaches a maximum that, times 1.10, is the constant to within half a
+ /// percent, so a constant lowered (or raised) by 1% fails here.
+ ///
+ [Fact]
+ public void LipschitzConstantPinsSampledGradients()
+ {
+ var (sampled, _, _) = MeasureGradient(0x11B6, 100_000, 0);
+ Assert.True(sampled <= SNOISE_LIPSCHITZ, $"|snoise_grad| reaches {sampled:R}, above SNOISE_LIPSCHITZ {SNOISE_LIPSCHITZ:R}");
+ var refined = AscendGradientNorm(LipschitzWitness);
+ Assert.True(refined <= SNOISE_LIPSCHITZ, $"|snoise_grad| reaches {refined:R}, above SNOISE_LIPSCHITZ {SNOISE_LIPSCHITZ:R}");
+ var ratio = SNOISE_LIPSCHITZ / (1.10f * refined);
+ Assert.True(MathF.Abs(ratio - 1.0f) <= 0.005f, $"SNOISE_LIPSCHITZ is {ratio:R} times 1.10 x the re-measured {refined:R}; re-run MeasureLipschitz and re-pin");
+ }
+
+ ///
+ /// Mean of interval width over the true range (max - min of the 64 sampled values) per ball: how
+ /// loose the Lipschitz bound is, for a later affine cut to beat. Asserts only finiteness.
+ ///
+ [Fact]
+ public void TightnessReport()
+ {
+ double Ratio(Func bound, Func f, float radius, int seed)
+ {
+ var random = new System.Random(seed);
+ var sum = 0.0;
+ for (var b = 0; b < 2000; b++)
+ {
+ var centre = RandomCentre(random);
+ var interval = bound(centre, radius);
+ float lo = float.MaxValue, hi = float.MinValue;
+ for (var i = 0; i < 64; i++)
+ {
+ var v = f(PointInBall(random, centre, radius, i));
+ lo = MathF.Min(lo, v);
+ hi = MathF.Max(hi, v);
+ }
+
+ sum += (interval.y - interval.x) / (double)(hi - lo);
+ }
+
+ return sum / 2000;
+ }
+
+ var snoiseRatios = new[] { 0.05f, 0.25f, 1.0f }.Select(r => (r, Ratio(iv_snoise_ball, snoise, r, 0x7167))).ToArray();
+ var fbmRatios = new[] { 0.05f, 0.25f }.Select(r => (r, Ratio((c, rr) => iv_fbm_ball(c, rr, 4, 2.0f, 0.5f), x => fbm_grad(x, 4, 2.0f, 0.5f).w, r, 0x7168))).ToArray();
+ output.WriteLine("IV-REPORT tightness (mean width / sampled range): " +
+ string.Join(", ", snoiseRatios.Select(t => $"snoise r={t.r:R} {t.Item2:F2}")) + "; " +
+ string.Join(", ", fbmRatios.Select(t => $"fbm4 r={t.r:R} {t.Item2:F2}")));
+ Assert.All(snoiseRatios.Concat(fbmRatios), t => Assert.True(double.IsFinite(t.Item2)));
+ }
+
+
+ // ---- The saving: envelope-then-noise fields, marched over screen tiles ----
+
+ // Combined cost in snoise units: one envelope evaluation, pointwise or as an interval probe, costs this
+ // much. An estimate (a handful of pulses or one height-map fetch against ~150 flops of snoise); both
+ // counts are printed, so any other weight can be applied to them.
+ internal const double EnvelopeCost = 0.2;
+
+ internal sealed class Counts
+ {
+ public long Snoise;
+ public long Envelope;
+
+ public double Cost => Snoise + EnvelopeCost * Envelope;
+ }
+
+ ///
+ /// One screen tile of one frame: a pinhole camera (origin and orthonormal basis; camera space looks
+ /// down +z and the ray of slope m is the points (m z, z)), the tile's central slope, its size N, the
+ /// focal length in pixels and the frame time. Every tile draws its own camera, so a probe result is
+ /// never carried to another frame.
+ ///
+ internal class Tile
+ {
+ public float3 Origin;
+ public float3 Right = new(1.0f, 0.0f, 0.0f);
+ public float3 Up = new(0.0f, 1.0f, 0.0f);
+ public float3 Forward = new(0.0f, 0.0f, 1.0f);
+ public float2 Slope;
+ public int N;
+ public float Focal;
+ public float Time;
+
+ public float Half => N / (2.0f * Focal);
+ public float Footprint => N / (MathF.Sqrt(2.0f) * Focal);
+ public float2 SlopesX => new(Slope.x - Half, Slope.x + Half);
+ public float2 SlopesY => new(Slope.y - Half, Slope.y + Half);
+
+ // Pixel (i, j) of the tile at sub-pixel offset (u, v) in [0, 1): its slope lies in SlopesX x SlopesY
+ // in float32, because (i + u) / f - N / (2f) is at most N / f - N / (2f) = N / (2f) exactly.
+ public float2 PixelSlope(int i, int j, float u, float v) =>
+ Slope + new float2((i + u) / Focal - Half, (j + v) / Focal - Half);
+
+ public float3 World(float3 c) => Origin + Right * c.x + Up * c.y + Forward * c.z;
+
+ public float3 Camera(float3 w)
+ {
+ var d = w - Origin;
+ return new float3(dot(d, Right), dot(d, Up), dot(d, Forward));
+ }
+
+ public void Look(float3 forward)
+ {
+ Forward = normalize(forward);
+ Right = normalize(cross(new float3(0.0f, 1.0f, 0.0f), Forward));
+ Up = cross(Forward, Right);
+ }
+ }
+
+ ///
+ /// A field with an analytic envelope and the dense grid it is marched on. Density is the pointwise
+ /// field, envelope-gated (noise only where the envelope leaves the density undecided), counted. Bound
+ /// is the one composition of the envelope's interval over a tile slice's box and the noise ball over the
+ /// slice's iv_frustum_ball: an interval enclosing Density at every point of every ray of the tile in the
+ /// slice. A slice is provably empty iff Bound(...).y <= 0. gated is the march's form: it returns the
+ /// envelope's verdict without snoise when the envelope alone proves the slice empty, and attempts the
+ /// noise ball only when the best outcome any centre value could give would prove it; ungated, the ball
+ /// is always evaluated.
+ ///
+ internal abstract class Scenario
+ {
+ public required string Name { get; init; }
+ public float[] Grid { get; protected init; } = Array.Empty();
+ public float Focal { get; init; } = 935.0f; // 1080 rows at a 60 degree vertical field of view
+ public float Extinction { get; init; } = 1.0f;
+
+ public abstract Tile DrawTile(System.Random random, int n);
+
+ public abstract float Density(Tile tile, float2 slope, float z, bool lod, Counts counts);
+
+ public abstract float2 Bound(Tile tile, float z0, float z1, bool gated, Counts counts);
+
+ // Cell boundaries near + (j / J)^2 (far - near): Aetheria's quadratic sample spacing
+ // (Assets/Shaders/Raymarching/CloudShader.shader:116).
+ protected static float[] Quadratic(float near, float far, int cells) =>
+ Enumerable.Range(0, cells + 1).Select(j => near + (float)j / cells * ((float)j / cells) * (far - near)).ToArray();
+
+ protected static float[] Cells(float step, float end) =>
+ Enumerable.Range(0, (int)MathF.Ceiling(end / step) + 1).Select(j => j * step).ToArray();
+ }
+
+ // iv_snoise_ball, with the snoise skipped where L r >= 2: there the bound is [-1, 1] whatever snoise(c)
+ // is (|snoise| < 1), so the result is the same and costs nothing.
+ private static float2 NoiseBall(float3 centre, float radius, Counts counts)
+ {
+ if (SNOISE_LIPSCHITZ * radius >= 2.0f)
+ return new float2(-1.0f, 1.0f);
+ counts.Snoise++;
+ return iv_snoise_ball(centre, radius);
+ }
+
+ // The best interval any centre value could give a ball of this radius, toward low or high values; the
+ // gate on whether a probe can prove anything. snoise(c) > -1, so a real ball's hi is at least -1 + L r,
+ // and likewise its lo is at most 1 - L r.
+ private static float2 BestBall(float radius, bool low)
+ {
+ var e = SNOISE_LIPSCHITZ * radius;
+ if (e >= 2.0f)
+ return new float2(-1.0f, 1.0f);
+ return low ? new float2(-1.0f, min(-1.0f + e, 1.0f)) : new float2(max(1.0f - e, -1.0f), 1.0f);
+ }
+
+ // Least and greatest distance along one axis from a box to a point, exact in float32: a coordinate in
+ // [lo, hi] is at least `near` and at most `far` from c, and squares, sums and sqrt are monotone.
+ private static (float Near, float Far) AxisDistance(float2 range, float c) =>
+ (c < range.x ? range.x - c : c > range.y ? c - range.y : 0.0f, max(abs(range.x - c), abs(range.y - c)));
+
+ ///
+ /// The noise of the fog and the void (r1's field shape): coverage is two snoise of the warped point,
+ /// cross-faded over two half-period phases; n = coverage + detail / 2, the detail one snoise at four
+ /// times the frequency. The flow is a unit-bounded 2D field times FlowAmplitude, so a phase moves a point
+ /// by at most Warp = FlowAmplitude * Period / 2. The flow is not counted: the site reads it from a
+ /// texture and Aetheria from its flow map.
+ ///
+ internal sealed class WarpedNoise
+ {
+ public required float F0 { get; init; }
+ public required float Warp { get; init; }
+ public float Period { get; init; } = 8.0f;
+ public float FlowFrequency { get; init; } = 1.0f / 512.0f;
+ public bool FlowInXY { get; init; }
+
+ private float FlowAmplitude => Warp * 2.0f / Period;
+
+ private float3 Flow(float3 p)
+ {
+ var q = (FlowInXY ? new float2(p.x, p.y) : new float2(p.x, p.z)) * FlowFrequency;
+ var v = new float2(snoise(q + new float2(17.0f, 3.0f)), snoise(q + new float2(-5.0f, 41.0f)));
+ v = v / max(1.0f, length(v)) * FlowAmplitude;
+ return FlowInXY ? new float3(v.x, v.y, 0.0f) : new float3(v.x, 0.0f, v.y);
+ }
+
+ private (float W0, float Shift0, float Shift1) Phases(float time)
+ {
+ var phase0 = frac(time / Period);
+ var phase1 = frac(time / Period + 0.5f);
+ return (1.0f - abs(2.0f * phase0 - 1.0f), (phase0 - 0.5f) * Period, (phase1 - 0.5f) * Period);
+ }
+
+ public float Coverage(float3 p, float time, Counts counts)
+ {
+ var (w0, shift0, shift1) = Phases(time);
+ var flow = Flow(p);
+ counts.Snoise += 2;
+ return w0 * snoise((p + flow * shift0) * F0) + (1.0f - w0) * snoise((p + flow * shift1) * F0);
+ }
+
+ public float Detail(float3 p, Counts counts)
+ {
+ counts.Snoise++;
+ return snoise(p * (4.0f * F0) + 17.0f);
+ }
+
+ // Both phases' warped points lie in a ball that carries the warp, so one ball encloses both snoise;
+ // the cross-fade is enclosed term by term with the frame's own weights.
+ public float2 CoverageBound(float2 ball, float time)
+ {
+ var w0 = Phases(time).W0;
+ return iv_add(iv_scale(ball, w0), iv_scale(ball, 1.0f - w0));
+ }
+
+ public float2 CoverageBall(float3 centre, float radius, Counts counts) => NoiseBall(centre * F0, radius * F0, counts);
+
+ public float2 DetailBall(float3 centre, float radius, Counts counts) => NoiseBall(centre * (4.0f * F0) + 17.0f, radius * (4.0f * F0), counts);
+
+ public float CoverageRadius(float radius) => radius * F0;
+
+ public float DetailRadius(float radius) => radius * (4.0f * F0);
+ }
+
+ ///
+ /// Aetheria's nebula (Assets/Shaders/Volumetric.cginc:165-221, cloudDensity :91-102) over its zone
+ /// height map (Assets/Scripts/ServerShared/Zone.cs:369-397 GetHeight, :426 PowerPulse), in its units
+ /// (Assets/Resources/Settings.asset: ZoneDepth :97-98, DefaultEnvironment :386-411). The pre-distortion
+ /// SDF is s(p) = p.y + h(p.xz), h >= 0 the well map: the zone bowl PowerPulse(|xz| / 2R, 2) * 64 plus,
+ /// per well, PowerPulse(|xz - c_b| / r_b, 16) * d_b, with r_b = 500 M^0.25, d_b = 30 M^0.175 and
+ /// PowerPulse(x, e) = (1 - 4x^2)^e on [0, 1/2], 0 beyond. Noise displaces only below the safety distance:
+ /// fade = 1 - smoothstep(0.75 S, S, s), s' = s + A fade n(warp(p)), density = max(0, (F - s') / B); the
+ /// patch band and the 1e-9 fill are left out. The warp moves only the noise argument; the height is read
+ /// at the unwarped xz, as Aetheria reads _NebulaSurfaceHeight at pos.xz. Envelope coordinates are the
+ /// tile's camera frame (the camera yaws about y, so horizontal distances are the world's). Extinction
+ /// 0.5 is the environment's (Settings.asset:397, applied at CloudShader.shader:91).
+ ///
+ internal sealed class FogField : Scenario
+ {
+ public const float Safety = 30.0f; // SafetyDistance, Settings.asset:399
+ public const float FloorOffset = -20.0f; // :393
+ public const float FloorBlend = 10.0f; // :394
+ public const float Amplitude = 20.0f;
+ public const float NoiseScale = 414.2167f; // :409
+ public const float AetheriaWarp = 60.0f; // Flow GlobalAmplitude 15 x Period 8 / 2, :402-405
+ public const float ZoneRadius = 2000.0f;
+ private static readonly float2 FreeNoise = new(-1.5f, 1.5f);
+
+ public required WarpedNoise Noise { get; init; }
+ public (float2 Centre, float Scale, float Exponent, float Depth)[] Bowls { get; init; } = Array.Empty<(float2, float, float, float)>();
+ public required Func CameraAt { get; init; }
+ public required float2 SlopesY { get; init; }
+
+ public FogField()
+ {
+ // Near 0.3 and far 2048: the ARPG scene's Main Camera, which carries the cloud march
+ // (Aetheria d3075eb0 Assets/Scenes/ARPG.unity:31777 m_Name, :31830 far clip plane;
+ // Assets/Scripts/Zone Display/VolumeCloudRenderer.cs:36 RequireComponent(Camera)).
+ Grid = Quadratic(0.3f, 2048.0f, 256);
+ Extinction = 0.5f;
+ }
+
+ private static WarpedNoise AetheriaNoise(float warp) => new() { F0 = 1.0f / NoiseScale, Warp = warp };
+
+ /// (a) Height fog: no wells, the camera above the safety band at y = S + 30, rays level and up.
+ public static FogField HeightFog(float warp) => new()
+ {
+ Name = "(a) height fog", Noise = AetheriaNoise(warp), SlopesY = new float2(0.0f, 0.5f),
+ CameraAt = (random, _) => new float3(Uniform(random, -1000.0f, 1000.0f), Safety + 30.0f, Uniform(random, -1000.0f, 1000.0f)),
+ };
+
+ /// (b) Inside the fog: no wells, the camera at y = F - 10, rays level.
+ public static FogField InsideFog(float warp) => new()
+ {
+ Name = "(b) inside the fog", Noise = AetheriaNoise(warp), SlopesY = new float2(-0.05f, 0.05f),
+ CameraAt = (random, _) => new float3(Uniform(random, -1000.0f, 1000.0f), FloorOffset - 10.0f, Uniform(random, -1000.0f, 1000.0f)),
+ };
+
+ ///
+ /// (c) Aetheria-like: the zone bowl R = 2000 and four wells of mass 100, 1000, 10000 and 1000 at seeded
+ /// positions inside it; the camera at a seeded xz in the bowl, above the noise band at
+ /// y = max(0, S + 10 - h), gazing across with m_y in [-0.25, 0.1].
+ ///
+ public static FogField Wells(float warp)
+ {
+ var random = new System.Random(0xA37E);
+ var bowls = new List<(float2, float, float, float)> { (new float2(0.0f, 0.0f), 2.0f * ZoneRadius, 2.0f, 64.0f) };
+ foreach (var mass in new[] { 100.0f, 1000.0f, 10000.0f, 1000.0f })
+ {
+ var angle = Uniform(random, 0.0f, 2.0f * MathF.PI);
+ var radius = ZoneRadius * MathF.Sqrt(random.NextSingle());
+ bowls.Add((new float2(radius * MathF.Cos(angle), radius * MathF.Sin(angle)), 500.0f * MathF.Pow(mass, 0.25f), 16.0f, 30.0f * MathF.Pow(mass, 0.175f)));
+ }
+
+ return new FogField
+ {
+ Name = "(c) Aetheria wells", Noise = AetheriaNoise(warp), Bowls = bowls.ToArray(), SlopesY = new float2(-0.25f, 0.1f),
+ CameraAt = (r, field) =>
+ {
+ var angle = Uniform(r, 0.0f, 2.0f * MathF.PI);
+ var radius = ZoneRadius * MathF.Sqrt(r.NextSingle());
+ var xz = new float2(radius * MathF.Cos(angle), radius * MathF.Sin(angle));
+ var h = field.Height(field.Bowls.Select(b => b.Centre - xz).ToArray(), 0.0f, 0.0f);
+ return new float3(xz.x, max(0.0f, Safety + 10.0f - h), xz.y);
+ },
+ };
+ }
+
+ private sealed class FogTile : Tile
+ {
+ public float2[] Centres = Array.Empty(); // bowl centres in the camera's (x, z)
+ }
+
+ public override Tile DrawTile(System.Random random, int n)
+ {
+ var origin = CameraAt(random, this);
+ var yaw = Uniform(random, 0.0f, 2.0f * MathF.PI);
+ var tile = new FogTile
+ {
+ Origin = origin, N = n, Focal = Focal, Time = Uniform(random, 0.0f, 100.0f),
+ Slope = new float2(Uniform(random, -1.0f, 1.0f), Uniform(random, SlopesY.x, SlopesY.y)),
+ };
+ tile.Look(new float3(MathF.Sin(yaw), 0.0f, MathF.Cos(yaw)));
+ tile.Centres = Bowls.Select(b =>
+ {
+ var c = tile.Camera(new float3(b.Centre.x, origin.y, b.Centre.y));
+ return new float2(c.x, c.z);
+ }).ToArray();
+ return tile;
+ }
+
+ private static float Pulse(float distance, float scale, float exponent)
+ {
+ var x = distance / scale;
+ var t = 1.0f - 4.0f * x * x;
+ return t <= 0.0f ? 0.0f : pow(t, exponent);
+ }
+
+ // h at the camera-frame point (x, z), bowl centres in the same frame.
+ private float Height(float2[] centres, float x, float z)
+ {
+ var h = 0.0f;
+ for (var b = 0; b < Bowls.Length; b++)
+ {
+ var dx = x - centres[b].x;
+ var dz = z - centres[b].y;
+ h += Pulse(sqrt(dx * dx + dz * dz), Bowls[b].Scale, Bowls[b].Exponent) * Bowls[b].Depth;
+ }
+
+ return h;
+ }
+
+ // Each bowl's term over the box is [pulse(far), pulse(near)], pulse decreasing in distance; pow is
+ // not required to be correctly rounded, so both ends are widened by one ulp.
+ private float2 HeightBound(float2[] centres, float2 x, float2 z)
+ {
+ var h = iv_point(0.0f);
+ for (var b = 0; b < Bowls.Length; b++)
+ {
+ var (nx, fx) = AxisDistance(x, centres[b].x);
+ var (nz, fz) = AxisDistance(z, centres[b].y);
+ var near = Pulse(sqrt(nx * nx + nz * nz), Bowls[b].Scale, Bowls[b].Exponent);
+ var far = Pulse(sqrt(fx * fx + fz * fz), Bowls[b].Scale, Bowls[b].Exponent);
+ var pulse = new float2(far > 0.0f ? MathF.BitDecrement(far) : 0.0f, near > 0.0f ? MathF.BitIncrement(near) : 0.0f);
+ h = iv_add(h, iv_scale(pulse, Bowls[b].Depth));
+ }
+
+ return h;
+ }
+
+ private static float2 Compose(float2 s, float2 fade, float2 n)
+ {
+ var displaced = iv_add(s, iv_scale(iv_mul(fade, n), Amplitude));
+ return new float2(max(0.0f, (FloorOffset - displaced.y) / FloorBlend), max(0.0f, (FloorOffset - displaced.x) / FloorBlend));
+ }
+
+ public override float Density(Tile t, float2 slope, float z, bool lod, Counts counts)
+ {
+ var tile = (FogTile)t;
+ var x = slope.x * z;
+ var y = slope.y * z;
+ counts.Envelope++;
+ var s = tile.Origin.y + y + Height(tile.Centres, x, z);
+ var fade = 1.0f - smoothstep(0.75f * Safety, Safety, s);
+ if (Compose(iv_point(s), iv_point(fade), FreeNoise).y <= 0.0f)
+ return 0.0f;
+ var p = tile.World(new float3(x, y, z));
+ var coverage = Noise.Coverage(p, tile.Time, counts);
+ if (Compose(iv_point(s), iv_point(fade), iv_add(iv_point(coverage), iv_scale(new float2(-1.0f, 1.0f), 0.5f))).y <= 0.0f)
+ return 0.0f;
+ var n = coverage + Noise.Detail(p, counts) * 0.5f;
+ return max(0.0f, (FloorOffset - (s + fade * n * Amplitude)) / FloorBlend);
+ }
+
+ public override float2 Bound(Tile t, float z0, float z1, bool gated, Counts counts)
+ {
+ var tile = (FogTile)t;
+ var depth = new float2(z0, z1);
+ counts.Envelope++;
+ var height = HeightBound(tile.Centres, iv_mul(depth, tile.SlopesX), depth);
+ var s = iv_add(iv_add(iv_point(tile.Origin.y), iv_mul(depth, tile.SlopesY)), height);
+ var fade = iv_sub(iv_point(1.0f), iv_smoothstep(0.75f * Safety, Safety, s));
+ var envelope = Compose(s, fade, FreeNoise);
+ if (gated && envelope.y <= 0.0f)
+ return envelope;
+ var ball = iv_frustum_ball(tile.Slope, z0, z1, tile.Footprint, Noise.Warp);
+ if (gated)
+ {
+ var provable = false;
+ foreach (var low in new[] { true, false })
+ {
+ var best = iv_add(Noise.CoverageBound(BestBall(Noise.CoverageRadius(ball.w), low), tile.Time), iv_scale(BestBall(Noise.DetailRadius(ball.w), low), 0.5f));
+ provable |= Compose(s, fade, best).y <= 0.0f;
+ }
+
+ if (!provable)
+ return envelope;
+ }
+
+ var centre = tile.World(new float3(ball.x, ball.y, ball.z));
+ var noise = iv_add(Noise.CoverageBound(Noise.CoverageBall(centre, ball.w, counts), tile.Time), iv_scale(Noise.DetailBall(centre, ball.w, counts), 0.5f));
+ return Compose(s, fade, noise);
+ }
+ }
+
+ ///
+ /// (d) r1's uniform slab, kept as the worst case: the map's warped coverage field (F0 0.35, cutoff 0.25,
+ /// depth 12, period 4, extinction 2, flow in xy at 0.15, warp 0.05) with no envelope, so every cell is a
+ /// noise question. The dense sample evaluates the coverage first and the detail only where the coverage
+ /// passes the cutoff; ungated it would cost three snoise, as r1 counted.
+ ///
+ internal sealed class SlabField : Scenario
+ {
+ public const float F0 = 0.35f;
+ public const float Cutoff = 0.25f;
+ public const float Depth = 12.0f;
+ public const float SlabExtinction = 2.0f;
+ public const float SlabWarp = 0.05f;
+
+ public required WarpedNoise Noise { get; init; }
+
+ public SlabField() => Grid = Cells(0.125f, Depth);
+
+ public static SlabField Create(float warp) => new()
+ {
+ Name = "(d) uniform slab", Noise = new WarpedNoise { F0 = F0, Warp = warp, Period = 4.0f, FlowFrequency = 0.15f, FlowInXY = true },
+ };
+
+ public override Tile DrawTile(System.Random random, int n) => new()
+ {
+ Origin = new float3(Uniform(random, -50.0f, 50.0f), Uniform(random, -50.0f, 50.0f), 0.0f),
+ Slope = new float2(Uniform(random, -0.5f, 0.5f), Uniform(random, -0.5f, 0.5f)),
+ N = n, Focal = Focal, Time = Uniform(random, 0.0f, 100.0f),
+ };
+
+ public override float Density(Tile tile, float2 slope, float z, bool lod, Counts counts)
+ {
+ var p = tile.World(new float3(slope.x * z, slope.y * z, z));
+ var coverage = Noise.Coverage(p, tile.Time, counts);
+ if (coverage <= Cutoff)
+ return 0.0f;
+ counts.Snoise++;
+ var detail = 1.0f - abs(snoise(p * (4.0f * F0) + 17.0f));
+ return max(0.0f, coverage - Cutoff) * SlabExtinction * lerp(0.6f, 1.0f, detail);
+ }
+
+ public override float2 Bound(Tile tile, float z0, float z1, bool gated, Counts counts)
+ {
+ var maybe = new float2(0.0f, float.MaxValue);
+ var ball = iv_frustum_ball(tile.Slope, z0, z1, tile.Footprint, Noise.Warp);
+ if (gated && Noise.CoverageBound(BestBall(Noise.CoverageRadius(ball.w), low: true), tile.Time).y > Cutoff)
+ return maybe;
+ var coverage = Noise.CoverageBound(Noise.CoverageBall(tile.World(new float3(ball.x, ball.y, ball.z)), ball.w, counts), tile.Time);
+ return coverage.y <= Cutoff ? new float2(0.0f, 0.0f) : maybe;
+ }
+ }
+
+ ///
+ /// (e) The void: a cloud volume everywhere, carved by one negative pseudo-Gaussian sphere. Authored
+ /// knobs: hollow radius Rh (density exactly 0 inside) and ramp width S (the density rises to full over
+ /// [Rh, Rh + S]); everything else is derived from them. The brush (1 - (d / Rc)^2)^e with CARVE = 1.5
+ /// reaches full at Rc, so Rc = Rh + S, and Rh = Rc sqrt(1 - CARVE^(-1/e)) gives
+ /// e = -ln CARVE / ln(1 - (Rh / Rc)^2): density = K max(0, 1 - CARVE (1 - (d' / Rc)^2)^e), exactly 0 for
+ /// d' <= Rh and K for d' >= Rc, K = 1/30. The noise displaces the wall: fade =
+ /// smoothstep(Rh - S, Rh, d), d' = d + A fade n(warp(p)), A = 20, F0 = 0.01 (wavelengths 100 and 25),
+ /// warp D = 10. The sun orbits the origin at radius 175 with period 72 s and rests on the hollow's
+ /// bottom (c_v = sun + (0, Rh - 12.5, 0)); the camera looks at sun + (0, 60, 0) and, shipped, sits 120
+ /// behind the sun along the orbit tangent and 45 above it.
+ ///
+ /// The one shipped void (ruling operator-shipped-void) is Shipped(): (Rh, S) = (198, 50), so Rc = 248
+ /// and e = -ln 1.5 / ln(1 - (198 / 248)^2), which the constructor derives; no other void is shipped.
+ ///
+ /// The footprint-truncated field (lod) weights octave i by w_i(z) = 1 - smoothstep(l_i / 2m, l_i / m,
+ /// z theta), m = 4, theta = 1 / f. Deviation from the spec's [w_i(z1), w_i(z0)]: the bound takes each
+ /// weight over [w_i(z1), 1], because one Bound serves the full field (w = 1) and the truncated one
+ /// (w in [w_i(z1), w_i(z0)], w decreasing in z); [w_i(z1), 1] contains both and iv_mul is monotone in
+ /// its interval argument, so the composed interval encloses both fields. The slack is only in far
+ /// tiles, whose step the footprint already sets. At f = 935 every weight is 1 inside every grid;
+ /// LodFar() (f = 467) reaches weights in (0, 1) and at 0.
+ ///
+ /// The dense grid is the fixed-step reference's: cells of min(10, S / 4) out to
+ /// |camera - c_v| + Rc + 6 / K. The density is provably K only from d >= Rc + A 1.5 (the free noise
+ /// interval), so 6 / K leaves at least 4.9 of optical depth past it and every ray goes opaque before
+ /// the grid ends.
+ ///
+ internal sealed class VoidField : Scenario
+ {
+ public const float Carve = 1.5f;
+ public const float Amplitude = 20.0f;
+ public const float K = 1.0f / 30.0f;
+ public const float StepFootprints = 4.0f; // k
+ public const float OctaveFootprints = 4.0f; // m
+ private static readonly float2 FreeNoise = new(-1.5f, 1.5f);
+
+ public static readonly float[] HollowRadii = { 100.0f, 200.0f, 400.0f, 1000.0f };
+ public static readonly float[] Ramps = { 10.0f, 50.0f, 150.0f };
+ public static readonly (string Name, float Fraction)[] CameraOffsets = { ("centred", 0.0f), ("mid", 0.5f), ("low", 0.94f) };
+
+ public readonly float HollowRadius;
+ public readonly float Ramp;
+ public readonly float CarveRadius;
+ public readonly float Exponent;
+ public readonly float RampStart;
+ public readonly float FixedStep;
+ public readonly float CameraFraction; // negative: the shipped camera
+ public readonly WarpedNoise Noise = new() { F0 = 0.01f, Warp = 10.0f };
+
+ private float CoarseWavelength => 1.0f / Noise.F0;
+ private float FineWavelength => 1.0f / (4.0f * Noise.F0);
+ private float Theta => 1.0f / Focal;
+
+ [System.Diagnostics.CodeAnalysis.SetsRequiredMembers]
+ public VoidField(float hollowRadius, float ramp, float cameraFraction, string camera, float focal = 935.0f)
+ {
+ Focal = focal;
+ HollowRadius = hollowRadius;
+ Ramp = ramp;
+ CarveRadius = hollowRadius + ramp;
+ Exponent = (float)(-Math.Log(Carve) / Math.Log(1.0 - (double)hollowRadius / CarveRadius * hollowRadius / CarveRadius));
+ RampStart = hollowRadius - ramp;
+ FixedStep = min(10.0f, ramp / 4.0f);
+ CameraFraction = cameraFraction;
+ Name = $"(e) void Rh={hollowRadius:R} camera={camera} S={ramp:R}" + (focal == 935.0f ? "" : $" f={focal:R}");
+ var reach = cameraFraction < 0.0f ? length(new float2(120.0f, 45.0f - (hollowRadius - 12.5f))) : cameraFraction * hollowRadius;
+ Grid = Cells(FixedStep, reach + CarveRadius + 6.0f / K);
+ }
+
+ public static VoidField Shipped() => new(198.0f, 50.0f, -1.0f, "shipped");
+
+ ///
+ /// The configuration that reaches the octave weights: Rh = 4000, S = 50, the low camera, f = 467
+ /// (1080 rows at a 98 degree vertical field of view). The fine weight is below 1 past z = 1460 and 0
+ /// past 2920, the coarse below 1 past 5840, and the grid reaches about 7990. Printed, not asserted.
+ ///
+ public static VoidField LodFar() => new(4000.0f, 50.0f, 0.94f, "low", 467.0f);
+
+ /// The coarse and fine octave weights at depth z (1 for the full field).
+ public (float Coarse, float Fine) Weights(float z) => (Weight(CoarseWavelength, z * Theta), Weight(FineWavelength, z * Theta));
+
+ /// The coarse and fine octave weight bounds Bound uses over [z0, z1].
+ public (float2 Coarse, float2 Fine) WeightBounds(float z0, float z1) => (WeightBound(CoarseWavelength, z0, z1), WeightBound(FineWavelength, z0, z1));
+
+ private sealed class VoidTile : Tile
+ {
+ public float3 Cavity; // c_v in the camera frame
+ }
+
+ public override Tile DrawTile(System.Random random, int n)
+ {
+ var time = Uniform(random, 0.0f, 72.0f);
+ var angle = 2.0f * MathF.PI * time / 72.0f;
+ var sun = new float3(175.0f * MathF.Cos(angle), 0.0f, 175.0f * MathF.Sin(angle));
+ var tangent = new float3(-MathF.Sin(angle), 0.0f, MathF.Cos(angle));
+ var cavity = sun + new float3(0.0f, HollowRadius - 12.5f, 0.0f);
+ var target = sun + new float3(0.0f, 60.0f, 0.0f);
+ float3 origin;
+ float3 forward;
+ if (CameraFraction < 0.0f)
+ {
+ origin = sun - tangent * 120.0f + new float3(0.0f, 45.0f, 0.0f);
+ forward = target - origin;
+ }
+ else
+ {
+ // The shipped camera's direction from the centre, at a fraction of this hollow's radius;
+ // centred, it looks along the orbit, slightly down.
+ origin = cavity + normalize(tangent * -120.0f + new float3(0.0f, 45.0f - (198.0f - 12.5f), 0.0f)) * (CameraFraction * HollowRadius);
+ forward = CameraFraction == 0.0f ? tangent * -1.0f + new float3(0.0f, -0.2f, 0.0f) : target - origin;
+ }
+
+ var tile = new VoidTile
+ {
+ Origin = origin, N = n, Focal = Focal, Time = time,
+ Slope = new float2((random.Next(1920 / n) * n + n * 0.5f - 960.0f) / Focal, (random.Next(1080 / n) * n + n * 0.5f - 540.0f) / Focal),
+ };
+ tile.Look(forward);
+ tile.Cavity = tile.Camera(cavity);
+ return tile;
+ }
+
+ private static float Weight(float wavelength, float footprint) =>
+ 1.0f - smoothstep(wavelength / (2.0f * OctaveFootprints), wavelength / OctaveFootprints, footprint);
+
+ private float2 WeightBound(float wavelength, float z0, float z1) =>
+ new(iv_sub(iv_point(1.0f), iv_smoothstep(wavelength / (2.0f * OctaveFootprints), wavelength / OctaveFootprints, new float2(z0 * Theta, z1 * Theta))).x, 1.0f);
+
+ /// The footprint step: max(k z theta, min(S / 4, half the finest live wavelength)).
+ public float LodStep(float z)
+ {
+ var footprint = z * Theta;
+ var fine = Weight(FineWavelength, footprint) > 0.0f ? FineWavelength : Weight(CoarseWavelength, footprint) > 0.0f ? CoarseWavelength : float.PositiveInfinity;
+ return max(StepFootprints * footprint, min(Ramp / 4.0f, fine / 2.0f));
+ }
+
+ private float Profile(float displaced)
+ {
+ if (displaced <= HollowRadius)
+ return 0.0f;
+ var x = displaced / CarveRadius;
+ var g = displaced < CarveRadius ? pow(max(0.0f, 1.0f - x * x), Exponent) : 0.0f;
+ return K * max(0.0f, 1.0f - Carve * g);
+ }
+
+ private float2 Compose(float2 d, float2 fade, float2 n)
+ {
+ var displaced = iv_add(d, iv_scale(iv_mul(fade, n), Amplitude));
+ return new float2(displaced.x >= CarveRadius ? K : 0.0f, displaced.y <= HollowRadius ? 0.0f : K);
+ }
+
+ public override float Density(Tile t, float2 slope, float z, bool lod, Counts counts)
+ {
+ var tile = (VoidTile)t;
+ var x = slope.x * z;
+ var y = slope.y * z;
+ counts.Envelope++;
+ var dx = x - tile.Cavity.x;
+ var dy = y - tile.Cavity.y;
+ var dz = z - tile.Cavity.z;
+ var d = sqrt(dx * dx + dy * dy + dz * dz);
+ var fade = smoothstep(RampStart, HollowRadius, d);
+ var free = Compose(iv_point(d), iv_point(fade), FreeNoise);
+ if (free.x == free.y)
+ return free.x;
+ var p = tile.World(new float3(x, y, z));
+ var wc = lod ? Weight(CoarseWavelength, z * Theta) : 1.0f;
+ var wd = lod ? Weight(FineWavelength, z * Theta) : 1.0f;
+ var coverage = wc > 0.0f ? Noise.Coverage(p, tile.Time, counts) * wc : 0.0f;
+ var partial = Compose(iv_point(d), iv_point(fade), iv_add(iv_point(coverage), iv_scale(new float2(-wd, wd), 0.5f)));
+ if (partial.x == partial.y)
+ return partial.x;
+ var detail = wd > 0.0f ? Noise.Detail(p, counts) * 0.5f * wd : 0.0f;
+ return Profile(d + fade * (coverage + detail) * Amplitude);
+ }
+
+ ///
+ /// The field with no envelope gate: both noise terms always evaluated, the profile applied to the
+ /// displaced distance. Density's shortcuts go through Compose; this does not, so it is the value a
+ /// wrong Compose cannot also corrupt. Not counted.
+ ///
+ public float Exact(Tile t, float2 slope, float z, bool lod)
+ {
+ var tile = (VoidTile)t;
+ var counts = new Counts();
+ var (x, y) = (slope.x * z, slope.y * z);
+ var (dx, dy, dz) = (x - tile.Cavity.x, y - tile.Cavity.y, z - tile.Cavity.z);
+ var d = sqrt(dx * dx + dy * dy + dz * dz);
+ var fade = smoothstep(RampStart, HollowRadius, d);
+ var p = tile.World(new float3(x, y, z));
+ var wc = lod ? Weight(CoarseWavelength, z * Theta) : 1.0f;
+ var wd = lod ? Weight(FineWavelength, z * Theta) : 1.0f;
+ var n = Noise.Coverage(p, tile.Time, counts) * wc + Noise.Detail(p, counts) * 0.5f * wd;
+ return Profile(d + fade * n * Amplitude);
+ }
+
+ public override float2 Bound(Tile t, float z0, float z1, bool gated, Counts counts)
+ {
+ var tile = (VoidTile)t;
+ var depth = new float2(z0, z1);
+ counts.Envelope++;
+ var (nx, fx) = AxisDistance(iv_mul(depth, tile.SlopesX), tile.Cavity.x);
+ var (ny, fy) = AxisDistance(iv_mul(depth, tile.SlopesY), tile.Cavity.y);
+ var (nz, fz) = AxisDistance(depth, tile.Cavity.z);
+ var d = new float2(sqrt(nx * nx + ny * ny + nz * nz), sqrt(fx * fx + fy * fy + fz * fz));
+ var fade = iv_smoothstep(RampStart, HollowRadius, d);
+ var envelope = Compose(d, fade, FreeNoise);
+ if (gated && envelope.y <= 0.0f)
+ return envelope;
+ var ball = iv_frustum_ball(tile.Slope, z0, z1, tile.Footprint, Noise.Warp);
+ var wc = WeightBound(CoarseWavelength, z0, z1);
+ var wd = WeightBound(FineWavelength, z0, z1);
+ if (gated)
+ {
+ var provable = false;
+ foreach (var low in new[] { true, false })
+ {
+ var best = iv_add(iv_mul(Noise.CoverageBound(BestBall(Noise.CoverageRadius(ball.w), low), tile.Time), wc), iv_mul(iv_scale(BestBall(Noise.DetailRadius(ball.w), low), 0.5f), wd));
+ provable |= Compose(d, fade, best).y <= 0.0f;
+ }
+
+ if (!provable)
+ return envelope;
+ }
+
+ var centre = tile.World(new float3(ball.x, ball.y, ball.z));
+ var noise = iv_add(iv_mul(Noise.CoverageBound(Noise.CoverageBall(centre, ball.w, counts), tile.Time), wc), iv_mul(iv_scale(Noise.DetailBall(centre, ball.w, counts), 0.5f), wd));
+ return Compose(d, fade, noise);
+ }
+ }
+
+ // One cell's transmittance: the density at the cell's mid-depth, over the ray's length through it.
+ private static float Integrate(Scenario field, Tile tile, float2 slope, float z0, float z1, bool lod, Counts counts) =>
+ exp(-(field.Density(tile, slope, (z0 + z1) * 0.5f, lod, counts) * field.Extinction) * ((z1 - z0) * sqrt(slope.x * slope.x + slope.y * slope.y + 1.0f)));
+
+ ///
+ /// The tile pre-pass on cell-index ranges: probe the slice of the range; provably empty skips it and
+ /// doubles the next range; otherwise a range of more than one cell halves, and a single cell is dense.
+ /// It starts with the whole grid, runs to the grid's end (a pre-pass knows no ray's transmittance) and
+ /// yields the dense cells in order. Its probes are counted once per tile.
+ ///
+ /// The body (ruling operator-analytic-wall-body): when a single dense cell's own bound is a positive
+ /// point (the cell is provably constant), one more probe over [grid[j], grid end]; if that is a
+ /// positive point too, every point of every ray of the tile over the rest of the grid has exactly that
+ /// density, so BodyStart = j, BodyDensity = its value, and every remaining cell is appended dense
+ /// without probing (adding cells is always sound, so the mask keeps its meaning). BodyStart is -1 when
+ /// no body is proven.
+ ///
+ private static (List Dense, int BodyStart, float BodyDensity) MarchTile(Scenario field, Tile tile, Counts probes)
+ {
+ var grid = field.Grid;
+ var cells = grid.Length - 1;
+ var dense = new List();
+ var j = 0;
+ var span = cells;
+ while (j < cells)
+ {
+ var end = Math.Min(j + span, cells);
+ var bound = field.Bound(tile, grid[j], grid[end], gated: true, probes);
+ if (bound.y <= 0.0f)
+ {
+ j = end;
+ span = Math.Min(span * 2, cells);
+ }
+ else if (end - j > 1)
+ {
+ span = (end - j + 1) / 2;
+ }
+ else
+ {
+ if (bound.x == bound.y && bound.x > 0.0f)
+ {
+ var rest = field.Bound(tile, grid[j], grid[cells], gated: true, probes);
+ if (rest.x == rest.y && rest.x > 0.0f)
+ {
+ dense.AddRange(Enumerable.Range(j, cells - j));
+ return (dense, j, rest.x);
+ }
+ }
+
+ dense.Add(j);
+ j++;
+ }
+ }
+
+ return (dense, -1, 0.0f);
+ }
+
+ // The body's optical depth from z to the grid's end, at the density the field's pointwise envelope gives
+ // at z (one envelope evaluation; no probe in the main pass). MarchLod takes exp(-BodyDepth) as one step.
+ private static float BodyDepth(Scenario field, Tile tile, float2 slope, float z, bool lod, Counts counts) =>
+ field.Density(tile, slope, z, lod, counts) * field.Extinction * ((field.Grid[^1] - z) * sqrt(slope.x * slope.x + slope.y * slope.y + 1.0f));
+
+ // The footprint-aware march, the void's shipped one: steps of LodStep(z) inside the unmasked runs only,
+ // octave weights applied, early-out below cutoff; at the body's start, one analytic step to the grid's
+ // end, and stop. AtBody is the transmittance when the march reaches the body's start, before that step,
+ // or the final transmittance when there is no body.
+ private static (float Transmittance, float AtBody) MarchLod(VoidField field, Tile tile, float2 slope, List cells, int bodyStart, float cutoff, Counts counts, ref long steps, ref long bodySteps)
+ {
+ var grid = field.Grid;
+ var body = bodyStart >= 0 ? grid[bodyStart] : float.PositiveInfinity;
+ var transmittance = 1.0f;
+ var k = 0;
+ while (k < cells.Count && transmittance >= cutoff)
+ {
+ var z = grid[cells[k]];
+ var end = cells[k] + 1;
+ for (k++; k < cells.Count && cells[k] == end; k++)
+ end++;
+ while (z < grid[end] && transmittance >= cutoff)
+ {
+ if (z >= body)
+ {
+ var atBody = transmittance;
+ transmittance *= exp(-BodyDepth(field, tile, slope, z, lod: true, counts));
+ steps++;
+ bodySteps++;
+ return (transmittance, atBody);
+ }
+
+ var next = min(min(z + field.LodStep(z), grid[end]), z < body ? body : grid[end]);
+ transmittance *= Integrate(field, tile, slope, z, next, lod: true, counts);
+ steps++;
+ z = next;
+ }
+ }
+
+ return (transmittance, transmittance);
+ }
+
+ // The accuracy oracle before the body: the optical depth over every grid cell below `stop` (never the
+ // mask's, so it trusts no pre-pass), by midpoint quadrature at step h, accumulated in double. Not counted.
+ private static double ConvergedDepth(Scenario field, Tile tile, float2 slope, int stop, float h, bool lod)
+ {
+ var grid = field.Grid;
+ var counts = new Counts();
+ var dir = Math.Sqrt(slope.x * (double)slope.x + slope.y * (double)slope.y + 1.0);
+ var tau = 0.0;
+ for (var c = 0; c < stop; c++)
+ {
+ var samples = Math.Max(1, (int)Math.Ceiling((grid[c + 1] - grid[c]) / h));
+ var dz = (grid[c + 1] - grid[c]) / (double)samples;
+ for (var s = 0; s < samples; s++)
+ tau += field.Density(tile, slope, (float)(grid[c] + (s + 0.5) * dz), lod, counts) * (double)field.Extinction * dz * dir;
+ }
+
+ return tau;
+ }
+
+ internal sealed class RunStats
+ {
+ public readonly Counts Dense = new();
+ public readonly Counts Masked = new();
+ public readonly Counts Probes = new();
+ public readonly Counts Lod = new();
+ public readonly Counts Oracle = new();
+ public long Rays;
+ public long Tiles;
+ public long DenseSteps;
+ public long MaskedSteps;
+ public long LodSteps;
+ public long BodySteps;
+ public long BodyTiles;
+ public long Unfinished;
+ public float MaxDifference;
+ public float MaxLodDifference;
+ public float MaxOracleDifference;
+ public long OracleEmpty;
+ public long OracleEmptySkipped;
+ public double MaxLodDepthError;
+ public double MaxRefDepthError;
+
+ public double PerRay(double value) => value / Rays;
+
+ // Dense cost over the tile march's: each ray's own dense cells plus the tile's probes, shared.
+ public double Ratio => Dense.Cost / (Masked.Cost + Probes.Cost);
+
+ public double LodRatio => Dense.Cost / (Lod.Cost + Probes.Cost);
+
+ // The oracle mask (the cells some ray of the tile samples nonzero) is a ceiling no pre-pass can move.
+ public double OracleCeiling => Dense.Cost / Oracle.Cost;
+
+ public double Efficiency => Oracle.Cost / (Masked.Cost + Probes.Cost);
+
+ // The tile march over the tile march plus its probes: what the probes cost, not what they cull.
+ public double ProbeOverhead => Masked.Cost / (Masked.Cost + Probes.Cost);
+
+ // The oracle-empty cells the pre-pass proves empty; a pre-pass that skips nothing scores 0.
+ public double CullFraction => (double)OracleEmptySkipped / OracleEmpty;
+ }
+
+ ///
+ /// Draws `tiles` tiles, runs the pre-pass once per tile, draws the tile's N^2 ray slopes (one per pixel,
+ /// jittered inside it), builds the oracle mask (the cells whose mid-depth density is nonzero for at least
+ /// one of those slopes), and marches each ray densely over every cell, masked over the tile's dense cells
+ /// and over the oracle cells, all with the early-out at 0.02 and the same Integrate. With lod, the
+ /// footprint-aware march as well, and the agreement at the body start (or the grid's end): the optical
+ /// depth of a cutoff-0 footprint march and of the fixed-step cells, each against ConvergedDepth at
+ /// h = 0.5. The oracle mask and the agreement marches count into throwaway Counts, the oracle march into
+ /// its own, so the cost rows do not move.
+ ///
+ private static RunStats Run(Scenario field, int n, int tiles, int seed, bool lod = false)
+ {
+ var random = new System.Random(seed);
+ var stats = new RunStats();
+ var scratch = new Counts();
+ var grid = field.Grid;
+ for (var t = 0; t < tiles; t++)
+ {
+ var tile = field.DrawTile(random, n);
+ var (cells, bodyStart, _) = MarchTile(field, tile, stats.Probes);
+ Assert.True(field is VoidField || bodyStart < 0, $"{field.Name}: a body reported at cell {bodyStart}; only the void's density can be provably constant");
+ stats.Tiles++;
+ stats.BodyTiles += bodyStart >= 0 ? 1 : 0;
+ var slopes = new List();
+ for (var j = 0; j < n; j++)
+ for (var i = 0; i < n; i++)
+ slopes.Add(tile.PixelSlope(i, j, random.NextSingle(), random.NextSingle()));
+ var inMask = new HashSet(cells);
+ var oracle = new List();
+ for (var c = 0; c < grid.Length - 1; c++)
+ {
+ var mid = (grid[c] + grid[c + 1]) * 0.5f;
+ if (slopes.Any(s => field.Density(tile, s, mid, false, scratch) > 0.0f))
+ {
+ oracle.Add(c);
+ }
+ else
+ {
+ stats.OracleEmpty++;
+ stats.OracleEmptySkipped += inMask.Contains(c) ? 0 : 1;
+ }
+ }
+
+ var stop = bodyStart >= 0 ? bodyStart : grid.Length - 1;
+ foreach (var slope in slopes)
+ {
+ var dense = 1.0f;
+ for (var c = 0; c < grid.Length - 1 && dense >= 0.02f; c++, stats.DenseSteps++)
+ dense *= Integrate(field, tile, slope, grid[c], grid[c + 1], lod: false, stats.Dense);
+ var masked = 1.0f;
+ for (var c = 0; c < cells.Count && masked >= 0.02f; c++, stats.MaskedSteps++)
+ masked *= Integrate(field, tile, slope, grid[cells[c]], grid[cells[c] + 1], lod: false, stats.Masked);
+ var exact = 1.0f;
+ for (var c = 0; c < oracle.Count && exact >= 0.02f; c++)
+ exact *= Integrate(field, tile, slope, grid[oracle[c]], grid[oracle[c] + 1], lod: false, stats.Oracle);
+ stats.Rays++;
+ stats.MaxDifference = MathF.Max(stats.MaxDifference, MathF.Abs(dense - masked));
+ stats.MaxOracleDifference = MathF.Max(stats.MaxOracleDifference, MathF.Abs(dense - exact));
+ if (dense >= 0.02f)
+ stats.Unfinished++;
+ if (lod)
+ {
+ var (marched, _) = MarchLod((VoidField)field, tile, slope, cells, bodyStart, 0.02f, stats.Lod, ref stats.LodSteps, ref stats.BodySteps);
+ stats.MaxLodDifference = MathF.Max(stats.MaxLodDifference, MathF.Abs(dense - marched));
+ long ignoredSteps = 0, ignoredBody = 0;
+ var (_, atBody) = MarchLod((VoidField)field, tile, slope, cells, bodyStart, 0.0f, scratch, ref ignoredSteps, ref ignoredBody);
+ var reference = 1.0;
+ for (var c = 0; c < stop; c++)
+ reference *= Integrate(field, tile, slope, grid[c], grid[c + 1], lod: false, scratch);
+ var tauConv = ConvergedDepth(field, tile, slope, stop, 0.5f, lod: true);
+ var tauLod = -Math.Log(Math.Max(atBody, 1.0e-30f));
+ var tauRef = -Math.Log(Math.Max(reference, 1.0e-300));
+ var scale = Math.Max(tauConv, 1.0e-3);
+ stats.MaxLodDepthError = Math.Max(stats.MaxLodDepthError, Math.Abs(tauLod - tauConv) / scale);
+ stats.MaxRefDepthError = Math.Max(stats.MaxRefDepthError, Math.Abs(tauRef - tauConv) / scale);
+ }
+ }
+ }
+
+ return stats;
+ }
+
+ private void Report(string label, RunStats s) => output.WriteLine(
+ $"IV-REPORT {label}: dense snoise/ray {s.PerRay(s.Dense.Snoise):F2} env/ray {s.PerRay(s.Dense.Envelope):F2} (samples/ray {s.PerRay(s.DenseSteps):F2}); "
+ + $"tile snoise/ray {s.PerRay(s.Masked.Snoise + (double)s.Probes.Snoise):F2} env/ray {s.PerRay(s.Masked.Envelope + (double)s.Probes.Envelope):F2} "
+ + $"(probes/tile env {(double)s.Probes.Envelope / s.Tiles:F1} snoise {(double)s.Probes.Snoise / s.Tiles:F2}; dense cells/ray {s.PerRay(s.MaskedSteps):F2}); "
+ + $"cost ratio {s.Ratio:F2}x, oracle ceiling {s.OracleCeiling:F2}x, efficiency {s.Efficiency:F3}, cull fraction {s.CullFraction:F3} ({s.OracleEmptySkipped}/{s.OracleEmpty}), "
+ + $"probe overhead {s.ProbeOverhead:F3}; max |dT| {s.MaxDifference:R}, oracle {s.MaxOracleDifference:R}; unfinished rays {s.Unfinished}");
+
+ private void ReportVoid(string label, RunStats s) => output.WriteLine(
+ $"IV-REPORT {label}: steps/px ref {s.PerRay(s.DenseSteps):F2}, masked {s.PerRay(s.MaskedSteps):F2}, LOD {s.PerRay(s.LodSteps):F2} ({(double)s.DenseSteps / s.LodSteps:F2}x fewer); "
+ + $"snoise/px ref {s.PerRay(s.Dense.Snoise):F2}, masked {s.PerRay(s.Masked.Snoise):F2}, LOD {s.PerRay(s.Lod.Snoise):F2}; "
+ + $"cost/px ref {s.PerRay(s.Dense.Cost):F2}, masked {s.PerRay(s.Masked.Cost + s.Probes.Cost):F2}, LOD {s.PerRay(s.Lod.Cost + s.Probes.Cost):F2} incl. probes (ref/LOD {s.LodRatio:F2}x, ref/masked {s.Ratio:F2}x); "
+ + $"probes/tile env {(double)s.Probes.Envelope / s.Tiles:F1} snoise {(double)s.Probes.Snoise / s.Tiles:F2}; body tiles {s.BodyTiles}/{s.Tiles}, analytic steps/px {s.PerRay(s.BodySteps):F2}; "
+ + $"depth error at the body start vs converged: LOD {s.MaxLodDepthError:G4}, fixed-step {s.MaxRefDepthError:G4} (ratio {s.MaxLodDepthError / s.MaxRefDepthError:F2}); "
+ + $"end-of-ray max |dT| masked {s.MaxDifference:R}, LOD {s.MaxLodDifference:R}; unfinished rays {s.Unfinished}");
+
+ ///
+ /// The analytic body's soundness pin (ruling operator-analytic-wall-body). Over 2,000 seeded tiles of
+ /// the (e) grid, the shipped void and LOD-far at N in {1, 4, 8, 16}, wherever the pre-pass reports a
+ /// body start b: (1) the gated bound over [grid[b], grid end] is a positive point equal to the reported
+ /// density, every cell from b to the end is in the mask, and at 64 random (pixel, sub-pixel, cell >= b,
+ /// depth in cell) points the full and the truncated Density and the ungated Exact field all equal it
+ /// exactly, no tolerance; (2) BodyDepth, the optical depth MarchLod's one closed-form step uses, equals
+ /// the summed fixed-step depth -ln Integrate of the body cells (no early-out, same cells) within 1e-4
+ /// relative, compared in the depth domain so LOD-far's deep bodies do not underflow. (3) No tile of (a)-(d), at either warp, reports a body: their bounds are never a positive
+ /// point.
+ ///
+ [Fact]
+ public void AnalyticBodyMatchesFineMarch()
+ {
+ var random = new System.Random(0xB0D1);
+ int[] sizes = { 1, 4, 8, 16 };
+ var voids = new List();
+ foreach (var rh in VoidField.HollowRadii)
+ foreach (var offset in VoidField.CameraOffsets)
+ foreach (var ramp in VoidField.Ramps)
+ voids.Add(new VoidField(rh, ramp, offset.Fraction, offset.Name));
+ voids.Add(VoidField.Shipped());
+ voids.Add(VoidField.LodFar());
+ var counts = new Counts();
+ var bodies = new Dictionary();
+ for (var t = 0; t < 2000; t++)
+ {
+ var field = t % 4 == 0 ? voids[^2] : t % 4 == 1 ? voids[^1] : voids[random.Next(voids.Count)];
+ var tile = field.DrawTile(random, sizes[random.Next(sizes.Length)]);
+ var (cells, b, density) = MarchTile(field, tile, counts);
+ if (b < 0)
+ continue;
+ bodies[field.Name] = bodies.GetValueOrDefault(field.Name) + 1;
+ var grid = field.Grid;
+ var last = grid.Length - 1;
+ var rest = field.Bound(tile, grid[b], grid[last], gated: true, counts);
+ Assert.True(rest.x == rest.y && rest.x > 0.0f && rest.x == density, $"{field.Name}: body at cell {b} with density {density:R}, but the bound over the rest is {rest}");
+ Assert.Equal(Enumerable.Range(b, last - b), cells.Skip(cells.Count - (last - b)));
+ for (var i = 0; i < 64; i++)
+ {
+ var slope = tile.PixelSlope(random.Next(tile.N), random.Next(tile.N), random.NextSingle(), random.NextSingle());
+ var c = random.Next(b, last);
+ var z = Uniform(random, grid[c], grid[c + 1]);
+ foreach (var lod in new[] { false, true })
+ {
+ var pointwise = field.Density(tile, slope, z, lod, counts);
+ var exact = field.Exact(tile, slope, z, lod);
+ Assert.True(pointwise == density && exact == density, $"{field.Name}: body from cell {b} claims {density:R}; at slope {slope}, depth {z:R} (cell {c}, lod {lod}) Density is {pointwise:R} and the ungated field {exact:R}");
+ }
+ }
+
+ var ray = tile.PixelSlope(random.Next(tile.N), random.Next(tile.N), random.NextSingle(), random.NextSingle());
+ var depth = BodyDepth(field, tile, ray, grid[b], lod: true, counts);
+ var sum = 0.0;
+ for (var c = b; c < last; c++)
+ sum += -Math.Log(Integrate(field, tile, ray, grid[c], grid[c + 1], lod: true, counts));
+ Assert.True(Math.Abs(depth - sum) <= 1.0e-4 * sum, $"{field.Name}: the body's depth from cell {b}, BodyDepth {depth:R}, differs from the fixed-step cells' summed depth {sum:R} over cells {b}..{last}");
+ }
+
+ output.WriteLine("IV-REPORT analytic body: tiles with a body " + string.Join(", ", bodies.OrderBy(kv => kv.Key).Select(kv => $"{kv.Key} {kv.Value}")));
+ Assert.True(bodies.ContainsKey(VoidField.Shipped().Name) && bodies.ContainsKey(VoidField.LodFar().Name), "no body proven at the shipped void or at LOD-far");
+
+ var scenarios = new (float Warp, Func Make)[]
+ {
+ (FogField.AetheriaWarp, FogField.HeightFog), (FogField.AetheriaWarp, FogField.InsideFog),
+ (FogField.AetheriaWarp, FogField.Wells), (SlabField.SlabWarp, SlabField.Create),
+ };
+ foreach (var (warp, make) in scenarios)
+ foreach (var w in new[] { 0.0f, warp })
+ {
+ var field = make(w);
+ for (var t = 0; t < 250; t++)
+ {
+ var (_, b, _) = MarchTile(field, field.DrawTile(random, sizes[random.Next(sizes.Length)]), counts);
+ Assert.True(b < 0, $"{field.Name} warp {w:R}: a body reported at cell {b}");
+ }
+ }
+ }
+
+ ///
+ /// The proving ground (docs/cultmath-interval-ground-cut.md, "Pass 3" and "Pass 4"). Scenarios
+ /// (a) height fog, (b) inside the fog, (c) Aetheria's wells and (d) r1's uniform slab, each over 64
+ /// tiles at N in {1, 4, 8, 16} and warp in {0, D}: the tile march integrates exactly the dense march's
+ /// cells, so the transmittance agrees within 1e-3 (a difference would be a skipped cell that was not
+ /// empty), the oracle march (each ray over the oracle mask: the cells whose mid-depth density is nonzero
+ /// for at least one of the tile's rays) agrees with the dense march exactly, and no (a)-(d) tile reports
+ /// a body. (e) the void: the grid over hollow radius, camera offset and ramp width (32 tiles at N = 8),
+ /// LOD-far (32 tiles) and the one shipped void (200 tiles), three marches each, the footprint-aware one
+ /// taking the analytic body step. 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; the oracle ceiling, the efficiency against it and the cull fraction are
+ /// printed; (e) at the shipped void the masked fixed-step march agrees with the reference exactly, the
+ /// footprint-aware march's optical depth at the body start is off a converged march (h = 0.5, every
+ /// cell) by at most twice the fixed-step march's error, and it takes at most half the reference's
+ /// steps per pixel. A shortfall skips naming saving-2x-scenarios and every shortfall; the fields, the
+ /// grids and the envelope cost are not tuned toward it.
+ ///
+ [Fact]
+ public void IntervalSkipHalvesEvaluations()
+ {
+ var ratios = new Dictionary();
+ var cullFloor = double.MaxValue;
+ var scenarios = new (string Key, float Warp, Func Make)[]
+ {
+ ("a", FogField.AetheriaWarp, FogField.HeightFog),
+ ("b", FogField.AetheriaWarp, FogField.InsideFog),
+ ("c", FogField.AetheriaWarp, FogField.Wells),
+ ("d", SlabField.SlabWarp, SlabField.Create),
+ };
+ foreach (var (key, warp, make) in scenarios)
+ foreach (var n in new[] { 1, 4, 8, 16 })
+ foreach (var w in new[] { 0.0f, warp })
+ {
+ var field = make(w);
+ var stats = Run(field, n, 64, 0x5A7E + n);
+ Report($"{field.Name} N={n} warp={w:R}", stats);
+ Assert.True(stats.MaxDifference <= 1.0e-3f, $"{field.Name} N={n} warp={w:R}: transmittance differs by {stats.MaxDifference:R}");
+ Assert.True(stats.MaxOracleDifference == 0.0f, $"{field.Name} N={n} warp={w:R}: the oracle march differs from the dense march by {stats.MaxOracleDifference:R}");
+ if (key == "d")
+ output.WriteLine($"IV-REPORT (d) N={n} warp={w:R}: ungated dense snoise/ray {3.0 * stats.PerRay(stats.DenseSteps):F2} (three per dense sample, as r1 counted)");
+ if (n == 8)
+ ratios[key] = Math.Min(ratios.GetValueOrDefault(key, double.MaxValue), key == "c" ? stats.ProbeOverhead : stats.Ratio);
+ if (n == 8 && key == "c")
+ {
+ Assert.True(stats.CullFraction >= 0.5, $"(c) warp={w:R}: the pre-pass proves {stats.CullFraction:F3} of the oracle-empty cells empty, under the 0.5 floor (ruling wells-overhead-plus-floor)");
+ cullFloor = Math.Min(cullFloor, stats.CullFraction);
+ }
+ }
+
+ foreach (var rh in VoidField.HollowRadii)
+ foreach (var offset in VoidField.CameraOffsets)
+ foreach (var ramp in VoidField.Ramps)
+ {
+ var field = new VoidField(rh, ramp, offset.Fraction, offset.Name);
+ ReportVoid($"{field.Name} N=8", Run(field, 8, 32, 0x701D, lod: true));
+ }
+
+ var lodFar = VoidField.LodFar();
+ ReportVoid($"{lodFar.Name} N=8 (LOD far)", Run(lodFar, 8, 32, 0x701D, lod: true));
+
+ var shipped = VoidField.Shipped();
+ var headline = Run(shipped, 8, 200, 0x5417, lod: true);
+ ReportVoid($"{shipped.Name} N=8 (headline, 200 tiles)", headline);
+ Assert.True(headline.MaxDifference == 0.0f, $"the masked fixed-step march differs from the reference by {headline.MaxDifference:R}");
+ Assert.True(headline.MaxLodDepthError <= 2.0 * headline.MaxRefDepthError, $"at the body start the footprint-aware march's optical depth is off the converged march's by {headline.MaxLodDepthError:G6} relative, more than twice the fixed-step march's {headline.MaxRefDepthError:G6}");
+ var steps = (double)headline.DenseSteps / headline.LodSteps;
+
+ var shortfalls = new List();
+ if (ratios["a"] < 2.0)
+ shortfalls.Add($"(a) saves {ratios["a"]:F2}x combined cost, under 2x");
+ if (ratios["c"] < 0.90)
+ shortfalls.Add($"(c) probe overhead {ratios["c"]:F3}, under 0.90");
+ if (steps < 2.0)
+ shortfalls.Add($"(e) takes {steps:F2}x fewer steps per pixel at the shipped void, under 2x");
+ output.WriteLine($"IV-REPORT contracts at N=8 (lower of warp 0 and D): (a) {ratios["a"]:F2}x cost; (c) probe overhead {ratios["c"]:F3} (at least 0.90), cull fraction {cullFloor:F3} (at least 0.5, hard); "
+ + $"(e) {steps:F2}x fewer steps/px, depth error at the body start LOD {headline.MaxLodDepthError:G4} vs fixed-step {headline.MaxRefDepthError:G4} (at most 2x)");
+ if (shortfalls.Count > 0)
+ Assert.Skip("saving-2x-scenarios: " + string.Join("; ", shortfalls));
+ }
+}
diff --git a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl
new file mode 100644
index 00000000..b6b16665
--- /dev/null
+++ b/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
diff --git a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl.meta b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl.meta
new file mode 100644
index 00000000..023d218c
--- /dev/null
+++ b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl.meta
@@ -0,0 +1,7 @@
+fileFormatVersion: 2
+guid: 18b031192f344246ac523cea014b8f24
+ShaderIncludeImporter:
+ externalObjects: {}
+ userData:
+ assetBundleName:
+ assetBundleVariant:
diff --git a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl
index af2d3741..da3362a6 100644
--- a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl
+++ b/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)
{