GameCult/CultLib
Watch
9
0
Fork
You've already forked CultLib
 0
Code Issues Pull requests 1 Projects Releases Packages Wiki Activity Actions

Affine forms: math.Affine.cs (af_* ops, af_snoise, af_fbm, af_frustum_*), SNOISE2_LIPSCHITZ and iv_snoise_ball(float2), AffineTests

Browse source 
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Meta Vens 2026-10-03 21:41:08 +02:00
commit ad7e3079d4
5 changed files with 706 additions and 1 deletions
Download patch file Download diff file Expand all files Collapse all files

224
packages/cultmath/src/CultMath/math.Affine.cs Normal file
View file

@ -0,0 +1,224 @@
namespace CultMath;
public static partial class math
{
// Reduced affine forms with one shared symbol (design.md, "Affine forms"): float3(x0, a, e) is the
// set { x0 + a eps + e delta : delta in [-1, 1] } at the region's shared eps in [-1, 1], with
// e >= 0. Every form over one region shares that eps (for a screen-tile slice, the depth
// parameter), and delta is private to the form. This is Messine's AF1 as Gamito and Maddock use it
// for gradient noise (t0 + t1 e1 + t2 e2), in one register, as an interval is a float2.
//
// Every af_* op returns a form whose set, at each eps, contains f evaluated pointwise in float32 at
// every point of its operands' sets at that same eps. Rounding rule: each op folds its own float32
// rounding into e as (S) 2^-20, where S is the sum of the magnitudes the op's roundings are
// relative to (named in each comment). With u = 2^-24, every op below needs at most 8u S (its
// component roundings, the pointwise result's own rounding, and the roundings of computing e), so
// 2^-20 = 16u covers it twice over; the factor of two is the allowance for a device whose + and *
// are faithful (within one ulp) rather than correctly rounded, and fma contraction only removes
// roundings. Products also add 2^-126, the most any underflow to a subnormal can lose. The
// enclosure tests evaluate the form's set in double and f in float32 with no tolerance
// (AffineTests.EveryOpEnclosesItsPoints). No op takes a compound type or a reference parameter: each is a
// plain function of float values, so it mirrors into the HLSL common subset.
// An upper bound on the spectral norm of the Hessian of snoise(float3). Provenance:
// NoiseBoundTests.MeasureHessian (slow, explicit) takes the largest Hessian norm (central
// differences of snoise_grad, largest absolute eigenvalue) over 1e6 seeded points, refines the 1e4
// largest by ascent, and multiplies the refined maximum by 1.10;
// NoiseBoundTests.HessianConstantPinsSampledCurvature pins this value against a fast re-measure and
// the ascent from MeasureHessian's witness. Empirical with a margin, not a proof; the enclosure
// tests are the defence. A change to the snoise kernel must re-run MeasureHessian and re-pin it.
public const float SNOISE_HESSIAN = 68.94f;
/// <summary>The point form (x, 0, 0): the value x at every eps. Exact.</summary>
public static float3 af_point(float value) => new(value, 0.0f, 0.0f);
/// <summary>
/// The shared symbol itself, (x0, a, 0): the value x0 + a eps. It is how a region's parameter enters
/// (for a tile slice, depth z = z_m + h eps), so it is exact by definition: eps is defined through it.
/// </summary>
public static float3 af_symbol(float x0, float a) => new(x0, a, 0.0f);
/// <summary>
/// Every value of the interval [lo, hi], uncorrelated with eps: ((lo + hi) / 2, 0, (hi - lo) / 2).
/// Rounding: x0 and e each round once, and a value of the interval is within e of x0 only up to
/// those; S = |x0| + e.
/// </summary>
public static float3 af_from_iv(float2 interval)
{
var x0 = (interval.x + interval.y) * 0.5f;
var e = (interval.y - interval.x) * 0.5f;
return new float3(x0, 0.0f, e + (abs(x0) + e) * 9.5367431640625e-7f);
}
/// <summary>
/// The interval of every value of the form over every eps: x0 -+ (|a| + e), the bridge to iv_*. The
/// endpoints round inward by at most u (|x0| + |a| + e) each, so each is moved outward by
/// (|x0| + |a| + e) 2^-20.
/// </summary>
public static float2 af_range(float3 x)
{
var r = abs(x.y) + x.z;
var w = (abs(x.x) + r) * 9.5367431640625e-7f;
return new float2(x.x - r - w, x.x + r + w);
}
/// <summary>
/// Encloses x + y: (x0 + y0, a_x + a_y, e_x + e_y). The shared symbol keeps the two forms
/// correlated, so af_sub(x, x) has a = 0. Rounding: S = |x0| + |a| + e of the result.
/// </summary>
public static float3 af_add(float3 x, float3 y)
{
var x0 = x.x + y.x;
var a = x.y + y.y;
var e = x.z + y.z;
return new float3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7f);
}
/// <summary>Encloses x - y: (x0 - y0, a_x - a_y, e_x + e_y). Rounding as af_add.</summary>
public static float3 af_sub(float3 x, float3 y)
{
var x0 = x.x - y.x;
var a = x.y - y.y;
var e = x.z + y.z;
return new float3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7f);
}
/// <summary>Encloses -x: (-x0, -a, e). Exact.</summary>
public static float3 af_neg(float3 x) => new(-x.x, -x.y, x.z);
/// <summary>
/// Encloses x * s: (x0 s, a s, e |s|). Rounding: S = |x0| + |a| + e of the result, plus 2^-126
/// for underflow.
/// </summary>
public static float3 af_scale(float3 x, float s)
{
var x0 = x.x * s;
var a = x.y * s;
var e = x.z * abs(s);
return new float3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7f + 1.17549435e-38f);
}
/// <summary>Encloses x + v for v in the interval, uncorrelated with eps: af_add(x, af_from_iv(interval)).</summary>
public static float3 af_add_iv(float3 x, float2 interval) => af_add(x, af_from_iv(interval));
/// <summary>
/// Encloses x * y, AF1's product: (x0 y0, x0 a_y + y0 a_x, |x0| e_y + |y0| e_x + (|a_x| + e_x)(|a_y| + e_y)).
/// The last term condenses the eps^2 term a_x a_y and every cross term with an error symbol into e,
/// so a form multiplied by itself is wider than its square: there is no af_sqr until a consumer
/// names one. Rounding: a can cancel, so S = |x0| + |x0 a_y| + |y0 a_x| + e, plus 2^-126 for
/// underflow.
/// </summary>
public static float3 af_mul(float3 x, float3 y)
{
var x0 = x.x * y.x;
var p = x.x * y.y;
var q = y.x * x.y;
var a = p + q;
var e = abs(x.x) * y.z + abs(y.x) * x.z + (abs(x.y) + x.z) * (abs(y.y) + y.z);
return new float3(x0, a, e + (abs(x0) + abs(p) + abs(q) + e) * 9.5367431640625e-7f + 1.17549435e-38f);
}
/// <summary>
/// Encloses snoise(p) for every p = centre + axis eps + w with |w| &lt;= radius, from one snoise_grad at
/// the centre (value n, gradient g). With R = |axis| + radius it returns the narrower of two sound
/// forms: the Lipschitz form (n, 0, L R), L = SNOISE_LIPSCHITZ, and the centred form
/// (n, g . axis, |g| radius + M R^2 / 2), M = SNOISE_HESSIAN. The centred form is Taylor's theorem:
/// snoise(c + d) = n + g . d + r with |r| &lt;= M |d|^2 / 2, d = axis eps + w, g . w &lt;= |g| radius.
/// radius must be at least 0.
///
/// Both forms carry an allowance for the float32 evaluation of snoise, which is not the exact
/// function: 2^-14 + 2^-12 R + L 2^-20 (|centre|_1 + R). The skew and the corner offsets round the
/// input point by at most about 2u |p|_1 before any kernel term is formed (the hashing is exact on
/// integers below 2^24), so a float32 snoise of p is the exact one at a point within 2u |p|_1, which
/// moves the value by at most L 2u |p|_1; for both p and the centre, with the caller's scaling of the
/// point by a frequency (one rounding) and a factor of two for a GPU, 16u = 2^-20. The kernel's
/// absolute rounding is at most 105 x four corners x terms of magnitude below 0.5^4 x 1.8 over about
/// 30 operations, under 2^-15, doubled to 2^-14; the gradient's error times |d| &lt;= R is under 2^-13 R,
/// doubled to 2^-12 R. AffineTests.SnoiseFormEnclosesBall checks the result with no tolerance on the
/// CPU, and the fixture's enclosure check measures it on a GPU. Rounding of the form itself: S =
/// |n| + |a| + e. The Lipschitz form is chosen when its half-width L R is at most the centred form's
/// |a| + |g| radius + M R^2 / 2, so the total width is never more than the Lipschitz form's.
/// </summary>
public static float3 af_snoise(float3 centre, float3 axis, float radius)
{
var g = snoise_grad(centre);
var gradient = new float3(g.x, g.y, g.z);
var reach = length(axis) + radius;
var lipschitz = SNOISE_LIPSCHITZ * reach;
var a = dot(gradient, axis);
var remainder = length(gradient) * radius + 0.5f * SNOISE_HESSIAN * reach * reach;
var allowance = 6.103515625e-5f + 2.44140625e-4f * reach
+ SNOISE_LIPSCHITZ * 9.5367431640625e-7f * (abs(centre.x) + abs(centre.y) + abs(centre.z) + reach);
var centred = abs(a) + remainder < lipschitz;
var slope = centred ? a : 0.0f;
var e = (centred ? remainder : lipschitz) + allowance;
return new float3(g.w, slope, e + (abs(g.w) + abs(slope) + e) * 9.5367431640625e-7f);
}
/// <summary>
/// Encloses the fBm value, sum of a_i snoise(p f_i) over the octaves, for every p = centre + axis eps
/// + w with |w| &lt;= radius: the octave sum of af_snoise(centre f_i, axis f_i, radius |f_i|) scaled by
/// a_i, with frequency f_i and amplitude a_i compounded exactly as fbm_grad compounds them (its .w is
/// what AffineTests.SnoiseFormEnclosesBall encloses). The rounding of p f_i, inside fbm_grad, and of
/// the scaled centre and axis here, moves each point by at most 3u (|centre f_i|_1 + R f_i), inside
/// af_snoise's allowance for its input point. octaves is clamped to [0, 16] as there.
/// </summary>
public static float3 af_fbm(float3 centre, float3 axis, float radius, int octaves, float lacunarity, float gain)
{
octaves = clamp(octaves, 0, 16);
var amplitude = 1.0f;
var frequency = 1.0f;
var sum = new float3(0.0f, 0.0f, 0.0f);
for (var i = 0; i < octaves; i++)
{
sum = af_add(sum, af_scale(af_snoise(centre * frequency, axis * frequency, radius * abs(frequency)), amplitude));
frequency *= lacunarity;
amplitude *= gain;
}
return sum;
}
/// <summary>
/// The axis of a screen tile's slice over the depth segment [z0, z1], for af_frustum_ball's centre:
/// (m_c, 1) h with h = (z1 - z0) / 2, in the camera frame of iv_frustum_ball (a pinhole at the
/// origin looking down +z, the ray of slope m the points (m z, z)). Its z component is h itself, the
/// symbol's coefficient for depth: z = z_m + h eps. Its rounding is af_frustum_ball's.
/// </summary>
public static float3 af_frustum_axis(float2 centreSlope, float z0, float z1)
{
var h = (z1 - z0) * 0.5f;
return new float3(centreSlope.x * h, centreSlope.y * h, h);
}
/// <summary>
/// The affine slice of a screen tile over [z0, z1]: a centre (xyz) and radius (w) such that every point
/// of every ray of the tile at depth z in [z0, z1] is centre + axis eps(z) + w', |w'| &lt;= radius, with
/// axis = af_frustum_axis(centreSlope, z0, z1), z_m = (z0 + z1) / 2, h = (z1 - z0) / 2 (both float32,
/// as computed here) and eps(z) = clamp((z - z_m) / h, -1, 1) (0 where h is 0). The centre is
/// (m_c z_m, z_m) and the radius z1 footprintPerDepth + warpVariation, where footprintPerDepth is
/// iv_frustum_ball's and warpVariation is what the caller's warp can add once its value at the centre
/// is applied to the centre (design.md, "Affine forms", the centred warp). Derivation: a point of
/// slope m_c + e (|e| &lt;= footprintPerDepth) at depth z is (m_c z_m, z_m) + (m_c, 1)(z - z_m) +
/// (e z, 0), and (m_c, 1)(z - z_m) = axis eps(z) + (m_c, 1) r(z), where the clamp residual
/// |r(z)| = |z - z_m - h eps(z)| is at most u (|z_m| + h), one rounding each of z_m and h. A caller
/// taking depth as af_symbol(z_m, h) adds the same residual, 2^-23 (|z_m| + h) covers it, to that
/// form's e. Rounding: the centre's products and the axis's products round once each (u |c|_1 and
/// u |axis|_1 against |eps| &lt;= 1), the residual is at most u (|c|_1 + |axis|_1) since
/// |(m_c, 1)| &lt;= |m_c|_1 + 1, and the radius rounds at most four times; 2^-20 (|c|_1 + |axis|_1 +
/// radius) covers it twice over (AffineTests.FrustumSliceEnclosesEveryRayPoint, extreme families
/// included). The caller rotates and translates the centre into world space, rotates the axis, and
/// scales all three by its noise frequency; the rounding of that transform is the caller's.
/// </summary>
public static float4 af_frustum_ball(float2 centreSlope, float z0, float z1, float footprintPerDepth, float warpVariation)
{
var zm = (z0 + z1) * 0.5f;
var h = (z1 - z0) * 0.5f;
var cx = centreSlope.x * zm;
var cy = centreSlope.y * zm;
var spread = (abs(centreSlope.x) + abs(centreSlope.y) + 1.0f) * h;
var radius = z1 * footprintPerDepth + warpVariation;
radius += (abs(cx) + abs(cy) + abs(zm) + spread + radius) * 9.5367431640625e-7f;
return new float4(cx, cy, zm, radius);
}
}

