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>
This commit is contained in:
Meta Vens 2026-10-16 04:49:13 +02:00
commit 5ccbb09487
4 changed files with 179 additions and 1 deletions

View file

@ -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

View file

@ -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);

View file

@ -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:

View file

@ -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