diff --git a/packages/cultmath/README.md b/packages/cultmath/README.md
index aba9c1a5..81d15859 100644
--- a/packages/cultmath/README.md
+++ b/packages/cultmath/README.md
@@ -39,7 +39,8 @@ shader semantics.
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. Every one returns an interval
+ 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
diff --git a/packages/cultmath/docs/design.md b/packages/cultmath/docs/design.md
index aa348744..77688dc5 100644
--- a/packages/cultmath/docs/design.md
+++ b/packages/cultmath/docs/design.md
@@ -409,12 +409,26 @@ 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)`, and the two vectors are never parallel, so the bound is strict. The
-caller rotates and translates the centre into world space and scales centre and
-radius by its noise frequency. `NoiseBoundTests.TileBallEnclosesEveryRaySegment`
-checks 2,000 seeded tiles x 64 points, tile corners at both depth ends with
-full-length flows included, with no tolerance, and pins the ball to the
-derivation, so a ball grown past it fails too. A moving camera needs nothing
+`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. At the shipped scales (`|c|_1` near 600, `r` near 5) it adds
+`6e-4`. 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.
diff --git a/packages/cultmath/shaders/CultMath.Interval.hlsl b/packages/cultmath/shaders/CultMath.Interval.hlsl
index 732ef3a0..b6b16665 100644
--- a/packages/cultmath/shaders/CultMath.Interval.hlsl
+++ b/packages/cultmath/shaders/CultMath.Interval.hlsl
@@ -106,13 +106,17 @@ float2 cultmath_iv_fbm_ball(float3 centre, float radius, int octaves, float lacu
// 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. Derivation in math.Interval.cs.
+// 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;
- return float4(centreSlope.x * zm, centreSlope.y * zm, zm, radius);
+ radius += (abs(cx) + abs(cy) + abs(zm) + radius) * 9.5367431640625e-7;
+ return float4(cx, cy, zm, radius);
}
#endif
diff --git a/packages/cultmath/src/CultMath/math.Interval.cs b/packages/cultmath/src/CultMath/math.Interval.cs
index 547c4acc..91e99f87 100644
--- a/packages/cultmath/src/CultMath/math.Interval.cs
+++ b/packages/cultmath/src/CultMath/math.Interval.cs
@@ -175,16 +175,25 @@ public static partial class math
/// 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
- /// (NoiseBoundTests.TileBallEnclosesEveryRaySegment). The two vectors are never parallel, so the
- /// triangle inequality is strict. The caller rotates and translates the centre into world space (the
- /// radius is unchanged) and scales both by its noise frequency.
+ /// 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;
- return new float4(centreSlope.x * zm, centreSlope.y * zm, zm, radius);
+ radius += (abs(cx) + abs(cy) + abs(zm) + radius) * 9.5367431640625e-7f;
+ return new float4(cx, cy, zm, radius);
}
}
diff --git a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs
index b1c39d14..d895828e 100644
--- a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs
+++ b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs
@@ -13,7 +13,11 @@ namespace CultMath.Tests;
/// 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.
+/// 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: w(z0) for w(z1) in VoidField.WeightBound must fail the truncated
+/// check of EnvelopeBoundEnclosesDensity on the LOD-far configuration; CarveRadius for HollowRadius in
+/// the lo test of VoidField.Compose must fail AnalyticBodyMatchesFineMarch.
///
public sealed class NoiseBoundTests
{
@@ -55,7 +59,7 @@ public sealed class NoiseBoundTests
for (var b = 0; b < 2000; b++)
{
var centre = RandomCentre(random);
- var radius = LogUniform(random, 1.0e-3f, 2.0f);
+ var radius = LogUniform(random, 1.0e-3f, 1.0e3f);
var bound = iv_snoise_ball(centre, radius);
for (var i = 0; i < 64; i++)
{
@@ -72,7 +76,7 @@ public sealed class NoiseBoundTests
for (var b = 0; b < 2000; b++)
{
var centre = RandomCentre(random);
- var radius = LogUniform(random, 1.0e-3f, 2.0f);
+ 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);
@@ -114,53 +118,134 @@ public sealed class NoiseBoundTests
// ---- 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 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: 2,000
- /// seeded tiles (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, the first eight at
- /// a corner of the footprint at either end of the segment with a full-length flow. It also pins the
- /// ball to the derivation: centre and radius equal the formula evaluated in double to one part in
- /// 1e5, so a ball that grows (sound but wasteful) fails too. 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.
+ /// 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);
- int[] sizes = { 1, 4, 8, 16 };
- float[] depths = { 12.0f, 1000.0f, 2400.0f };
- for (var t = 0; t < 2000; t++)
+ 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)
{
- 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);
- 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 centre = new float3(ball.x, ball.y, ball.z);
-
- var zm = (z0 + (double)z1) / 2.0;
- 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;
- Assert.True(Math.Abs(ball.w - radius) <= 1.0e-5 * radius, $"radius {ball.w:R}, derivation {radius:R}");
- Assert.True(Math.Abs(ball.x - slope.x * zm) + Math.Abs(ball.y - slope.y * zm) + Math.Abs(ball.z - zm) <= 1.0e-5 * zm, $"centre {centre}, derivation ({slope.x * zm:R}, {slope.y * zm:R}, {zm:R})");
-
- for (var i = 0; i < 64; i++)
+ misses[family] = 0;
+ pins[family] = 0;
+ for (var t = 0; t < 2000; t++)
{
- 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 flow = UnitVector(random) * (i < 8 ? warp : warp * MathF.Sqrt(random.NextSingle()));
- var p = new float3(m.x * z, m.y * z, z) + flow;
- Assert.True(length(p - centre) <= ball.w, $"iv_frustum_ball({slope}, {z0:R}, {z1:R}, {footprint:R}, {warp:R}) = {ball} misses {p} (slope {m}, depth {z:R}), distance {length(p - centre):R}");
+ 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));
}
///
diff --git a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl
index 732ef3a0..b6b16665 100644
--- a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl
+++ b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl
@@ -106,13 +106,17 @@ float2 cultmath_iv_fbm_ball(float3 centre, float radius, int octaves, float lacu
// 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. Derivation in math.Interval.cs.
+// 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;
- return float4(centreSlope.x * zm, centreSlope.y * zm, zm, radius);
+ radius += (abs(cx) + abs(cy) + abs(zm) + radius) * 9.5367431640625e-7;
+ return float4(cx, cy, zm, radius);
}
#endif