2
packages/cultmath/src/CultMath/math.Affine.cs.meta Normal file
View file

@ -0,0 +1,2 @@
fileFormatVersion: 2
guid: 86b87ffb6c8a4d72b067a64f9f8a06ce

21
packages/cultmath/src/CultMath/math.Interval.cs
View file

@ -20,6 +20,14 @@ public static partial class math
// snoise kernel must re-run MeasureLipschitz and re-pin this constant.
public const float SNOISE_LIPSCHITZ = 10.099261f;
// Lipschitz constant of snoise(float2) in its input, with SNOISE_LIPSCHITZ's provenance rule:
// NoiseBoundTests.MeasureLipschitz2 (slow, explicit) takes the largest gradient length (central
// differences over exact float steps) over 1e6 seeded points, refines the 1e4 largest by ascent and
// multiplies by 1.10; NoiseBoundTests.Lipschitz2ConstantPins pins it against a fast re-measure and
// the ascent from MeasureLipschitz2's witness. A change to the 2D snoise kernel must re-run
// MeasureLipschitz2 and re-pin this constant.
public const float SNOISE2_LIPSCHITZ = 9.19f;
/// <summary>The point interval [x, x].</summary>
public static float2 iv_point(float value) => new(value, value);
@ -142,6 +150,19 @@ public static partial class math
return new(max(n - e, -1.0f), min(n + e, 1.0f));
}
/// <summary>
/// Encloses snoise(x) for every 2D x within distance radius of centre: [n - L2 r, n + L2 r] intersected
/// with [-1, 1], n = snoise(centre), L2 = SNOISE2_LIPSCHITZ. radius must be at least 0. One snoise
/// evaluation, value only. A consumer's 2D flow built from snoise(float2) varies over a ball by at
/// most the width this gives (NoiseBoundTests.CentredWarpStaysEnclosed).
/// </summary>
public static float2 iv_snoise_ball(float2 centre, float radius)
{
var n = snoise(centre);
var e = SNOISE2_LIPSCHITZ * radius;
return new(max(n - e, -1.0f), min(n + e, 1.0f));
}
/// <summary>
/// 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,

