diff --git a/packages/cultmath/src/CultMath/math.Affine.cs b/packages/cultmath/src/CultMath/math.Affine.cs
new file mode 100644
index 00000000..1d2dd70b
--- /dev/null
+++ b/packages/cultmath/src/CultMath/math.Affine.cs
@@ -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;
+
+ /// The point form (x, 0, 0): the value x at every eps. Exact.
+ public static float3 af_point(float value) => new(value, 0.0f, 0.0f);
+
+ ///
+ /// 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.
+ ///
+ public static float3 af_symbol(float x0, float a) => new(x0, a, 0.0f);
+
+ ///
+ /// 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.
+ ///
+ 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);
+ }
+
+ ///
+ /// 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.
+ ///
+ 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);
+ }
+
+ ///
+ /// 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.
+ ///
+ 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);
+ }
+
+ /// Encloses x - y: (x0 - y0, a_x - a_y, e_x + e_y). Rounding as af_add.
+ 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);
+ }
+
+ /// Encloses -x: (-x0, -a, e). Exact.
+ public static float3 af_neg(float3 x) => new(-x.x, -x.y, x.z);
+
+ ///
+ /// Encloses x * s: (x0 s, a s, e |s|). Rounding: S = |x0| + |a| + e of the result, plus 2^-126
+ /// for underflow.
+ ///
+ 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);
+ }
+
+ /// Encloses x + v for v in the interval, uncorrelated with eps: af_add(x, af_from_iv(interval)).
+ public static float3 af_add_iv(float3 x, float2 interval) => af_add(x, af_from_iv(interval));
+
+ ///
+ /// 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.
+ ///
+ 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);
+ }
+
+ ///
+ /// Encloses snoise(p) for every p = centre + axis eps + w with |w| <= 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| <= M |d|^2 / 2, d = axis eps + w, g . w <= |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| <= 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.
+ ///
+ 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);
+ }
+
+ ///
+ /// Encloses the fBm value, sum of a_i snoise(p f_i) over the octaves, for every p = centre + axis eps
+ /// + w with |w| <= 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.
+ ///
+ 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;
+ }
+
+ ///
+ /// 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.
+ ///
+ 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);
+ }
+
+ ///
+ /// 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'| <= 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| <= 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| <= 1), the residual is at most u (|c|_1 + |axis|_1) since
+ /// |(m_c, 1)| <= |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.
+ ///
+ 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);
+ }
+}
diff --git a/packages/cultmath/src/CultMath/math.Affine.cs.meta b/packages/cultmath/src/CultMath/math.Affine.cs.meta
new file mode 100644
index 00000000..34ba0047
--- /dev/null
+++ b/packages/cultmath/src/CultMath/math.Affine.cs.meta
@@ -0,0 +1,2 @@
+fileFormatVersion: 2
+guid: 86b87ffb6c8a4d72b067a64f9f8a06ce
diff --git a/packages/cultmath/src/CultMath/math.Interval.cs b/packages/cultmath/src/CultMath/math.Interval.cs
index 91e99f87..deedea32 100644
--- a/packages/cultmath/src/CultMath/math.Interval.cs
+++ b/packages/cultmath/src/CultMath/math.Interval.cs
@@ -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;
+
/// The point interval [x, x].
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));
}
+ ///
+ /// 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).
+ ///
+ 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));
+ }
+
///
/// 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,
diff --git a/packages/cultmath/tests/CultMath.Tests/AffineTests.cs b/packages/cultmath/tests/CultMath.Tests/AffineTests.cs
new file mode 100644
index 00000000..38ccd071
--- /dev/null
+++ b/packages/cultmath/tests/CultMath.Tests/AffineTests.cs
@@ -0,0 +1,458 @@
+using CultMath;
+using Xunit;
+using static CultMath.math;
+
+namespace CultMath.Tests;
+
+///
+/// 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.
+///
+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> Build, Func F);
+
+ private static Func Form(float3 form) => (eps, value) => InForm(form, eps, value);
+
+ private static Func 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 Witness) Harness(Op op, int seed, int forms = 10_000, Func? draw = null)
+ {
+ draw ??= DrawForm;
+ var random = new System.Random(seed);
+ var misses = 0;
+ var witness = new List();
+ 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();
+ 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 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",
+ };
+
+ ///
+ /// 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.
+ ///
+ [Fact]
+ public void EveryRoundingTermIsCaught()
+ {
+ var problems = new List();
+ 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),
+ };
+
+ ///
+ /// 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|.
+ ///
+ [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);
+ }
+
+ ///
+ /// 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.
+ ///
+ [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;
+ }
+
+ ///
+ /// 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.
+ ///
+ [Fact]
+ public void SnoiseFormEnclosesBall()
+ {
+ var random = new System.Random(0xAF10);
+ var (misses, remainderMisses, skipped, centredCases) = (0, 0, 0, 0);
+ var witness = new List();
+ 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");
+ }
+
+ ///
+ /// 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.
+ ///
+ [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");
+ }
+
+ ///
+ /// 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 <= w <= r + 2^-19 (r + |c|_1 + |axis|_1), r = z1 footprintPerDepth + warp in double.
+ ///
+ [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();
+ var witness = new List();
+ 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));
+ }
+}
diff --git a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs
index 2b00ca88..fc0ba0ac 100644
--- a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs
+++ b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs
@@ -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 };