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

CultMath intervals: enclosure harness, one-ulp shrink mutants, noise bounds, Lipschitz pin, saving

Browse source 
IntervalTests runs every op over 10,000 seeded intervals x 16 points with no tolerance, and
kills each op's one-ulp-shrunk lo and hi; exp's widening is an equivalent mutant on this CPU,
and smoothstep's is killed by a constructed non-monotone step. NoiseBoundTests encloses snoise
and fbm over 2,000 balls x 64 points and the bounded domain warp, pins SNOISE_LIPSCHITZ against
1e5 points and the ascent from MeasureLipschitz's witness, prints tightness, and prints the
march saving. The saving is 0.93x, not 2x: that test skips with the fork named.
MeasureLipschitz and MeasureMonotonicity are explicit, slow tests.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Meta Vens 2026-10-02 19:56:30 +02:00
commit cd7913bfcc
2 changed files with 739 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

312
packages/cultmath/tests/CultMath.Tests/IntervalTests.cs Normal file
View file

@ -0,0 +1,312 @@
using CultMath;
using Xunit;
using static CultMath.math;
namespace CultMath.Tests;
/// <summary>
/// intervals-enclose for every iv_* op (design.md, "Intervals"): for each op, intervals with seeded random
/// centres and widths and the points inside them; f(x), evaluated pointwise in float32 by the same C# math
/// the op mirrors, must lie in op(I) with no tolerance.
///
/// Stryker mutates operators, not float bounds, so the mutants that matter most here are written by hand
/// and run by <see cref="EveryOneUlpShrinkIsCaught"/>: for every op, the lo bound raised by one ulp and the
/// hi bound lowered by one ulp must each fail the harness. Two pairs are equivalent mutants on this CPU and
/// are not in that list:
/// - iv_exp's one-ulp widening, both sides: MathF.Exp is monotone over every float on the Linux runtime
/// (MeasureMonotonicity prints 0 violations), so removing the widening encloses here. It is kept because
/// exp is not required to be correctly rounded, and another libm or a GPU need not be monotone.
/// - iv_smoothstep's one-ulp widening: the generic harness never meets one of smoothstep's 1-ulp
/// non-monotone steps; <see cref="SmoothstepWideningCoversItsFloatNonMonotonicity"/> constructs one on
/// each side and kills both mutants.
/// </summary>
public sealed class IntervalTests
{
private readonly ITestOutputHelper output;
public IntervalTests(ITestOutputHelper output) => this.output = output;
private const int IntervalCount = 10000;
private const int PointsPerInterval = 16;
// One row of the shared harness: how many interval arguments the op takes, how to draw its scalar
// arguments, the op itself and the pointwise function it encloses. A new op is one entry here.
internal sealed record Op(
string Name,
int Intervals,
Func<System.Random, float[]> Scalars,
Func<float2[], float[], float2> Bound,
Func<float[], float[], float> Pointwise,
bool NonNegativeDomain = false);
private static float[] None(System.Random random) => Array.Empty<float>();
private static float Uniform(System.Random random, float lo, float hi) => lo + (hi - lo) * random.NextSingle();
internal static readonly Op[] Ops =
{
new("iv_point", 0, r => new[] { Uniform(r, -20.0f, 20.0f) }, (i, s) => iv_point(s[0]), (x, s) => s[0]),
new("iv_add", 2, None, (i, s) => iv_add(i[0], i[1]), (x, s) => x[0] + x[1]),
new("iv_sub", 2, None, (i, s) => iv_sub(i[0], i[1]), (x, s) => x[0] - x[1]),
new("iv_neg", 1, None, (i, s) => iv_neg(i[0]), (x, s) => -x[0]),
new("iv_mul", 2, None, (i, s) => iv_mul(i[0], i[1]), (x, s) => x[0] * x[1]),
new("iv_scale", 1, r => new[] { Uniform(r, -4.0f, 4.0f) }, (i, s) => iv_scale(i[0], s[0]), (x, s) => x[0] * s[0]),
new("iv_abs", 1, None, (i, s) => iv_abs(i[0]), (x, s) => abs(x[0])),
new("iv_min", 2, None, (i, s) => iv_min(i[0], i[1]), (x, s) => min(x[0], x[1])),
new("iv_max", 2, None, (i, s) => iv_max(i[0], i[1]), (x, s) => max(x[0], x[1])),
new("iv_sqr", 1, None, (i, s) => iv_sqr(i[0]), (x, s) => x[0] * x[0]),
new("iv_sqrt", 1, None, (i, s) => iv_sqrt(i[0]), (x, s) => sqrt(x[0]), NonNegativeDomain: true),
new("iv_exp", 1, None, (i, s) => iv_exp(i[0]), (x, s) => exp(x[0])),
new("iv_clamp", 1, r => new[] { Uniform(r, -25.0f, 25.0f), Uniform(r, -25.0f, 25.0f) },
(i, s) => iv_clamp(i[0], s[0], s[1]), (x, s) => clamp(x[0], s[0], s[1])),
new("iv_saturate", 1, None, (i, s) => iv_saturate(i[0]), (x, s) => saturate(x[0])),
new("iv_smoothstep", 1, r => new[] { Uniform(r, -25.0f, 25.0f), Uniform(r, -25.0f, 25.0f) },
(i, s) => iv_smoothstep(s[0], s[1], i[0]), (x, s) => smoothstep(s[0], s[1], x[0])),
new("iv_lerp", 2, r => new[] { Uniform(r, -0.5f, 1.5f) }, (i, s) => iv_lerp(i[0], i[1], s[0]), (x, s) => lerp(x[0], x[1], s[0])),
};
// Centres in [-20, 20], widths log-uniform in [1e-4, 10]. Every eighth interval is a point
// interval: a bound that is only reached when an argument is a point (iv_lerp's, where x appears
// twice) is still reached, so a one-ulp shrink of it is still caught.
private static float2 RandomInterval(System.Random random, int index)
{
var centre = Uniform(random, -20.0f, 20.0f);
if (index % 8 == 7)
return iv_point(centre);
var width = MathF.Pow(10.0f, Uniform(random, -4.0f, 1.0f));
return new float2(centre - width * 0.5f, centre + width * 0.5f);
}
// The points inside one argument: both endpoints, zero when the interval holds it, and the rest
// uniform. Endpoints and zero are where a tight bound is reached.
private static float[] PointsInside(System.Random random, float2 interval, bool nonNegative)
{
var lo = nonNegative ? max(interval.x, 0.0f) : interval.x;
var hi = nonNegative ? max(interval.y, 0.0f) : interval.y;
var points = new float[PointsPerInterval];
points[0] = lo;
points[1] = hi;
var next = 2;
if (lo < 0.0f && hi > 0.0f)
points[next++] = 0.0f;
for (; next < points.Length; next++)
points[next] = min(max(lo + (hi - lo) * random.NextSingle(), lo), hi);
return points;
}
/// <summary>Runs one op over the seeded set; returns the violations, at most a few, described.</summary>
internal static List<string> Violations(Op op, Func<float2, float2>? mutate = null, int seed = 0x1A7E)
{
var random = new System.Random(seed);
var violations = new List<string>();
var count = 0;
for (var k = 0; k < IntervalCount; k++)
{
var intervals = new float2[op.Intervals];
for (var a = 0; a < op.Intervals; a++)
intervals[a] = RandomInterval(random, k + a);
var scalars = op.Scalars(random);
var bound = op.Bound(intervals, scalars);
if (mutate is not null)
bound = mutate(bound);
var points = intervals.Select(i => PointsInside(random, i, op.NonNegativeDomain)).ToArray();
for (var p = 0; p < PointsPerInterval; p++)
{
// Binary ops visit the four endpoint corners first (p = 0..3), then pair points by index.
var x = new float[op.Intervals];
for (var a = 0; a < op.Intervals; a++)
x[a] = op.Intervals == 2 && p < 4 ? points[a][a == 0 ? p & 1 : p >> 1] : points[a][p];
var value = op.Pointwise(x, scalars);
if (!(value >= bound.x && value <= bound.y))
{
if (count++ < 4)
violations.Add($"{op.Name}({string.Join(", ", intervals)}; {string.Join(", ", scalars)}) = {bound}, f({string.Join(", ", x)}) = {value:R}");
}
}
}
if (count > violations.Count)
violations.Add($"... {count} violations in all");
return violations;
}
[Fact]
public void EveryOpEnclosesItsPoints()
{
var failures = Ops.SelectMany(op => Violations(op)).ToList();
Assert.True(failures.Count == 0, string.Join(Environment.NewLine, failures));
}
internal static float NextUp(float v) => v == 0.0f ? asfloat(1u) : asfloat(v > 0.0f ? asuint(v) + 1u : asuint(v) - 1u);
private static float NextDown(float v) => -NextUp(-v);
/// <summary>
/// The hand-mutation list: every op's lo raised by one ulp, and its hi lowered by one ulp, must each
/// fail the harness. The two equivalent pairs are named in the class comment.
/// </summary>
[Fact]
public void EveryOneUlpShrinkIsCaught()
{
var survivors = new List<string>();
foreach (var op in Ops.Where(o => o.Name is not "iv_exp" and not "iv_smoothstep"))
{
if (Violations(op, b => new float2(NextUp(b.x), b.y)).Count == 0)
survivors.Add(op.Name + " lo + 1 ulp");
if (Violations(op, b => new float2(b.x, NextDown(b.y))).Count == 0)
survivors.Add(op.Name + " hi - 1 ulp");
}
Assert.True(survivors.Count == 0, "surviving one-ulp shrinks: " + string.Join(", ", survivors));
}
[Fact]
public void MulCoversEverySignCase()
{
// Both positive, both negative, mixed both ways, and a factor straddling zero on either side.
Assert.Equal(new float2(2.0f, 12.0f), iv_mul(new float2(1.0f, 3.0f), new float2(2.0f, 4.0f)));
Assert.Equal(new float2(2.0f, 12.0f), iv_mul(new float2(-3.0f, -1.0f), new float2(-4.0f, -2.0f)));
Assert.Equal(new float2(-12.0f, -2.0f), iv_mul(new float2(1.0f, 3.0f), new float2(-4.0f, -2.0f)));
Assert.Equal(new float2(-12.0f, -2.0f), iv_mul(new float2(-3.0f, -1.0f), new float2(2.0f, 4.0f)));
Assert.Equal(new float2(-8.0f, 12.0f), iv_mul(new float2(-2.0f, 3.0f), new float2(2.0f, 4.0f)));
Assert.Equal(new float2(-12.0f, 8.0f), iv_mul(new float2(-2.0f, 3.0f), new float2(-4.0f, -2.0f)));
Assert.Equal(new float2(-15.0f, 12.0f), iv_mul(new float2(-2.0f, 3.0f), new float2(-5.0f, 4.0f)));
}
[Fact]
public void SqrtClampsANegativeLoAtZero()
{
Assert.Equal(new float2(0.0f, 3.0f), iv_sqrt(new float2(-4.0f, 9.0f)));
Assert.Equal(new float2(0.0f, 0.0f), iv_sqrt(new float2(-4.0f, -1.0f)));
Assert.Equal(new float2(2.0f, 3.0f), iv_sqrt(new float2(4.0f, 9.0f)));
}
[Fact]
public void AbsAndSqrStraddlingZeroStartAtZero()
{
Assert.Equal(new float2(0.0f, 3.0f), iv_abs(new float2(-3.0f, 2.0f)));
Assert.Equal(new float2(0.0f, 5.0f), iv_abs(new float2(-2.0f, 5.0f)));
Assert.Equal(new float2(2.0f, 5.0f), iv_abs(new float2(-5.0f, -2.0f)));
Assert.Equal(new float2(0.0f, 9.0f), iv_sqr(new float2(-3.0f, 2.0f)));
Assert.Equal(new float2(4.0f, 25.0f), iv_sqr(new float2(-5.0f, -2.0f)));
}
[Fact]
public void MinAndMaxAreThePointwiseMinAndMaxOfTheEnds()
{
var a = new float2(-1.0f, 4.0f);
var b = new float2(2.0f, 3.0f);
Assert.Equal(new float2(-1.0f, 3.0f), iv_min(a, b));
Assert.Equal(new float2(2.0f, 4.0f), iv_max(a, b));
}
[Fact]
public void SmoothstepBoundsAreItsOrderedEndpointsWidenedByOneUlp()
{
var rising = iv_smoothstep(0.0f, 1.0f, new float2(0.25f, 0.5f));
Assert.Equal(NextDown(smoothstep(0.0f, 1.0f, 0.25f)), rising.x);
Assert.Equal(NextUp(smoothstep(0.0f, 1.0f, 0.5f)), rising.y);
// minimum > maximum: smoothstep falls, so the endpoints swap.
var falling = iv_smoothstep(1.0f, 0.0f, new float2(0.25f, 0.5f));
Assert.Equal(NextDown(smoothstep(1.0f, 0.0f, 0.5f)), falling.x);
Assert.Equal(NextUp(smoothstep(1.0f, 0.0f, 0.25f)), falling.y);
// The widening never leaves [0, 1].
Assert.Equal(new float2(0.0f, 1.0f), iv_smoothstep(0.0f, 1.0f, new float2(-2.0f, 3.0f)));
}
/// <summary>
/// smoothstep's float32 evaluation steps down by one ulp in places (MeasureMonotonicity). Three
/// consecutive floats t0 &lt; t1 &lt; t2 with f(t1) below both neighbours put f(t1) outside the ordered
/// endpoints of [t0, t2]; three with f(t1) above both put it outside on the other side. The widened
/// bound encloses both; without the widening it does not.
/// </summary>
[Fact]
public void SmoothstepWideningCoversItsFloatNonMonotonicity()
{
float F(float t) => smoothstep(0.0f, 1.0f, t);
float2? dip = null, peak = null;
for (var t = 1.0e-4f; t < 0.01f && (dip is null || peak is null); t = NextUp(t))
{
var t1 = NextUp(t);
var t2 = NextUp(t1);
if (dip is null && F(t1) < F(t) && F(t1) < F(t2))
dip = new float2(t, t2);
if (peak is null && F(t1) > F(t) && F(t1) > F(t2))
peak = new float2(t, t2);
}
Assert.NotNull(dip);
Assert.NotNull(peak);
foreach (var (interval, low) in new[] { (dip!.Value, true), (peak!.Value, false) })
{
var mid = F(NextUp(interval.x));
var unwidened = new float2(min(F(interval.x), F(interval.y)), max(F(interval.x), F(interval.y)));
Assert.False(mid >= unwidened.x && mid <= unwidened.y, $"the {(low ? "dip" : "peak")} at {interval} is inside the endpoints");
var bound = iv_smoothstep(0.0f, 1.0f, interval);
Assert.True(mid >= bound.x && mid <= bound.y, $"iv_smoothstep{interval} = {bound} misses {mid:R}");
}
}
[Fact]
public void LerpIsTheNaturalExtensionAndTightForPointArguments()
{
Assert.Equal(new float2(lerp(1.0f, 3.0f, 0.25f), lerp(1.0f, 7.0f, 0.25f)), iv_lerp(iv_point(1.0f), new float2(3.0f, 7.0f), 0.25f));
Assert.Equal(new float2(-1.0f + (2.0f - 1.0f) * 0.5f, 1.0f + (4.0f + 1.0f) * 0.5f), iv_lerp(new float2(-1.0f, 1.0f), new float2(2.0f, 4.0f), 0.5f));
}
[Fact]
public void ExpIsWidenedByOneUlpOnEachSide()
{
var bound = iv_exp(new float2(-1.0f, 2.0f));
Assert.Equal(NextDown(exp(-1.0f)), bound.x);
Assert.Equal(NextUp(exp(2.0f)), bound.y);
Assert.Equal(0.0f, iv_exp(new float2(-200.0f, 0.0f)).x);
}
/// <summary>
/// The monotonicity the widening rules rest on, measured over every float: exp over [-87, 88.7], the
/// range where it neither underflows to a subnormal nor overflows, and smoothstep(0, 1, t) over [0, 1],
/// which covers every smoothstep because t = saturate((x - minimum) / (maximum - minimum)) is itself
/// monotone in x. Prints the number of steps down and the largest, in ulps. Slow: about a minute.
/// </summary>
[Fact(Explicit = true)]
[Trait("Category", "Slow")]
public void MeasureMonotonicity()
{
long expDown = 0, expWorst = 0;
var previous = exp(-87.0f);
for (var x = NextUp(-87.0f); x <= 88.7f; x = NextUp(x))
{
var e = exp(x);
if (e < previous)
{
expDown++;
expWorst = Math.Max(expWorst, (long)asuint(previous) - asuint(e));
}
previous = e;
}
long smoothDown = 0, smoothWorst = 0;
var runningMax = 0.0f;
for (var bits = 0u; bits <= asuint(1.0f); bits++)
{
var s = smoothstep(0.0f, 1.0f, asfloat(bits));
if (s < runningMax)
{
smoothDown++;
smoothWorst = Math.Max(smoothWorst, (long)asuint(runningMax) - asuint(s));
}
else
{
runningMax = s;
}
}
output.WriteLine($"IV-REPORT monotonicity: exp steps down {expDown} times, worst {expWorst} ulp; smoothstep steps down {smoothDown} times, worst {smoothWorst} ulp");
Assert.True(expWorst <= 1, $"exp steps down by {expWorst} ulp; one ulp of widening is not enough");
Assert.True(smoothWorst <= 1, $"smoothstep steps down by {smoothWorst} ulp; one ulp of widening is not enough");
}
}