458
packages/cultmath/tests/CultMath.Tests/AffineTests.cs Normal file
View file

@ -0,0 +1,458 @@
using CultMath;
using Xunit;
using static CultMath.math;
namespace CultMath.Tests;
/// <summary>
/// The affine forms of math.Affine.cs (design.md, "Affine forms"). A form (x0, a, e) is the set
/// x0 + a eps + e delta at the region's shared eps; every check evaluates that set in double and the
/// pointwise function in float32, with no tolerance. Operand values are float32 values drawn inside
/// their operands' sets at one shared eps, endpoints (eps and delta at -1 and 1) included.
/// </summary>
public sealed class AffineTests
{
private readonly ITestOutputHelper output;
public AffineTests(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 float Signed(System.Random random, float magnitude) => random.Next(2) == 0 ? -magnitude : magnitude;
// Whether value lies in the form's set at eps, in double.
private static bool InForm(float3 form, float eps, double value) =>
Math.Abs(value - (form.x + (double)form.y * eps)) <= form.z;
// A form with its centre in [-20, 20] and |a| and e log-uniform over [1e-5, 20]; one in eight has
// a = 0 and one in eight e = 0.
private static float3 DrawForm(System.Random random)
{
var x0 = Uniform(random, -20.0f, 20.0f);
var a = random.Next(8) == 0 ? 0.0f : Signed(random, LogUniform(random, 1.0e-5f, 20.0f));
var e = random.Next(8) == 0 ? 0.0f : LogUniform(random, 1.0e-5f, 20.0f);
return new float3(x0, a, e);
}
private static float2 DrawInterval(System.Random random)
{
var lo = Uniform(random, -20.0f, 20.0f);
return new float2(lo, lo + (random.Next(8) == 0 ? 0.0f : LogUniform(random, 1.0e-5f, 40.0f)));
}
// A float32 value in the form's set at eps, delta at an endpoint when asked; NaN when rounding leaves
// no float of the set near the drawn point (a set narrower than an ulp).
private static float Member(System.Random random, float3 form, float eps, bool endpoint)
{
var centre = form.x + (double)form.y * eps;
var delta = endpoint ? (random.Next(2) == 0 ? -1.0 : 1.0) : Uniform(random, -1.0f, 1.0f);
var value = (float)(centre + form.z * delta);
for (var k = 0; k < 4 && !InForm(form, eps, value); k++)
value = value > centre ? MathF.BitDecrement(value) : MathF.BitIncrement(value);
return InForm(form, eps, value) ? value : float.NaN;
}
private static float MemberOf(System.Random random, float2 interval, bool endpoint) =>
endpoint ? (random.Next(2) == 0 ? interval.x : interval.y) : Math.Clamp(Uniform(random, interval.x, interval.y), interval.x, interval.y);
private readonly record struct Operands(float3 X, float3 Y, float2 Iv, float S);
private readonly record struct Values(float X, float Y, float V, float S);
// One op under the harness: Build gives the result's membership test (eps, value) from the operand
// forms, F the pointwise float32 function of operand values.
private sealed record Op(string Name, Func<Operands, Func<float, double, bool>> Build, Func<Values, float> F);
private static Func<float, double, bool> Form(float3 form) => (eps, value) => InForm(form, eps, value);
private static Func<float, double, bool> Range(float2 range) => (_, value) => value >= range.x && value <= range.y;
private static readonly Op[] Ops =
{
new("af_point", o => Form(af_point(o.X.x)), v => v.S),
new("af_symbol", o => Form(af_symbol(o.X.x, o.X.y)), v => v.S),
new("af_from_iv", o => Form(af_from_iv(o.Iv)), v => v.V),
new("af_range", o => Range(af_range(o.X)), v => v.X),
new("af_add", o => Form(af_add(o.X, o.Y)), v => v.X + v.Y),
new("af_sub", o => Form(af_sub(o.X, o.Y)), v => v.X - v.Y),
new("af_neg", o => Form(af_neg(o.X)), v => -v.X),
new("af_scale", o => Form(af_scale(o.X, o.S)), v => v.X * v.S),
new("af_add_iv", o => Form(af_add_iv(o.X, o.Iv)), v => v.X + v.V),
new("af_mul", o => Form(af_mul(o.X, o.Y)), v => v.X * v.Y),
};
// Runs one op over `forms` seeded operand draws x 16 (eps, delta) draws; returns the misses and up to
// four witnesses. The first four draws per form put eps at an endpoint, the first eight delta.
// af_point and af_symbol are checked on their defining value: the point x0 itself, and x0 + a eps in
// double (the symbol defines eps, so it has no float32 evaluation to enclose).
private static (int Misses, List<string> Witness) Harness(Op op, int seed, int forms = 10_000, Func<System.Random, float3>? draw = null)
{
draw ??= DrawForm;
var random = new System.Random(seed);
var misses = 0;
var witness = new List<string>();
for (var t = 0; t < forms; t++)
{
var operands = new Operands(draw(random), draw(random), DrawInterval(random), random.Next(8) == 0 ? 0.0f : Uniform(random, -4.0f, 4.0f));
var contains = op.Build(operands);
for (var d = 0; d < 16; d++)
{
var eps = d < 4 ? (d % 2 == 0 ? -1.0f : 1.0f) : Uniform(random, -1.0f, 1.0f);
var endpoint = d < 8;
var x = Member(random, operands.X, eps, endpoint);
var y = Member(random, operands.Y, eps, endpoint);
if (float.IsNaN(x) || float.IsNaN(y))
continue;
var v = MemberOf(random, operands.Iv, endpoint);
double f = op.Name switch
{
"af_point" => operands.X.x,
"af_symbol" => operands.X.x + (double)operands.X.y * eps,
_ => op.F(new Values(x, y, v, operands.S)),
};
if (!contains(eps, f))
{
misses++;
if (witness.Count < 4)
witness.Add($"{op.Name}: X {operands.X}, Y {operands.Y}, iv {operands.Iv}, s {operands.S:R}; eps {eps:R}, values ({x:R}, {y:R}, {v:R}) give {f:R}");
}
}
}
return (misses, witness);
}
[Fact]
public void EveryOpEnclosesItsPoints()
{
var failures = new List<string>();
foreach (var op in Ops)
{
var (misses, witness) = Harness(op, 0xAF00);
if (misses > 0)
failures.Add($"{op.Name}: {misses} misses" + Environment.NewLine + string.Join(Environment.NewLine, witness));
}
Assert.True(failures.Count == 0, string.Join(Environment.NewLine, failures));
}
private static float3 Shrink(float3 form) => new(form.x, form.y, form.z > 0.0f ? MathF.BitDecrement(form.z) : 0.0f);
// Each op with its rounding absorption removed: the propagated error alone, computed in float32.
private static readonly Op[] Unabsorbed =
{
new("af_from_iv", o => Form(new float3((o.Iv.x + o.Iv.y) * 0.5f, 0.0f, (o.Iv.y - o.Iv.x) * 0.5f)), v => v.V),
new("af_range", o => Range(new float2(o.X.x - (MathF.Abs(o.X.y) + o.X.z), o.X.x + (MathF.Abs(o.X.y) + o.X.z))), v => v.X),
new("af_add", o => Form(new float3(o.X.x + o.Y.x, o.X.y + o.Y.y, o.X.z + o.Y.z)), v => v.X + v.Y),
new("af_sub", o => Form(new float3(o.X.x - o.Y.x, o.X.y - o.Y.y, o.X.z + o.Y.z)), v => v.X - v.Y),
new("af_scale", o => Form(new float3(o.X.x * o.S, o.X.y * o.S, o.X.z * MathF.Abs(o.S))), v => v.X * v.S),
new("af_add_iv", o => Form(new float3(o.X.x + (o.Iv.x + o.Iv.y) * 0.5f, o.X.y, o.X.z + (o.Iv.y - o.Iv.x) * 0.5f)), v => v.X + v.V),
new("af_mul", o => Form(new float3(o.X.x * o.Y.x, o.X.x * o.Y.y + o.Y.x * o.X.y,
MathF.Abs(o.X.x) * o.Y.z + MathF.Abs(o.Y.x) * o.X.z + (MathF.Abs(o.X.y) + o.X.z) * (MathF.Abs(o.Y.y) + o.Y.z))), v => v.X * v.Y),
};
// One-ulp mutants that survive, each with the reason. A rounding absorption bounds the worst case
// over every operand value, and 2^-20 S is 16u S, at least eight ulps of the larger of |x0|, |a| and
// e; one ulp less of e leaves the rest of that margin, which no sampled rounding can cross.
private static readonly Dictionary<string, string> OneUlpSurvivors = new()
{
["af_from_iv"] = "the absorption (|x0| + e) 2^-20 is at least eight ulps of e over the one rounding of x0 and e it covers",
["af_range"] = "each endpoint moves outward by (|x0| + |a| + e) 2^-20, eight ulps past the one rounding it covers",
["af_add"] = "the absorption 2^-20 S is eight ulps of S's largest term past the three roundings it covers",
["af_sub"] = "as af_add",
["af_scale"] = "as af_add",
["af_add_iv"] = "two absorptions, af_from_iv's and af_add's",
["af_mul"] = "the absorption 2^-20 (|x0| + |x0 a_y| + |y0 a_x| + e) is eight ulps past the products it covers",
};
/// <summary>
/// The rounding absorptions are load-bearing: with each op's absorption removed (Unabsorbed), the
/// harness finds a miss for every op. Each op's e shrunk by one ulp is run too; it survives for every
/// op, because the absorption is a worst-case bound many ulps wide, and each survivor is named with
/// its reason in OneUlpSurvivors. A one-ulp mutant that starts failing, or one that is not named,
/// fails this test, so the list cannot go stale.
/// </summary>
[Fact]
public void EveryRoundingTermIsCaught()
{
var problems = new List<string>();
foreach (var mutant in Unabsorbed)
{
var (misses, _) = Harness(mutant, 0xAF01);
output.WriteLine($"AF-REPORT {mutant.Name} without absorption: {misses} misses");
if (misses == 0)
problems.Add($"{mutant.Name}: removing the rounding absorption is not caught");
}
foreach (var op in Ops.Where(o => o.Name is not ("af_point" or "af_symbol" or "af_neg")))
{
var shrunk = op.Name == "af_range"
? op with { Build = o => { var r = af_range(o.X); return Range(new float2(MathF.BitIncrement(r.x), MathF.BitDecrement(r.y))); } }
: op with { Build = o => Form(Shrink(BuildForm(op.Name, o))) };
var (misses, _) = Harness(shrunk, 0xAF02);
output.WriteLine($"AF-REPORT {op.Name} with e one ulp short: {misses} misses");
var named = OneUlpSurvivors.ContainsKey(op.Name);
if (misses == 0 && !named)
problems.Add($"{op.Name}: e one ulp short survives and is not named");
if (misses > 0 && named)
problems.Add($"{op.Name}: e one ulp short is now caught; remove it from OneUlpSurvivors");
}
Assert.True(problems.Count == 0, string.Join(Environment.NewLine, problems));
}
private static float3 BuildForm(string name, Operands o) => name switch
{
"af_from_iv" => af_from_iv(o.Iv),
"af_add" => af_add(o.X, o.Y),
"af_sub" => af_sub(o.X, o.Y),
"af_scale" => af_scale(o.X, o.S),
"af_add_iv" => af_add_iv(o.X, o.Iv),
"af_mul" => af_mul(o.X, o.Y),
_ => throw new ArgumentException(name),
};
/// <summary>
/// The shared symbol correlates: for x = af_symbol(x0, a), af_sub(x, x) and af_add(x, af_neg(x)) have
/// a = 0 and e at most their own rounding, and af_sub(af_scale(x, 2), af_add(x, x)) has a = 0. A
/// mutant of af_add, af_sub or af_scale that condensed a into e would leave e of order |a|.
/// </summary>
[Fact]
public void SharedSymbolCorrelates()
{
var random = new System.Random(0xAF03);
for (var t = 0; t < 10_000; t++)
{
var x = af_symbol(Uniform(random, -20.0f, 20.0f), Signed(random, LogUniform(random, 1.0e-5f, 20.0f)));
foreach (var (name, form) in new[]
{
("af_sub(x, x)", af_sub(x, x)),
("af_add(x, af_neg(x))", af_add(x, af_neg(x))),
("af_sub(af_scale(x, 2), af_add(x, x))", af_sub(af_scale(x, 2.0f), af_add(x, x))),
})
{
var rounding = (MathF.Abs(x.x) + MathF.Abs(x.y)) * 1.0e-5f;
Assert.True(form.y == 0.0f && form.z <= rounding, $"{name} for x = {x} is {form}: a must be 0 and e at most {rounding:R}");
}
}
}
// Forms whose range straddles zero (|x0| below |a| + e) or is one-signed (|x0| above it).
private static float3 DrawSigned(System.Random random, bool straddle)
{
var a = Signed(random, LogUniform(random, 1.0e-3f, 10.0f));
var e = LogUniform(random, 1.0e-3f, 10.0f);
var reach = MathF.Abs(a) + e;
return new float3(Signed(random, straddle ? Uniform(random, 0.0f, 0.9f * reach) : reach * Uniform(random, 1.1f, 4.0f)), a, e);
}
/// <summary>
/// af_mul encloses the product over forms that straddle zero and forms that are one-signed, every
/// sign combination. The eps^2 term a_x a_y is what the condensed e carries for eps at +-1: a mutant
/// dropping |a_x a_y| from e misses there.
/// </summary>
[Fact]
public void MulEnclosesEverySignCase()
{
var mul = Ops.Single(o => o.Name == "af_mul");
var dropped = mul with
{
Build = o =>
{
var r = af_mul(o.X, o.Y);
return Form(new float3(r.x, r.y, r.z - MathF.Abs(o.X.y * o.Y.y)));
},
};
foreach (var (sx, sy) in new[] { (true, true), (true, false), (false, true), (false, false) })
{
var seed = 0xAF04 + (sx ? 1 : 0) + (sy ? 2 : 0);
var flip = 0;
float3 Draw(System.Random r) => DrawSigned(r, (flip++ % 2 == 0) ? sx : sy);
var (misses, witness) = Harness(mul, seed, 5_000, Draw);
Assert.True(misses == 0, $"af_mul, straddling ({sx}, {sy}): {misses} misses" + Environment.NewLine + string.Join(Environment.NewLine, witness));
flip = 0;
var (caught, _) = Harness(dropped, seed, 5_000, Draw);
Assert.True(caught > 0, $"af_mul without |a_x a_y|, straddling ({sx}, {sy}), is not caught");
}
}
// af_snoise's float32 allowance, as its doc comment states it, for the width comparison below.
private static double Allowance(float3 centre, double reach) =>
Math.ScaleB(1.0, -14) + Math.ScaleB(reach, -12)
+ SNOISE_LIPSCHITZ * Math.ScaleB(Math.Abs(centre.x) + Math.Abs(centre.y) + Math.Abs(centre.z) + reach, -20);
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);
}
// A float32 point of the set centre + axis eps + w, |w| <= radius, checked in double; the first eight
// w are on the surface (pulled in by 2^-20 so their rounding stays inside). NaN-free: a point whose
// rounding leaves the set is reported as null and skipped.
private static float3? PointOf(System.Random random, float3 centre, float3 axis, float radius, float eps, int index)
{
var w = UnitVector(random) * (index < 8 ? radius * (1.0f - 9.5367431640625e-7f) : radius * MathF.Cbrt(random.NextSingle()));
var p = centre + axis * eps + w;
var (dx, dy, dz) = (p.x - (centre.x + (double)axis.x * eps), p.y - (centre.y + (double)axis.y * eps), p.z - (centre.z + (double)axis.z * eps));
return Math.Sqrt(dx * dx + dy * dy + dz * dz) <= radius ? p : null;
}
/// <summary>
/// af_snoise and af_fbm over 2,000 (centre, axis, radius) each: centres in [-20, 20]^3, radius
/// log-uniform in [1e-4, 2], |axis| uniform in [0, 2] (an eighth exactly 0). 64 points
/// centre + axis eps + w per case, eps at -1 and 1 for the first eight; snoise (fbm_grad(...).w) in
/// float32 lies in the form at eps, no tolerance. Dropping the M R^2 / 2 remainder from a centred
/// form misses. af_snoise's total width is never more than the Lipschitz form's, L R plus the
/// allowance.
/// </summary>
[Fact]
public void SnoiseFormEnclosesBall()
{
var random = new System.Random(0xAF10);
var (misses, remainderMisses, skipped, centredCases) = (0, 0, 0, 0);
var witness = new List<string>();
for (var b = 0; b < 4000; b++)
{
var fbm = b >= 2000;
var centre = new float3(Uniform(random, -20.0f, 20.0f), Uniform(random, -20.0f, 20.0f), Uniform(random, -20.0f, 20.0f));
var axis = random.Next(8) == 0 ? new float3(0.0f, 0.0f, 0.0f) : UnitVector(random) * Uniform(random, 0.0f, 2.0f);
var radius = LogUniform(random, 1.0e-4f, 2.0f);
var (octaves, lacunarity, gain) = (random.Next(0, 7), Uniform(random, 1.5f, 2.5f), Uniform(random, 0.3f, 0.7f));
var form = fbm ? af_fbm(centre, axis, radius, octaves, lacunarity, gain) : af_snoise(centre, axis, radius);
var reach = (double)length(axis) + radius;
var remainder = 0.5 * SNOISE_HESSIAN * reach * reach;
var mutant = new float3(form.x, form.y, (float)(form.z - remainder));
if (!fbm)
{
var lipschitz = (SNOISE_LIPSCHITZ * reach + Allowance(centre, reach) + Math.ScaleB(Math.Abs(form.x), -20)) * (1.0 + Math.ScaleB(1.0, -18));
Assert.True(Math.Abs(form.y) + form.z <= lipschitz, $"af_snoise({centre}, {axis}, {radius:R}) = {form} is wider than the Lipschitz form's half-width {lipschitz:R}");
centredCases += form.y != 0.0f ? 1 : 0;
}
for (var i = 0; i < 64; i++)
{
var eps = i < 8 ? (i % 2 == 0 ? -1.0f : 1.0f) : Uniform(random, -1.0f, 1.0f);
var p = PointOf(random, centre, axis, radius, eps, i);
if (p is not { } point)
{
skipped++;
continue;
}
var value = fbm ? fbm_grad(point, octaves, lacunarity, gain).w : snoise(point);
if (!InForm(form, eps, value))
{
misses++;
if (witness.Count < 4)
witness.Add($"{(fbm ? $"af_fbm(..., {octaves}, {lacunarity:R}, {gain:R})" : "af_snoise")}({centre}, {axis}, {radius:R}) = {form} misses {value:R} at eps {eps:R}, {point}");
}
if (!fbm && form.y != 0.0f && !InForm(mutant, eps, value))
remainderMisses++;
}
}
output.WriteLine($"AF-REPORT af_snoise/af_fbm enclosure: {misses} misses, {skipped} points skipped (rounded out of the set), {centredCases} of 2000 af_snoise cases centred, {remainderMisses} misses with the remainder dropped");
Assert.True(misses == 0, string.Join(Environment.NewLine, witness));
Assert.True(skipped < 64 * 4000 / 100, $"{skipped} points rounded out of their sets");
Assert.True(remainderMisses > 0, "dropping M R^2 / 2 from the centred form is not caught");
}
/// <summary>
/// At radius 1e-3 and |axis| = 1e-3, af_snoise takes the centred form in every one of 2,000 seeded
/// cases (a is nonzero), and its mean total width is under half of the Lipschitz form's 2 L R. A
/// mutant that always returns the Lipschitz form has a = 0.
/// </summary>
[Fact]
public void SnoiseChoosesTheCentredFormWhereNarrower()
{
var random = new System.Random(0xAF11);
var sum = 0.0;
for (var b = 0; b < 2000; b++)
{
var centre = new float3(Uniform(random, -20.0f, 20.0f), Uniform(random, -20.0f, 20.0f), Uniform(random, -20.0f, 20.0f));
var axis = UnitVector(random) * 1.0e-3f;
var form = af_snoise(centre, axis, 1.0e-3f);
var lipschitz = 2.0 * SNOISE_LIPSCHITZ * 2.0e-3;
var width = 2.0 * (Math.Abs(form.y) + form.z);
Assert.True(form.y != 0.0f && width < lipschitz, $"af_snoise({centre}, {axis}, 1e-3) = {form}: not the centred form, or no narrower than 2 L R = {lipschitz:R}");
sum += width / lipschitz;
}
output.WriteLine($"AF-REPORT af_snoise at radius 1e-3: mean width over 2 L R {sum / 2000:F3}");
Assert.True(sum / 2000 < 0.5, $"mean width over 2 L R is {sum / 2000:F3}, not under half");
}
/// <summary>
/// af_frustum_ball and af_frustum_axis over iv_frustum_ball's families (NoiseBoundTests.DrawBallCase:
/// the r2 domain and the degenerate, thin, wide-slope, far, from-the-camera, tiny-focal and huge-warp
/// extremes), 2,000 tiles each, warp passed as the warp variation: every exact ray point (m z, z),
/// formed in double from the float32 slope and depth at the footprint's corners at both depth ends
/// (and 56 more for the r2 domain), moved by a full-length outward flow of length warp, is within
/// radius of centre + axis eps(z), eps(z) = clamp((z - z_m) / h, -1, 1) with z_m and h the centre's z
/// and the axis's z (0 when h is 0). The radius is pinned to the derivation both ways:
/// r &lt;= w &lt;= r + 2^-19 (r + |c|_1 + |axis|_1), r = z1 footprintPerDepth + warp in double.
/// </summary>
[Fact]
public void FrustumSliceEnclosesEveryRayPoint()
{
var random = new System.Random(0xAF12);
string[] families = { "r2 domain", "degenerate", "thin", "wide slope", "far", "from the camera", "tiny focal", "huge warp" };
var failures = new Dictionary<string, int>();
var witness = new List<string>();
foreach (var family in families)
{
failures[family] = 0;
for (var t = 0; t < 2000; t++)
{
var (n, focal, slope, z0, z1, warp) = NoiseBoundTests.DrawBallCase(random, family);
var half = n / (2.0f * focal);
var footprint = n / (MathF.Sqrt(2.0f) * focal);
var axis = af_frustum_axis(slope, z0, z1);
var ball = af_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 r = z1 * (double)footprint + warp;
var c1 = Math.Abs(cx) + Math.Abs(cy) + Math.Abs(cz);
var a1 = Math.Abs((double)axis.x) + Math.Abs((double)axis.y) + Math.Abs((double)axis.z);
if (!(r <= w && w <= r + Math.ScaleB(r + c1 + a1, -19)))
{
failures[family]++;
if (witness.Count < 8)
witness.Add($"{family} pin: af_frustum_ball({slope}, {z0:R}, {z1:R}, {footprint:R}, {warp:R}) = {ball}, derivation {r: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 eps = axis.z == 0.0f ? 0.0 : Math.Clamp((z - cz) / axis.z, -1.0, 1.0);
var (qx, qy, qz) = (cx + axis.x * eps, cy + axis.y * eps, cz + axis.z * eps);
var (dx, dy, dz) = (m.x * (double)z - qx, m.y * (double)z - qy, z - qz);
var away = Math.Sqrt(dx * dx + dy * dy + dz * dz);
var distance = away + warp;
if (!(distance <= w))
{
failures[family]++;
if (witness.Count < 8)
witness.Add($"{family}: af_frustum_ball({slope}, {z0:R}, {z1:R}, {footprint:R}, {warp:R}) = {ball}, axis {axis} misses slope {m}, depth {z:R} by {(distance - w) / w:R} r");
}
}
}
}
Assert.True(failures.Values.Sum() == 0, string.Join(", ", families.Select(f => $"{f} {failures[f]}")) + Environment.NewLine + string.Join(Environment.NewLine, witness));
}
}

2
packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs
View file

@ -121,7 +121,7 @@ public sealed class NoiseBoundTests
// 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)
internal 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 };