From 5f6bbfd4b51a84af0854b1845fd03d10fa7cdbc3 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Sat, 3 Oct 2026 11:37:30 +0200 Subject: [PATCH] CultMath intervals r3: iv_frustum_ball widened by its derived float32 rounding bound; extreme-family enclosure The radius gains (|c|_1 + radius) 2^-20 in C# and both HLSL copies; the doc's strictness claim leaves. TileBallEnclosesEveryRaySegment measures in double over the r2 domain and seven extreme families with outward corner flows, pins the radius two-sided, and counts misses per family. The noise ball enclosure tests draw radii to 1e3. Co-Authored-By: Claude Opus 5.5 --- packages/cultmath/README.md | 3 +- packages/cultmath/docs/design.md | 26 ++- .../cultmath/shaders/CultMath.Interval.hlsl | 8 +- .../cultmath/src/CultMath/math.Interval.cs | 19 +- .../tests/CultMath.Tests/NoiseBoundTests.cs | 165 +++++++++++++----- .../Shaders/CultMath.Interval.hlsl | 8 +- 6 files changed, 173 insertions(+), 56 deletions(-) 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