427
packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs Normal file
View file

@ -0,0 +1,427 @@
using CultMath;
using Xunit;
using static CultMath.math;
namespace CultMath.Tests;
/// <summary>
/// The noise bounds (design.md, "Intervals"): iv_snoise_ball and iv_fbm_ball enclose their functions over
/// a ball, a domain warp bounded by D is enclosed by enlarging the radius by D, SNOISE_LIPSCHITZ has
/// measured provenance, and an interval-culled march spends fewer snoise evaluations than a dense one for
/// the same transmittance. Lines starting IV-REPORT carry the numbers the cut report quotes.
/// </summary>
public sealed class NoiseBoundTests
{
private readonly ITestOutputHelper output;
public NoiseBoundTests(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 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);
}
// The first eight points of a ball are on its surface, where the bound is closest to binding; the
// rest are uniform in its volume.
private static float3 PointInBall(System.Random random, float3 centre, float radius, int index) =>
centre + UnitVector(random) * (index < 8 ? radius : radius * MathF.Cbrt(random.NextSingle()));
private static float3 RandomCentre(System.Random random) =>
new(Uniform(random, -20.0f, 20.0f), Uniform(random, -20.0f, 20.0f), Uniform(random, -20.0f, 20.0f));
private static bool Inside(float value, float2 bound) => value >= bound.x && value <= bound.y;
[Fact]
public void SnoiseBallEnclosesPoints()
{
var random = new System.Random(0xBA11);
for (var b = 0; b < 2000; b++)
{
var centre = RandomCentre(random);
var radius = LogUniform(random, 1.0e-3f, 2.0f);
var bound = iv_snoise_ball(centre, radius);
for (var i = 0; i < 64; i++)
{
var x = PointInBall(random, centre, radius, i);
Assert.True(Inside(snoise(x), bound), $"iv_snoise_ball({centre}, {radius:R}) = {bound} misses snoise({x}) = {snoise(x):R}");
}
}
}
[Fact]
public void FbmBallEnclosesPoints()
{
var random = new System.Random(0xFB11);
for (var b = 0; b < 2000; b++)
{
var centre = RandomCentre(random);
var radius = LogUniform(random, 1.0e-3f, 2.0f);
var octaves = random.Next(0, 7);
var lacunarity = Uniform(random, 1.5f, 2.5f);
var gain = Uniform(random, 0.3f, 0.7f);
var bound = iv_fbm_ball(centre, radius, octaves, lacunarity, gain);
for (var i = 0; i < 64; i++)
{
var x = PointInBall(random, centre, radius, i);
var value = fbm_grad(x, octaves, lacunarity, gain).w;
Assert.True(Inside(value, bound), $"iv_fbm_ball({centre}, {radius:R}, {octaves}, {lacunarity:R}, {gain:R}) = {bound} misses {value:R} at {x}");
}
}
}
/// <summary>
/// The site's march rests on this: a point x in the ball (c, r), displaced by a flow of length at most
/// D, lies in the ball (c, r + D), so snoise of the warped point lies in iv_snoise_ball(c, r + D). The
/// geometric check allows one part in 1e5 for the rounding of x + flow.
/// </summary>
[Fact]
public void WarpedPointsStayEnclosed()
{
var random = new System.Random(0x3A4B);
for (var b = 0; b < 2000; b++)
{
var centre = RandomCentre(random);
var radius = LogUniform(random, 1.0e-3f, 2.0f);
var warp = Uniform(random, 0.0f, 1.0f);
var bound = iv_snoise_ball(centre, radius + warp);
for (var i = 0; i < 64; i++)
{
var x = PointInBall(random, centre, radius, i);
var flow = UnitVector(random) * (warp * MathF.Sqrt(random.NextSingle()));
var warped = x + flow;
Assert.True(length(warped - centre) <= (radius + warp) * (1.0f + 1.0e-5f), $"warped point {warped} leaves the ball ({centre}, {radius + warp:R})");
Assert.True(Inside(snoise(warped), bound), $"iv_snoise_ball({centre}, {radius + warp:R}) = {bound} misses snoise({warped}) = {snoise(warped):R}");
}
}
}
// ---- The Lipschitz constant ----
private static float GradientNorm(float3 p)
{
var g = snoise_grad(p);
return MathF.Sqrt(g.x * g.x + g.y * g.y + g.z * g.z);
}
// Hill-climbs |snoise_grad| from p: the ascent direction is the central difference of |grad|, the
// step grows on success and halves on failure. Returns the largest |grad| it reached.
private static float AscendGradientNorm(float3 p)
{
const float h = 1.0e-3f;
var step = 1.0e-2f;
var best = GradientNorm(p);
for (var k = 0; k < 200 && step > 1.0e-6f; k++)
{
var d = new float3(
GradientNorm(p + new float3(h, 0.0f, 0.0f)) - GradientNorm(p - new float3(h, 0.0f, 0.0f)),
GradientNorm(p + new float3(0.0f, h, 0.0f)) - GradientNorm(p - new float3(0.0f, h, 0.0f)),
GradientNorm(p + new float3(0.0f, 0.0f, h)) - GradientNorm(p - new float3(0.0f, 0.0f, h)));
if (dot(d, d) == 0.0f)
break;
var q = p + normalize(d) * step;
var value = GradientNorm(q);
if (value > best)
{
p = q;
best = value;
step *= 1.5f;
}
else
{
step *= 0.5f;
}
}
return best;
}
// Same climb on |snoise|, along the value's own analytic gradient.
private static float AscendAbsValue(float3 p)
{
var step = 1.0e-2f;
var best = MathF.Abs(snoise(p));
for (var k = 0; k < 200 && step > 1.0e-6f; k++)
{
var g = snoise_grad(p);
var d = new float3(g.x, g.y, g.z) * (g.w < 0.0f ? -1.0f : 1.0f);
if (dot(d, d) == 0.0f)
break;
var q = p + normalize(d) * step;
var value = MathF.Abs(snoise(q));
if (value > best)
{
p = q;
best = value;
step *= 1.5f;
}
else
{
step *= 0.5f;
}
}
return best;
}
// The largest |snoise_grad| over `samples` seeded points in [-256, 256]^3, the largest after refining
// the `refine` largest of them by ascent, and the start whose ascent reached it (the witness).
private static (float Sampled, float Refined, float3 Witness) MeasureGradient(int seed, int samples, int refine)
{
var random = new System.Random(seed);
var points = new float3[samples];
var norms = new float[samples];
for (var i = 0; i < samples; i++)
{
points[i] = new float3(Uniform(random, -256.0f, 256.0f), Uniform(random, -256.0f, 256.0f), Uniform(random, -256.0f, 256.0f));
norms[i] = GradientNorm(points[i]);
}
var order = Enumerable.Range(0, samples).OrderByDescending(i => norms[i]).Take(Math.Max(refine, 1)).ToArray();
var refined = norms[order[0]];
var witness = points[order[0]];
foreach (var i in order.Take(refine))
{
var value = AscendGradientNorm(points[i]);
if (value > refined)
{
refined = value;
witness = points[i];
}
}
return (norms[order[0]], refined, witness);
}
/// <summary>
/// The provenance of SNOISE_LIPSCHITZ: the largest |snoise_grad| over 1e6 seeded points, the 1e4
/// largest refined by gradient ascent, times 1.10. Prints the start whose ascent reached the maximum,
/// which LipschitzConstantPinsSampledGradients climbs from again. Also prints the largest |snoise|
/// reached the same way, which is what lets iv_snoise_ball intersect with [-1, 1]. Slow; run
/// explicitly after any change to the snoise kernel and re-pin the constant and the witness.
/// </summary>
[Fact(Explicit = true)]
[Trait("Category", "Slow")]
public void MeasureLipschitz()
{
var (sampled, refined, witness) = MeasureGradient(0x11B5, 1_000_000, 10_000);
var random = new System.Random(0xAB5);
var valueStarts = Enumerable.Range(0, 10_000)
.Select(_ => new float3(Uniform(random, -256.0f, 256.0f), Uniform(random, -256.0f, 256.0f), Uniform(random, -256.0f, 256.0f)));
var maxValue = valueStarts.Max(AscendAbsValue);
output.WriteLine($"IV-REPORT lipschitz: max |snoise_grad| sampled {sampled:R}, refined {refined:R}, x 1.10 = {refined * 1.10f:R}, witness ({witness.x:R}, {witness.y:R}, {witness.z:R}); max |snoise| refined {maxValue:R}");
Assert.True(refined <= SNOISE_LIPSCHITZ, $"measured {refined:R} exceeds SNOISE_LIPSCHITZ {SNOISE_LIPSCHITZ:R}");
Assert.True(maxValue < 1.0f, $"|snoise| reaches {maxValue:R}; iv_snoise_ball may not intersect with [-1, 1]");
}
// The start MeasureLipschitz printed as its witness.
private static readonly float3 LipschitzWitness = new(-241.98447f, -249.98927f, -60.058212f);
/// <summary>
/// Pins SNOISE_LIPSCHITZ: no |snoise_grad| over 1e5 seeded points exceeds it, and the ascent from
/// MeasureLipschitz's witness reaches a maximum that, times 1.10, is the constant to within half a
/// percent, so a constant lowered (or raised) by 1% fails here.
/// </summary>
[Fact]
public void LipschitzConstantPinsSampledGradients()
{
var (sampled, _, _) = MeasureGradient(0x11B6, 100_000, 0);
Assert.True(sampled <= SNOISE_LIPSCHITZ, $"|snoise_grad| reaches {sampled:R}, above SNOISE_LIPSCHITZ {SNOISE_LIPSCHITZ:R}");
var refined = AscendGradientNorm(LipschitzWitness);
Assert.True(refined <= SNOISE_LIPSCHITZ, $"|snoise_grad| reaches {refined:R}, above SNOISE_LIPSCHITZ {SNOISE_LIPSCHITZ:R}");
var ratio = SNOISE_LIPSCHITZ / (1.10f * refined);
Assert.True(MathF.Abs(ratio - 1.0f) <= 0.005f, $"SNOISE_LIPSCHITZ is {ratio:R} times 1.10 x the re-measured {refined:R}; re-run MeasureLipschitz and re-pin");
}
/// <summary>
/// Mean of interval width over the true range (max - min of the 64 sampled values) per ball: how
/// loose the Lipschitz bound is, for a later affine cut to beat. Asserts only finiteness.
/// </summary>
[Fact]
public void TightnessReport()
{
double Ratio(Func<float3, float, float2> bound, Func<float3, float> f, float radius, int seed)
{
var random = new System.Random(seed);
var sum = 0.0;
for (var b = 0; b < 2000; b++)
{
var centre = RandomCentre(random);
var interval = bound(centre, radius);
float lo = float.MaxValue, hi = float.MinValue;
for (var i = 0; i < 64; i++)
{
var v = f(PointInBall(random, centre, radius, i));
lo = MathF.Min(lo, v);
hi = MathF.Max(hi, v);
}
sum += (interval.y - interval.x) / (double)(hi - lo);
}
return sum / 2000;
}
var snoiseRatios = new[] { 0.05f, 0.25f, 1.0f }.Select(r => (r, Ratio(iv_snoise_ball, snoise, r, 0x7167))).ToArray();
var fbmRatios = new[] { 0.05f, 0.25f }.Select(r => (r, Ratio((c, rr) => iv_fbm_ball(c, rr, 4, 2.0f, 0.5f), x => fbm_grad(x, 4, 2.0f, 0.5f).w, r, 0x7168))).ToArray();
output.WriteLine("IV-REPORT tightness (mean width / sampled range): " +
string.Join(", ", snoiseRatios.Select(t => $"snoise r={t.r:R} {t.Item2:F2}")) + "; " +
string.Join(", ", fbmRatios.Select(t => $"fbm4 r={t.r:R} {t.Item2:F2}")));
Assert.All(snoiseRatios.Concat(fbmRatios), t => Assert.True(double.IsFinite(t.Item2)));
}
// ---- The saving ----
/// <summary>
/// The map's field shape (docs/cultmath-interval-ground-cut.md, "Density field"), in world units along
/// a ray through a slab z in [0, Depth]: coverage is two snoise of the domain-warped point, cross-faded
/// over two half-period phases; occupancy = max(0, coverage - Cutoff) * Extinction; one more snoise is
/// the detail term. The flow is a unit-bounded 2D field times FlowAmplitude, so a warp moves a point by
/// at most Warp = FlowAmplitude * Period / 2. Flow evaluations are not counted: the site reads flow from
/// a texture. Counted: three snoise per dense sample, one per interval probe.
/// </summary>
internal sealed class SavingField
{
public const float F0 = 0.35f;
public const float Cutoff = 0.25f;
public const float Depth = 12.0f;
public const float Period = 4.0f;
public const float Extinction = 2.0f;
public const float FlowFrequency = 0.15f;
public float Step0 { get; init; } = 0.25f;
public float Warp { get; init; } = 0.05f;
public float StepDense => Step0 * 0.5f;
public float StepMax => Step0 * 8.0f;
public float FlowAmplitude => Warp * 2.0f / Period;
public long Evaluations;
private float2 Flow(float3 p)
{
var q = new float2(p.x, p.y) * FlowFrequency;
var v = new float2(snoise(q + new float2(17.0f, 3.0f)), snoise(q + new float2(-5.0f, 41.0f)));
return v / max(1.0f, length(v)) * FlowAmplitude;
}
public float Density(float3 p, float time)
{
var flow = Flow(p);
var phase0 = frac(time / Period);
var phase1 = frac(time / Period + 0.5f);
var w0 = 1.0f - abs(2.0f * phase0 - 1.0f);
var coverage = w0 * snoise((p + new float3(flow * ((phase0 - 0.5f) * Period), 0.0f)) * F0)
+ (1.0f - w0) * snoise((p + new float3(flow * ((phase1 - 0.5f) * Period), 0.0f)) * F0);
var detail = 1.0f - abs(snoise(p * (4.0f * F0) + 17.0f));
Evaluations += 3;
return max(0.0f, coverage - Cutoff) * Extinction * lerp(0.6f, 1.0f, detail);
}
// Transmittance through the dense grid: cells of StepDense in z, one sample at each midpoint.
public float MarchDense(float3 origin, float3 direction, float time)
{
var transmittance = 1.0f;
for (var z = 0.0f; z < Depth && transmittance >= 0.02f; z += StepDense)
transmittance *= Integrate(origin, direction, time, z, StepDense);
return transmittance;
}
private float Integrate(float3 origin, float3 direction, float time, float z, float step) =>
exp(-Density(origin + direction * ((z + step * 0.5f) / direction.z), time) * (step / direction.z));
// The map's adaptive march: probe a segment with one iv_snoise_ball; provably empty skips and grows,
// maybe-dense halves down to StepDense, where the cell is integrated exactly as the dense march does.
// Steps are StepDense times powers of two, so every integrated cell is a cell of the dense grid.
public float MarchIntervals(float3 origin, float3 direction, float time)
{
var transmittance = 1.0f;
var z = 0.0f;
var step = Step0;
while (z < Depth && transmittance >= 0.02f)
{
while (z + step > Depth)
step *= 0.5f;
var centre = origin + direction * ((z + step * 0.5f) / direction.z);
var radius = step * 0.5f / direction.z + Warp;
Evaluations++;
if (iv_snoise_ball(centre * F0, radius * F0).y < Cutoff)
{
z += step;
step = min(step * 2.0f, StepMax);
}
else if (step > StepDense)
{
step *= 0.5f;
}
else
{
transmittance *= Integrate(origin, direction, time, z, step);
z += step;
}
}
return transmittance;
}
public (long Dense, long Intervals, float MaxDifference) Run(int rays, int seed)
{
var random = new System.Random(seed);
long dense = 0, intervals = 0;
var maxDifference = 0.0f;
for (var r = 0; r < rays; r++)
{
var origin = new float3(Uniform(random, -50.0f, 50.0f), Uniform(random, -50.0f, 50.0f), 0.0f);
var direction = normalize(new float3(Uniform(random, -0.5f, 0.5f), Uniform(random, -0.5f, 0.5f), 1.0f));
var time = Uniform(random, 0.0f, 100.0f);
Evaluations = 0;
var a = MarchDense(origin, direction, time);
dense += Evaluations;
Evaluations = 0;
var b = MarchIntervals(origin, direction, time);
intervals += Evaluations;
maxDifference = MathF.Max(maxDifference, MathF.Abs(a - b));
}
return (dense, intervals, maxDifference);
}
}
/// <summary>
/// The interval march must reach the dense march's transmittance within 1e-3 (it integrates the same
/// cells, so any difference is a skipped cell that was not empty) with at most half its snoise
/// evaluations. Also prints the ratio over a small grid of first steps and warp bounds (300 rays each).
/// The 2x saving does not hold at SNOISE_LIPSCHITZ as measured; that is an open fork for the campaign
/// (cultmath-tapes, interval-ops cut report), so a ratio below 2 skips with the fork named rather than
/// failing or being tuned away.
/// </summary>
[Fact]
public void IntervalSkipHalvesEvaluations()
{
var field = new SavingField();
var (dense, intervals, maxDifference) = field.Run(1000, 0x5A7E);
output.WriteLine($"IV-REPORT saving: dense {dense} snoise, interval march {intervals} snoise, ratio {(double)dense / intervals:F2}; max transmittance difference {maxDifference:R} (step0 {field.Step0:R}, warp {field.Warp:R}, L {SNOISE_LIPSCHITZ:R})");
foreach (var step0 in new[] { 0.0625f, 0.25f, 1.0f })
foreach (var warp in new[] { 0.0f, 0.05f })
{
var (d, i, m) = new SavingField { Step0 = step0, Warp = warp }.Run(300, 0x5A7F);
Assert.True(m <= 1.0e-3f, $"transmittance differs by {m:R} at step0 {step0:R}, warp {warp:R}");
output.WriteLine($"IV-REPORT saving grid: step0 {step0:R} warp {warp:R}: dense {d}, interval march {i}, ratio {(double)d / i:F2}");
}
Assert.True(maxDifference <= 1.0e-3f, $"transmittance differs by {maxDifference:R}");
if (dense < 2 * intervals)
Assert.Skip($"FORK (cultmath-tapes interval-ops): the interval march saves {(double)dense / intervals:F2}x, not 2x, at SNOISE_LIPSCHITZ {SNOISE_LIPSCHITZ:R}");
}
}