cultmath tests: band read families with an allowance check, a reach test, GradDouble

Four coordinate-band snoise_grad / snoise_hessian families with Verify AssertSpanAllowance,
SpanAllowanceCoversReadsAlongTheSpan, NoiseBoundTests.GradDouble and its pin, design.md
paragraph. The fixture and CultMath.glsl are regenerated in the next commit; this one compiles.

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

packages/cultmath/docs/design.md

@@ -243,6 +243,13 @@ the main families' reach is at least about 1e-3, where the second-order remainde
 allowance still passes. `GoldenFixtureMatchesCSharp` requires each of the two to hold cases that fail on C#'s own results once
 the allowance is removed.
 
+Four more appended families, the coordinate-band `snoise_grad` and `snoise_hessian` reads (256 cases each, the magnitude of the
+point log-uniform from 1e-4 to 1e2 and from 1e2 to 4e6), carry an allowance `check`: with `S = |x|_1` and the allowances
+`math.span_slack(math.span_rounding(1, x, zero stage, 0, 0, 0))` (the owners of the span stage's allowances in C#, HLSL and GLSL),
+the result must lie within the allowance of the kernel's exact value, gradient or Hessian in double, no tolerance. The two small
+families must also fail it in at least 16 cases once the absolute term is removed. `GoldenFixtureMatchesCSharp` applies the check
+to C#'s own results, and a device capture is judged by the same `Verify`.
+
 `GoldenFixtureMatchesCSharp` evaluates C# on the committed fixture's own
 arguments. Arguments are drawn only when the fixture is regenerated
 (`CULTMATH_WRITE_GLSL=1`), because the draw itself calls `MathF.Pow`, whose

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

@@ -408,6 +408,74 @@ static int4 clamp(int4 x, int4 a, int4 b) { if (a.x > b.x || a.y > b.y || a.z >
         return $"least enclosure slack {enclosure:E3} of the radius, least slack above R0 + d {lower:E3} of R0, below its widening {upper:F3} of 2^-19 (R0 + |c0|_1 + spread)";
     }
 
+    private const string SpanAllowanceRead =
+        "Read at the case's point x, with S = |x|_1 and (e_n, e_g, e_h) = math.span_slack(math.span_rounding(1, x, zero stage, 0, 0, 0)), " +
+        "against the kernel's exact values in double: |n - n_exact| <= e_n and |g - g_exact|_2 <= e_g for snoise_grad, " +
+        "||H - H_exact||_F <= e_h for snoise_hessian, no tolerance. " +
+        "Report the largest fraction of each allowance per decade of S and the largest absolute error at S < 1e-2.";
+
+    // SpanAllowanceRead applied to a family's results (C#'s own, or a device's): every read is within the
+    // allowance math.span_slack gives at its point, the owner's own numbers and no restated literal. A small
+    // family (bands 1e-4 to 1e2) must also fail the check once the absolute term is removed in at least 16
+    // cases, and hold at least 48 cases with S <= 1e-2: the allowance's absolute term is what decides there.
+    private static string AssertSpanAllowance(JsonArray cases, bool hessian, bool small)
+    {
+        var absolute = math.span_slack(0.0f);
+        var fractions = new SortedDictionary<int, (double Value, double Gradient)>();
+        double smallestError = 0.0;
+        int absoluteDecides = 0, nearOrigin = 0;
+        foreach (var entry in cases)
+        {
+            var sides = ((string)entry!).Split(" -> ");
+            double[] x = Floats(sides[0]), result = Floats(sides[1]);
+            var s = Math.Abs(x[0]) + Math.Abs(x[1]) + Math.Abs(x[2]);
+            var slack = math.span_slack(math.span_rounding(1.0f, new float3((float)x[0], (float)x[1], (float)x[2]), default, new float3(0.0f, 0.0f, 0.0f), 0.0f, 0.0f));
+            double valueError = 0.0, gradientError, valueFraction = 0.0, gradientFraction;
+            bool needsAbsolute;
+            if (hessian)
+            {
+                var exact = NoiseBoundTests.HessianDouble(x[0], x[1], x[2]);
+                var sum = 0.0;
+                for (var i = 0; i < 9; i++)
+                    sum += (result[i] - exact[i]) * (result[i] - exact[i]);
+                gradientError = Math.Sqrt(sum);
+                Assert.True(gradientError <= slack.z, $"snoise_hessian case {entry}: ||H - H_exact||_F is {gradientError:R}, past the allowance {slack.z:R}");
+                gradientFraction = gradientError / slack.z;
+                needsAbsolute = gradientError > slack.z - absolute.z;
+            }
+            else
+            {
+                var exact = NoiseBoundTests.GradDouble(x[0], x[1], x[2]);
+                valueError = Math.Abs(result[3] - NoiseBoundTests.ValueDouble(x[0], x[1], x[2]));
+                gradientError = Math.Sqrt((result[0] - exact[0]) * (result[0] - exact[0]) + (result[1] - exact[1]) * (result[1] - exact[1]) + (result[2] - exact[2]) * (result[2] - exact[2]));
+                Assert.True(valueError <= slack.x, $"snoise_grad case {entry}: |n - n_exact| is {valueError:R}, past the allowance {slack.x:R}");
+                Assert.True(gradientError <= slack.y, $"snoise_grad case {entry}: |g - g_exact|_2 is {gradientError:R}, past the allowance {slack.y:R}");
+                valueFraction = valueError / slack.x;
+                gradientFraction = gradientError / slack.y;
+                needsAbsolute = valueError > slack.x - absolute.x || gradientError > slack.y - absolute.y;
+            }
+
+            var decade = (int)Math.Floor(Math.Log10(s));
+            fractions.TryGetValue(decade, out var largest);
+            fractions[decade] = (Math.Max(largest.Value, valueFraction), Math.Max(largest.Gradient, gradientFraction));
+            if (needsAbsolute)
+                absoluteDecides++;
+            if (s < 1.0e-2)
+                smallestError = Math.Max(smallestError, Math.Max(valueError, gradientError));
+            if (s <= 1.0e-2)
+                nearOrigin++;
+        }
+
+        if (small)
+        {
+            Assert.True(absoluteDecides >= 16, $"only {absoluteDecides} of {cases.Count} cases fail once the absolute term is removed: the family cannot tell a device that needs it from one that does not");
+            Assert.True(nearOrigin >= 48, $"only {nearOrigin} of {cases.Count} cases have |x|_1 <= 1e-2");
+        }
+
+        var report = string.Join(", ", fractions.Select(f => hessian ? $"1e{f.Key}: {f.Value.Gradient:F3}" : $"1e{f.Key}: {f.Value.Value:F3}/{f.Value.Gradient:F3}"));
+        return $"largest fraction of {(hessian ? "e_h" : "e_n/e_g")} per decade of |x|_1: {report}; largest absolute error at |x|_1 < 1e-2 {smallestError:E3}; {absoluteDecides} of {cases.Count} cases need the absolute term";
+    }
+
     // A committed case's arguments, decoded into the C# method's parameter types in Hex's field order.
     private static object[] Arguments(MethodInfo method, string hex)
     {
@@ -575,7 +643,8 @@ static int4 clamp(int4 x, int4 a, int4 b) { if (a.x > b.x || a.y > b.y || a.z >
     // Every family the site's WebGL2 readback evaluates, in a fixed order. A new family is appended, so the
     // seeded stream and every earlier family's cases stay as they were. pcg3d and pcg4d are "exact"; the
     // float families are ulp-bounded, with the bound measured on the device by the consumer, and
-    // iv_frustum_ball also carries an enclosure check. phacelle is in CultMath.Phacelle.glsl, so its
+    // iv_frustum_ball also carries an enclosure check; the four coordinate-band snoise_grad and snoise_hessian
+    // families carry an allowance check, applied to C#'s own results here and to a device's by the device judge. phacelle is in CultMath.Phacelle.glsl, so its
     // evaluator concatenates that file after CultMath.glsl; the site, which does not vendor it, skips it.
     private static readonly Family[] Families = new Family[]
     {
@@ -720,8 +789,24 @@ static int4 clamp(int4 x, int4 a, int4 b) { if (a.x > b.x || a.y > b.y || a.z >
         // midpoint a Point3, one span in four a point.
         new Family("SpanNoise cultmath_span_fbm(SpanPath, SpanStage, int, float, float)", "ulp-bounded", SpanArguments),
         new Family("SpanNoise cultmath_span_ridged(SpanPath, SpanStage, int, float, float)", "ulp-bounded", SpanArguments),
+        // The kernel's reads over coordinate bands, the magnitude log-uniform in the band: the allowance of
+        // math.span_slack is checked against the exact kernel at every case (SpanAllowanceRead).
+        new Family("vec4 cultmath_snoise_grad(vec3)", "ulp-bounded, coordinate bands 1e-4 to 1e2", r => new object[] { BandPoint3(r, 1.0e-4, 1.0e2) }, SpanAllowanceRead, Verify: cases => AssertSpanAllowance(cases, hessian: false, small: true)),
+        new Family("vec4 cultmath_snoise_grad(vec3)", "ulp-bounded, coordinate bands 1e2 to 4e6", r => new object[] { BandPoint3(r, 1.0e2, 4.0e6) }, SpanAllowanceRead, Verify: cases => AssertSpanAllowance(cases, hessian: false, small: false)),
+        new Family("mat3 cultmath_snoise_hessian(vec3)", "ulp-bounded, coordinate bands 1e-4 to 1e2", r => new object[] { BandPoint3(r, 1.0e-4, 1.0e2) }, SpanAllowanceRead, Verify: cases => AssertSpanAllowance(cases, hessian: true, small: true)),
+        new Family("mat3 cultmath_snoise_hessian(vec3)", "ulp-bounded, coordinate bands 1e2 to 4e6", r => new object[] { BandPoint3(r, 1.0e2, 4.0e6) }, SpanAllowanceRead, Verify: cases => AssertSpanAllowance(cases, hessian: true, small: false)),
     }).ToArray();
 
+    // A point of magnitude m, m log-uniform in [lo, hi]: one component at +-m, the other two uniform in [-m, m],
+    // the axis of the large one random.
+    private static float3 BandPoint3(System.Random r, double lo, double hi)
+    {
+        var m = lo * Math.Pow(hi / lo, r.NextDouble());
+        var v = new double[] { (r.NextDouble() * 2.0 - 1.0) * m, (r.NextDouble() * 2.0 - 1.0) * m, (r.NextDouble() * 2.0 - 1.0) * m };
+        v[r.Next(3)] = r.Next(2) == 0 ? m : -m;
+        return new float3((float)v[0], (float)v[1], (float)v[2]);
+    }
+
     private static object[] SpanArguments(System.Random r)
     {
         var stage = SpanStageTests.Stage.Draw(r);

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

@@ -722,6 +722,44 @@ public sealed class NoiseBoundTests
         return total;
     }
 
+    // The gradient of snoise at v in double: 105 sum (m^4 p - 8 m^3 s x) over Corners, the derivative of
+    // ValueDouble (every corner's x moves one for one with v, dm = -2 x, ds = p).
+    internal static double[] GradDouble(double vx, double vy, double vz)
+    {
+        var (x, p, _) = Corners(vx, vy, vz);
+        var gradient = new double[3];
+        for (var c = 0; c < 4; c++)
+        {
+            var m = 0.5 - (x[c][0] * x[c][0] + x[c][1] * x[c][1] + x[c][2] * x[c][2]);
+            if (m <= 0.0)
+                continue;
+            var s = p[c][0] * x[c][0] + p[c][1] * x[c][1] + p[c][2] * x[c][2];
+            for (var i = 0; i < 3; i++)
+                gradient[i] += 105.0 * (m * m * m * m * p[c][i] - 8.0 * m * m * m * s * x[c][i]);
+        }
+
+        return gradient;
+    }
+
+    /// <summary>
+    /// GradDouble is the kernel's gradient: within float rounding of snoise_grad (L1 distance at most 2e-4)
+    /// at 300 seeded points of [-20, 20]^3, the way ExactReferenceIsTheKernel pins ValueDouble.
+    /// </summary>
+    [Fact]
+    public void GradDoubleIsTheKernelsGradient()
+    {
+        var random = new System.Random(0x6BAD);
+        double Coordinate() => random.NextDouble() * 40.0 - 20.0;
+        for (var k = 0; k < 300; k++)
+        {
+            var p = new float3((float)Coordinate(), (float)Coordinate(), (float)Coordinate());
+            var exact = GradDouble(p.x, p.y, p.z);
+            var kernel = snoise_grad(p);
+            var distance = Math.Abs(exact[0] - kernel.x) + Math.Abs(exact[1] - kernel.y) + Math.Abs(exact[2] - kernel.z);
+            Assert.True(distance <= 2.0e-4, $"GradDouble({p}) is {exact[0]:R}, {exact[1]:R}, {exact[2]:R}, snoise_grad {kernel.x:R}, {kernel.y:R}, {kernel.z:R}");
+        }
+    }
+
     private static double Delta(int a, int b) => a == b ? 1.0 : 0.0;
 
     // The symmetric third-derivative tensor of snoise at p, entry (i * 3 + j) * 3 + k, in closed form:

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

@@ -752,6 +752,54 @@ public sealed class SpanStageTests
         Assert.True(crossing >= 100, $"only {crossing} of 600 midpoints had an exact sign opposite to the float one: the check has no teeth");
     }
 
+    /// <summary>
+    /// The allowance a span computes covers the float kernel's read at every point of its path: the identity
+    /// stage, the midpoint within |q|_inf 1, a velocity of random direction and speed 1 to 1e3 (log-uniform),
+    /// a halfLength 1 to 1e3, so path points lie out to about 1e6. At 16 points of the straight path
+    /// q + v t the value, the gradient's length and the Hessian's Frobenius norm of the float kernel are
+    /// within math.span_slack(math.span_rounding(...)) of the exact kernel's, the reach term of S_o included.
+    /// Teeth: the largest value error exceeds the allowance computed with nu = 0 by more than 10 times.
+    /// </summary>
+    [Fact]
+    public void SpanAllowanceCoversReadsAlongTheSpan()
+    {
+        var random = new System.Random(0x5EA4);
+        double worstValue = 0.0, worstGradient = 0.0, worstHessian = 0.0, teeth = 0.0;
+        for (var k = 0; k < 200; k++)
+        {
+            var q = Vector(random, -1.0f, 1.0f);
+            var speed = LogUniform(random, 1.0f, 1.0e3f);
+            var velocity = UnitVector(random) * speed;
+            var halfLength = LogUniform(random, 1.0f, 1.0e3f);
+            var slack = span_slack(span_rounding(1.0f, q, SpanStage.Identity, None, speed, halfLength));
+            var midpointSlack = span_slack(span_rounding(1.0f, q, SpanStage.Identity, None, 0.0f, halfLength));
+            for (var j = 0; j < 16; j++)
+            {
+                var x = q + velocity * (halfLength * (j / 7.5f - 1.0f));
+                var value = NoiseBoundTests.ValueDouble(x.x, x.y, x.z);
+                var gradient = NoiseBoundTests.GradDouble(x.x, x.y, x.z);
+                var hessian = NoiseBoundTests.HessianDouble(x.x, x.y, x.z);
+                var read = snoise_grad(x);
+                var h = snoise_hessian(x);
+                var valueError = Math.Abs(read.w - value);
+                var gradientError = Math.Sqrt(Math.Pow(read.x - gradient[0], 2) + Math.Pow(read.y - gradient[1], 2) + Math.Pow(read.z - gradient[2], 2));
+                var hessianError = Math.Sqrt(Math.Pow(h._m00 - hessian[0], 2) + Math.Pow(h._m01 - hessian[1], 2) + Math.Pow(h._m02 - hessian[2], 2)
+                    + Math.Pow(h._m10 - hessian[3], 2) + Math.Pow(h._m11 - hessian[4], 2) + Math.Pow(h._m12 - hessian[5], 2)
+                    + Math.Pow(h._m20 - hessian[6], 2) + Math.Pow(h._m21 - hessian[7], 2) + Math.Pow(h._m22 - hessian[8], 2));
+                Assert.True(valueError <= slack.x, $"midpoint {q}, speed {speed:R}, halfLength {halfLength:R}, point {x}: the value error {valueError:R} is past the allowance {slack.x:R}");
+                Assert.True(gradientError <= slack.y, $"midpoint {q}, speed {speed:R}, halfLength {halfLength:R}, point {x}: the gradient error {gradientError:R} is past the allowance {slack.y:R}");
+                Assert.True(hessianError <= slack.z, $"midpoint {q}, speed {speed:R}, halfLength {halfLength:R}, point {x}: the Hessian error {hessianError:R} is past the allowance {slack.z:R}");
+                worstValue = Math.Max(worstValue, valueError / slack.x);
+                worstGradient = Math.Max(worstGradient, gradientError / slack.y);
+                worstHessian = Math.Max(worstHessian, hessianError / slack.z);
+                teeth = Math.Max(teeth, valueError / midpointSlack.x);
+            }
+        }
+
+        output.WriteLine($"IV-REPORT span allowance: largest fraction of e_n {worstValue:F3}, e_g {worstGradient:F3}, e_h {worstHessian:F3}; largest value error over the allowance without the reach {teeth:F1}");
+        Assert.True(teeth > 10.0, $"the largest value error is only {teeth:F1} times the allowance without the reach: the check has no teeth");
+    }
+
     /// <summary>
     /// The fold's acceleration term, alone, at the identity stage and under a transform (frequency 2.5, stretch
     /// 0.5, shift and offset): the midpoint velocity is zero and the acceleration points along the gradient