From bbf52697515274736b093220e0f90bf91c615433 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Fri, 2 Oct 2026 19:52:33 +0200 Subject: [PATCH 01/14] CultMath intervals: cultmath_iv_* over float2(lo, hi) in C# and HLSL, with snoise and fbm ball bounds math.Interval.cs and shaders/CultMath.Interval.hlsl (included from CultMath.hlsl after Phacelle; both tracked copies) add iv_point, iv_add, iv_sub, iv_neg, iv_mul, iv_scale, iv_abs, iv_min, iv_max, iv_sqr, iv_sqrt, iv_exp, iv_clamp, iv_saturate, iv_smoothstep, iv_lerp, iv_snoise_ball and iv_fbm_ball, and SNOISE_LIPSCHITZ = 10.099261 (1.10 x the refined max |snoise_grad| that NoiseBoundTests.MeasureLipschitz prints). iv_exp and iv_smoothstep widen by one ulp through asfloat(uint), added beside asuint as the HLSL intrinsic. The mirror test now compares 52 functions; the Unity package and the NuGet content files carry the new include. Co-Authored-By: Claude Opus 5.5 --- .../cultmath/scripts/build-unity-package.ps1 | 4 +- .../cultmath/shaders/CultMath.Interval.hlsl | 107 ++++++++++++ .../shaders/CultMath.Interval.hlsl.meta | 7 + packages/cultmath/shaders/CultMath.hlsl | 1 + .../cultmath/src/CultMath/CultMath.csproj | 1 + .../cultmath/src/CultMath/math.Interval.cs | 165 ++++++++++++++++++ .../src/CultMath/math.Interval.cs.meta | 2 + packages/cultmath/src/CultMath/math.cs | 3 + .../tests/CultMath.Tests/HlslMirrorTests.cs | 1 + .../HlslSourceCompatibilityTests.cs | 4 +- .../Shaders/CultMath.Interval.hlsl | 107 ++++++++++++ .../Shaders/CultMath.Interval.hlsl.meta | 7 + .../Shaders/CultMath.hlsl | 1 + 13 files changed, 408 insertions(+), 2 deletions(-) create mode 100644 packages/cultmath/shaders/CultMath.Interval.hlsl create mode 100644 packages/cultmath/shaders/CultMath.Interval.hlsl.meta create mode 100644 packages/cultmath/src/CultMath/math.Interval.cs create mode 100644 packages/cultmath/src/CultMath/math.Interval.cs.meta create mode 100644 packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl create mode 100644 packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl.meta diff --git a/packages/cultmath/scripts/build-unity-package.ps1 b/packages/cultmath/scripts/build-unity-package.ps1 index 228d735c..2209ec66 100644 --- a/packages/cultmath/scripts/build-unity-package.ps1 +++ b/packages/cultmath/scripts/build-unity-package.ps1 @@ -80,7 +80,9 @@ $requiredTemplateFiles = @( "Shaders\CultMath.hlsl", "Shaders\CultMath.hlsl.meta", "Shaders\CultMath.Phacelle.hlsl", - "Shaders\CultMath.Phacelle.hlsl.meta" + "Shaders\CultMath.Phacelle.hlsl.meta", + "Shaders\CultMath.Interval.hlsl", + "Shaders\CultMath.Interval.hlsl.meta" ) foreach ($relativePath in $requiredTemplateFiles) { if (-not (Test-Path -LiteralPath (Join-Path $templateRoot $relativePath) -PathType Leaf)) { diff --git a/packages/cultmath/shaders/CultMath.Interval.hlsl b/packages/cultmath/shaders/CultMath.Interval.hlsl new file mode 100644 index 00000000..ed6f83db --- /dev/null +++ b/packages/cultmath/shaders/CultMath.Interval.hlsl @@ -0,0 +1,107 @@ +#ifndef CULTMATH_INTERVAL_HLSL +#define CULTMATH_INTERVAL_HLSL + +// Interval arithmetic over float2(lo, hi), lo <= hi, both finite. Every cultmath_iv_* function returns +// an interval containing f(x) for every x in its input intervals, mirroring math.Interval.cs bit for +// bit; see that file and docs/design.md ("Intervals") for each op's enclosure and its ulp widening. +// Included by CultMath.hlsl after its own definitions of cultmath_smoothstep and cultmath_snoise. + +// Lipschitz constant of cultmath_snoise(float3); provenance in math.Interval.cs (SNOISE_LIPSCHITZ). +static const float CULTMATH_SNOISE_LIPSCHITZ = 10.099261; + +float2 cultmath_iv_point(float value) { return float2(value, value); } +float2 cultmath_iv_add(float2 a, float2 b) { return float2(a.x + b.x, a.y + b.y); } +float2 cultmath_iv_sub(float2 a, float2 b) { return float2(a.x - b.y, a.y - b.x); } +float2 cultmath_iv_neg(float2 a) { return float2(-a.y, -a.x); } + +float2 cultmath_iv_mul(float2 a, float2 b) +{ + float p0 = a.x * b.x; + float p1 = a.x * b.y; + float p2 = a.y * b.x; + float p3 = a.y * b.y; + return float2(min(min(p0, p1), min(p2, p3)), max(max(p0, p1), max(p2, p3))); +} + +float2 cultmath_iv_scale(float2 a, float s) { return s >= 0.0 ? float2(a.x * s, a.y * s) : float2(a.y * s, a.x * s); } + +float2 cultmath_iv_abs(float2 a) +{ + if (a.x >= 0.0) + return a; + if (a.y <= 0.0) + return float2(-a.y, -a.x); + return float2(0.0, max(-a.x, a.y)); +} + +float2 cultmath_iv_min(float2 a, float2 b) { return float2(min(a.x, b.x), min(a.y, b.y)); } +float2 cultmath_iv_max(float2 a, float2 b) { return float2(max(a.x, b.x), max(a.y, b.y)); } + +float2 cultmath_iv_sqr(float2 a) +{ + float l = a.x * a.x; + float h = a.y * a.y; + if (a.x >= 0.0) + return float2(l, h); + if (a.y <= 0.0) + return float2(h, l); + return float2(0.0, max(l, h)); +} + +float2 cultmath_iv_sqrt(float2 a) { return float2(sqrt(max(a.x, 0.0)), sqrt(max(a.y, 0.0))); } + +// Widened by one ulp on each side: exp is not required to be correctly rounded. +float2 cultmath_iv_exp(float2 a) +{ + float lo = exp(a.x); + float hi = exp(a.y); + return float2(lo > 0.0 ? asfloat(asuint(lo) - 1u) : 0.0, asfloat(asuint(hi) + 1u)); +} + +float2 cultmath_iv_clamp(float2 a, float minimum, float maximum) +{ + return float2(cultmath_clamp(a.x, minimum, maximum), cultmath_clamp(a.y, minimum, maximum)); +} + +float2 cultmath_iv_saturate(float2 a) { return float2(cultmath_saturate(a.x), cultmath_saturate(a.y)); } + +// The endpoints, ordered, widened by one ulp inside [0, 1]: the float32 evaluation of smoothstep is +// monotone only to within one ulp. +float2 cultmath_iv_smoothstep(float minimum, float maximum, float2 a) +{ + float s0 = cultmath_smoothstep(minimum, maximum, a.x); + float s1 = cultmath_smoothstep(minimum, maximum, a.y); + float lo = min(s0, s1); + float hi = max(s0, s1); + return float2(lo > 0.0 ? asfloat(asuint(lo) - 1u) : 0.0, hi < 1.0 ? asfloat(asuint(hi) + 1u) : 1.0); +} + +// The natural interval extension of x + (y - x) * amount, which is how cultmath_lerp evaluates. +float2 cultmath_iv_lerp(float2 a, float2 b, float amount) { return cultmath_iv_add(a, cultmath_iv_scale(cultmath_iv_sub(b, a), amount)); } + +// [n - L r, n + L r] intersected with [-1, 1], n = cultmath_snoise(centre). One snoise, value only. +float2 cultmath_iv_snoise_ball(float3 centre, float radius) +{ + float n = cultmath_snoise(centre); + float e = CULTMATH_SNOISE_LIPSCHITZ * radius; + return float2(max(n - e, -1.0), min(n + e, 1.0)); +} + +// Octave sum of cultmath_iv_snoise_ball, compounding frequency and amplitude as cultmath_fbm_grad does. +float2 cultmath_iv_fbm_ball(float3 centre, float radius, int octaves, float lacunarity, float gain) +{ + octaves = clamp(octaves, 0, 16); + float amplitude = 1.0; + float frequency = 1.0; + float2 sum = float2(0.0, 0.0); + for (int i = 0; i < octaves; i++) + { + sum = cultmath_iv_add(sum, cultmath_iv_scale(cultmath_iv_snoise_ball(centre * frequency, radius * abs(frequency)), amplitude)); + frequency *= lacunarity; + amplitude *= gain; + } + + return sum; +} + +#endif diff --git a/packages/cultmath/shaders/CultMath.Interval.hlsl.meta b/packages/cultmath/shaders/CultMath.Interval.hlsl.meta new file mode 100644 index 00000000..a22ec627 --- /dev/null +++ b/packages/cultmath/shaders/CultMath.Interval.hlsl.meta @@ -0,0 +1,7 @@ +fileFormatVersion: 2 +guid: 4d8a899561654ac990cab0bf48658cc1 +ShaderIncludeImporter: + externalObjects: {} + userData: + assetBundleName: + assetBundleVariant: diff --git a/packages/cultmath/shaders/CultMath.hlsl b/packages/cultmath/shaders/CultMath.hlsl index af2d3741..da3362a6 100644 --- a/packages/cultmath/shaders/CultMath.hlsl +++ b/packages/cultmath/shaders/CultMath.hlsl @@ -504,6 +504,7 @@ CultCellular cultmath_cellular(float3 p) } #include "CultMath.Phacelle.hlsl" +#include "CultMath.Interval.hlsl" float cultmath_value_noise(float2 position) { diff --git a/packages/cultmath/src/CultMath/CultMath.csproj b/packages/cultmath/src/CultMath/CultMath.csproj index 0b1ee09f..2a82f8db 100644 --- a/packages/cultmath/src/CultMath/CultMath.csproj +++ b/packages/cultmath/src/CultMath/CultMath.csproj @@ -24,6 +24,7 @@ + diff --git a/packages/cultmath/src/CultMath/math.Interval.cs b/packages/cultmath/src/CultMath/math.Interval.cs new file mode 100644 index 00000000..1430e385 --- /dev/null +++ b/packages/cultmath/src/CultMath/math.Interval.cs @@ -0,0 +1,165 @@ +namespace CultMath; + +public static partial class math +{ + // Interval arithmetic over float2(lo, hi), lo <= hi, both finite (design.md, "Intervals"). Every + // iv_* function returns an interval that contains f(x) for every x in its input intervals, where + // f(x) is the same function evaluated pointwise in float32, mirrored by cultmath_iv_* in + // shaders/CultMath.Interval.hlsl. Each bound is computed with the pointwise function's own float32 + // operations; because IEEE round-to-nearest is monotone, an op that is monotone on its region + // encloses its float32 evaluation with no widening. Where float32 evaluation is not monotone, or + // the libm is not required to be, the bound is widened by exactly one ulp, and the doc comment + // says so. There is no empty interval: a caller that needs "provably empty" tests hi < cutoff. + + // Lipschitz constant of snoise(float3) in its input: an upper bound on |snoise_grad(p).xyz|. Provenance: + // NoiseBoundTests.MeasureLipschitz (slow, explicit) takes the largest |snoise_grad| over 1e6 seeded + // points, refines the 1e4 largest by gradient ascent, and multiplies the refined maximum by 1.10; + // NoiseBoundTests.LipschitzConstantPinsSampledGradients pins this value against a fast re-measure. + // Empirical with a margin, not a proof; the enclosure tests are the defence. A change to the + // snoise kernel must re-run MeasureLipschitz and re-pin this constant. + public const float SNOISE_LIPSCHITZ = 10.099261f; + + /// The point interval [x, x]. + public static float2 iv_point(float value) => new(value, value); + + /// Encloses x + y for x in a, y in b. Monotone: no widening. + public static float2 iv_add(float2 a, float2 b) => new(a.x + b.x, a.y + b.y); + + /// Encloses x - y for x in a, y in b. Monotone: no widening. + public static float2 iv_sub(float2 a, float2 b) => new(a.x - b.y, a.y - b.x); + + /// Encloses -x for x in a. Exact. + public static float2 iv_neg(float2 a) => new(-a.y, -a.x); + + /// + /// Encloses x * y for x in a, y in b: the least and greatest of the four corner products. The + /// product is monotone in each factor for a fixed sign of the other, so no widening. + /// + public static float2 iv_mul(float2 a, float2 b) + { + var p0 = a.x * b.x; + var p1 = a.x * b.y; + var p2 = a.y * b.x; + var p3 = a.y * b.y; + return new(min(min(p0, p1), min(p2, p3)), max(max(p0, p1), max(p2, p3))); + } + + /// Encloses x * s for x in a. Monotone: no widening. + public static float2 iv_scale(float2 a, float s) => s >= 0.0f ? new(a.x * s, a.y * s) : new(a.y * s, a.x * s); + + /// Encloses abs(x) for x in a; an interval straddling zero gives [0, max(-lo, hi)]. Exact. + public static float2 iv_abs(float2 a) + { + if (a.x >= 0.0f) + return a; + if (a.y <= 0.0f) + return new(-a.y, -a.x); + return new(0.0f, max(-a.x, a.y)); + } + + /// Encloses min(x, y) for x in a, y in b. Monotone: no widening. + public static float2 iv_min(float2 a, float2 b) => new(min(a.x, b.x), min(a.y, b.y)); + + /// Encloses max(x, y) for x in a, y in b. Monotone: no widening. + public static float2 iv_max(float2 a, float2 b) => new(max(a.x, b.x), max(a.y, b.y)); + + /// + /// Encloses x * x for x in a; an interval straddling zero gives [0, max(lo * lo, hi * hi)]. Tighter + /// than iv_mul(a, a), which cannot know both factors are the same x. Monotone on each side of zero: + /// no widening. + /// + public static float2 iv_sqr(float2 a) + { + var l = a.x * a.x; + var h = a.y * a.y; + if (a.x >= 0.0f) + return new(l, h); + if (a.y <= 0.0f) + return new(h, l); + return new(0.0f, max(l, h)); + } + + /// + /// Encloses sqrt(x) for x in a with x >= 0 (sqrt's domain): lo is clamped at 0, so [-4, 9] gives [0, 3] + /// and an interval wholly below zero gives [0, 0]. sqrt is correctly rounded, hence monotone: no + /// widening. + /// + public static float2 iv_sqrt(float2 a) => new(sqrt(max(a.x, 0.0f)), sqrt(max(a.y, 0.0f))); + + /// + /// Encloses exp(x) for x in a. Widened by one ulp on each side, because exp is not required to be + /// correctly rounded, so a libm need not be monotone. Where exp(hi) overflows to infinity, hi is NaN, + /// which a cull test hi < cutoff reads as "not provably below", the safe side. + /// + public static float2 iv_exp(float2 a) + { + var lo = exp(a.x); + var hi = exp(a.y); + return new(lo > 0.0f ? asfloat(asuint(lo) - 1u) : 0.0f, asfloat(asuint(hi) + 1u)); + } + + /// Encloses clamp(x, minimum, maximum) for x in a. Monotone: no widening. + public static float2 iv_clamp(float2 a, float minimum, float maximum) => + new(clamp(a.x, minimum, maximum), clamp(a.y, minimum, maximum)); + + /// Encloses saturate(x) for x in a. Monotone: no widening. + public static float2 iv_saturate(float2 a) => new(saturate(a.x), saturate(a.y)); + + /// + /// Encloses smoothstep(minimum, maximum, x) for x in a. smoothstep is monotone in x (increasing when + /// minimum < maximum, decreasing when minimum > maximum), so the bound is its value at the two + /// endpoints, ordered. Its float32 evaluation t * t * (3 - 2t) is not monotone, by at most one ulp + /// (IntervalTests.MeasureMonotonicity), so the bound is widened by one ulp inside [0, 1]. + /// + public static float2 iv_smoothstep(float minimum, float maximum, float2 a) + { + var s0 = smoothstep(minimum, maximum, a.x); + var s1 = smoothstep(minimum, maximum, a.y); + var lo = min(s0, s1); + var hi = max(s0, s1); + return new(lo > 0.0f ? asfloat(asuint(lo) - 1u) : 0.0f, hi < 1.0f ? asfloat(asuint(hi) + 1u) : 1.0f); + } + + /// + /// Encloses lerp(x, y, amount) = x + (y - x) * amount for x in a, y in b, by evaluating that expression + /// with iv_sub, iv_scale and iv_add. Each step is monotone, so the result encloses the float32 lerp + /// with no widening, for any amount. It is the natural interval extension: x appears twice, so for + /// non-point a the width is (width(a) + width(b)) * |amount| + width(a), wider than the true range. + /// + public static float2 iv_lerp(float2 a, float2 b, float amount) => iv_add(a, iv_scale(iv_sub(b, a), amount)); + + /// + /// Encloses snoise(x) for every x within distance radius of centre: [n - L r, n + L r] intersected + /// with [-1, 1], where n = snoise(centre) and L = SNOISE_LIPSCHITZ. radius must be at least 0. One + /// snoise evaluation, value only. A domain warp that moves x by at most D is enclosed by passing + /// radius + D (NoiseBoundTests.WarpedPointsStayEnclosed). + /// + public static float2 iv_snoise_ball(float3 centre, float radius) + { + var n = snoise(centre); + var e = SNOISE_LIPSCHITZ * radius; + return new(max(n - e, -1.0f), min(n + e, 1.0f)); + } + + /// + /// Encloses fbm_grad(x, octaves, lacunarity, gain).w 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, with frequency f_i and + /// amplitude a_i compounded exactly as fbm_grad compounds them. octaves is clamped to [0, 16] like + /// fbm_grad. + /// + public static float2 iv_fbm_ball(float3 centre, float radius, int octaves, float lacunarity, float gain) + { + octaves = clamp(octaves, 0, 16); + var amplitude = 1.0f; + var frequency = 1.0f; + var sum = new float2(0.0f, 0.0f); + for (var i = 0; i < octaves; i++) + { + sum = iv_add(sum, iv_scale(iv_snoise_ball(centre * frequency, radius * abs(frequency)), amplitude)); + frequency *= lacunarity; + amplitude *= gain; + } + + return sum; + } +} diff --git a/packages/cultmath/src/CultMath/math.Interval.cs.meta b/packages/cultmath/src/CultMath/math.Interval.cs.meta new file mode 100644 index 00000000..eb3b76d1 --- /dev/null +++ b/packages/cultmath/src/CultMath/math.Interval.cs.meta @@ -0,0 +1,2 @@ +fileFormatVersion: 2 +guid: 3eb757caae684b388921e28e3c87ba3a diff --git a/packages/cultmath/src/CultMath/math.cs b/packages/cultmath/src/CultMath/math.cs index 23916b6b..d6c85847 100644 --- a/packages/cultmath/src/CultMath/math.cs +++ b/packages/cultmath/src/CultMath/math.cs @@ -388,6 +388,7 @@ public static partial class math public static float hash(float3 value) => hash(dot(value, new float3(127.1f, 311.7f, 74.7f))); public static uint asuint(float value) => (uint)BitConverter.SingleToInt32Bits(value); + public static float asfloat(uint value) => BitConverter.Int32BitsToSingle((int)value); // PCG hashes from Jarzynski and Olano, "Hash Functions for GPU Rendering" (JCGT 9(3), 2020). Integer-only, so // dxc and C# agree bit for bit. pcg3d/pcg4d carry uint bit patterns in int3/int4; float overloads hash IEEE bits. @@ -599,6 +600,8 @@ public static partial class math // src/noise3D.glsl at 22434e04d7 (kernel radius^2 0.5, scale 105, permute (34x+10)x); see // THIRD-PARTY-NOTICES.md. Continuous across simplex-cell boundaries // (NoiseGradTests.SnoiseAndItsGradientAreContinuousAcrossSimplexCellBoundaries). + // SNOISE_LIPSCHITZ (math.Interval.cs) is measured from this kernel: a change here must re-run + // NoiseBoundTests.MeasureLipschitz and re-pin it. public static float snoise(float3 value) { const float cx = 1.0f / 6.0f; diff --git a/packages/cultmath/tests/CultMath.Tests/HlslMirrorTests.cs b/packages/cultmath/tests/CultMath.Tests/HlslMirrorTests.cs index 6b69a908..9e2dc9af 100644 --- a/packages/cultmath/tests/CultMath.Tests/HlslMirrorTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/HlslMirrorTests.cs @@ -65,6 +65,7 @@ public sealed class HlslMirrorTests [Theory] [InlineData("CultMath.hlsl")] [InlineData("CultMath.Phacelle.hlsl")] + [InlineData("CultMath.Interval.hlsl")] public void UnityPackageShaderIsIdenticalToTheCanonicalShader(string name) { var root = GetCultMathRoot(); diff --git a/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs b/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs index 2cb8ae49..07f026ff 100644 --- a/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs @@ -31,6 +31,8 @@ public sealed class HlslSourceCompatibilityTests Assert.NotNull(ShaderFunction(shader, "cultmath_fbm_grad", typeof(float3), typeof(int), typeof(float), typeof(float))); Assert.NotNull(ShaderFunction(shader, "cultmath_ridged_grad", typeof(float3), typeof(int), typeof(float), typeof(float))); Assert.NotNull(ShaderFunction(shader, "cultmath_phacelle", typeof(float3), typeof(float3), typeof(float), typeof(float))); + Assert.NotNull(ShaderFunction(shader, "cultmath_iv_snoise_ball", typeof(float3), typeof(float))); + Assert.NotNull(ShaderFunction(shader, "cultmath_iv_fbm_ball", typeof(float3), typeof(float), typeof(int), typeof(float), typeof(float))); } // HLSL functions carry no access modifier, so they compile as private instance methods. @@ -131,7 +133,7 @@ public sealed class HlslSourceCompatibilityTests } Assert.True(mismatches.Count == 0, $"{mismatches.Count} mismatches:{Environment.NewLine}{string.Join(Environment.NewLine, mismatches.Take(12))}"); - Assert.Equal(34, compared.Count); + Assert.Equal(52, compared.Count); } /// diff --git a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl new file mode 100644 index 00000000..ed6f83db --- /dev/null +++ b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl @@ -0,0 +1,107 @@ +#ifndef CULTMATH_INTERVAL_HLSL +#define CULTMATH_INTERVAL_HLSL + +// Interval arithmetic over float2(lo, hi), lo <= hi, both finite. Every cultmath_iv_* function returns +// an interval containing f(x) for every x in its input intervals, mirroring math.Interval.cs bit for +// bit; see that file and docs/design.md ("Intervals") for each op's enclosure and its ulp widening. +// Included by CultMath.hlsl after its own definitions of cultmath_smoothstep and cultmath_snoise. + +// Lipschitz constant of cultmath_snoise(float3); provenance in math.Interval.cs (SNOISE_LIPSCHITZ). +static const float CULTMATH_SNOISE_LIPSCHITZ = 10.099261; + +float2 cultmath_iv_point(float value) { return float2(value, value); } +float2 cultmath_iv_add(float2 a, float2 b) { return float2(a.x + b.x, a.y + b.y); } +float2 cultmath_iv_sub(float2 a, float2 b) { return float2(a.x - b.y, a.y - b.x); } +float2 cultmath_iv_neg(float2 a) { return float2(-a.y, -a.x); } + +float2 cultmath_iv_mul(float2 a, float2 b) +{ + float p0 = a.x * b.x; + float p1 = a.x * b.y; + float p2 = a.y * b.x; + float p3 = a.y * b.y; + return float2(min(min(p0, p1), min(p2, p3)), max(max(p0, p1), max(p2, p3))); +} + +float2 cultmath_iv_scale(float2 a, float s) { return s >= 0.0 ? float2(a.x * s, a.y * s) : float2(a.y * s, a.x * s); } + +float2 cultmath_iv_abs(float2 a) +{ + if (a.x >= 0.0) + return a; + if (a.y <= 0.0) + return float2(-a.y, -a.x); + return float2(0.0, max(-a.x, a.y)); +} + +float2 cultmath_iv_min(float2 a, float2 b) { return float2(min(a.x, b.x), min(a.y, b.y)); } +float2 cultmath_iv_max(float2 a, float2 b) { return float2(max(a.x, b.x), max(a.y, b.y)); } + +float2 cultmath_iv_sqr(float2 a) +{ + float l = a.x * a.x; + float h = a.y * a.y; + if (a.x >= 0.0) + return float2(l, h); + if (a.y <= 0.0) + return float2(h, l); + return float2(0.0, max(l, h)); +} + +float2 cultmath_iv_sqrt(float2 a) { return float2(sqrt(max(a.x, 0.0)), sqrt(max(a.y, 0.0))); } + +// Widened by one ulp on each side: exp is not required to be correctly rounded. +float2 cultmath_iv_exp(float2 a) +{ + float lo = exp(a.x); + float hi = exp(a.y); + return float2(lo > 0.0 ? asfloat(asuint(lo) - 1u) : 0.0, asfloat(asuint(hi) + 1u)); +} + +float2 cultmath_iv_clamp(float2 a, float minimum, float maximum) +{ + return float2(cultmath_clamp(a.x, minimum, maximum), cultmath_clamp(a.y, minimum, maximum)); +} + +float2 cultmath_iv_saturate(float2 a) { return float2(cultmath_saturate(a.x), cultmath_saturate(a.y)); } + +// The endpoints, ordered, widened by one ulp inside [0, 1]: the float32 evaluation of smoothstep is +// monotone only to within one ulp. +float2 cultmath_iv_smoothstep(float minimum, float maximum, float2 a) +{ + float s0 = cultmath_smoothstep(minimum, maximum, a.x); + float s1 = cultmath_smoothstep(minimum, maximum, a.y); + float lo = min(s0, s1); + float hi = max(s0, s1); + return float2(lo > 0.0 ? asfloat(asuint(lo) - 1u) : 0.0, hi < 1.0 ? asfloat(asuint(hi) + 1u) : 1.0); +} + +// The natural interval extension of x + (y - x) * amount, which is how cultmath_lerp evaluates. +float2 cultmath_iv_lerp(float2 a, float2 b, float amount) { return cultmath_iv_add(a, cultmath_iv_scale(cultmath_iv_sub(b, a), amount)); } + +// [n - L r, n + L r] intersected with [-1, 1], n = cultmath_snoise(centre). One snoise, value only. +float2 cultmath_iv_snoise_ball(float3 centre, float radius) +{ + float n = cultmath_snoise(centre); + float e = CULTMATH_SNOISE_LIPSCHITZ * radius; + return float2(max(n - e, -1.0), min(n + e, 1.0)); +} + +// Octave sum of cultmath_iv_snoise_ball, compounding frequency and amplitude as cultmath_fbm_grad does. +float2 cultmath_iv_fbm_ball(float3 centre, float radius, int octaves, float lacunarity, float gain) +{ + octaves = clamp(octaves, 0, 16); + float amplitude = 1.0; + float frequency = 1.0; + float2 sum = float2(0.0, 0.0); + for (int i = 0; i < octaves; i++) + { + sum = cultmath_iv_add(sum, cultmath_iv_scale(cultmath_iv_snoise_ball(centre * frequency, radius * abs(frequency)), amplitude)); + frequency *= lacunarity; + amplitude *= gain; + } + + return sum; +} + +#endif diff --git a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl.meta b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl.meta new file mode 100644 index 00000000..023d218c --- /dev/null +++ b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl.meta @@ -0,0 +1,7 @@ +fileFormatVersion: 2 +guid: 18b031192f344246ac523cea014b8f24 +ShaderIncludeImporter: + externalObjects: {} + userData: + assetBundleName: + assetBundleVariant: diff --git a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl index af2d3741..da3362a6 100644 --- a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl +++ b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl @@ -504,6 +504,7 @@ CultCellular cultmath_cellular(float3 p) } #include "CultMath.Phacelle.hlsl" +#include "CultMath.Interval.hlsl" float cultmath_value_noise(float2 position) { From cd7913bfccbf3deabdd85ca60ab46fdec49c4006 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Fri, 2 Oct 2026 19:56:30 +0200 Subject: [PATCH 02/14] CultMath intervals: enclosure harness, one-ulp shrink mutants, noise bounds, Lipschitz pin, saving 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 --- .../tests/CultMath.Tests/IntervalTests.cs | 312 +++++++++++++ .../tests/CultMath.Tests/NoiseBoundTests.cs | 427 ++++++++++++++++++ 2 files changed, 739 insertions(+) create mode 100644 packages/cultmath/tests/CultMath.Tests/IntervalTests.cs create mode 100644 packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs diff --git a/packages/cultmath/tests/CultMath.Tests/IntervalTests.cs b/packages/cultmath/tests/CultMath.Tests/IntervalTests.cs new file mode 100644 index 00000000..4f165eab --- /dev/null +++ b/packages/cultmath/tests/CultMath.Tests/IntervalTests.cs @@ -0,0 +1,312 @@ +using CultMath; +using Xunit; +using static CultMath.math; + +namespace CultMath.Tests; + +/// +/// 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 : 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; constructs one on +/// each side and kills both mutants. +/// +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 Scalars, + Func Bound, + Func Pointwise, + bool NonNegativeDomain = false); + + private static float[] None(System.Random random) => Array.Empty(); + + 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; + } + + /// Runs one op over the seeded set; returns the violations, at most a few, described. + internal static List Violations(Op op, Func? mutate = null, int seed = 0x1A7E) + { + var random = new System.Random(seed); + var violations = new List(); + 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); + + /// + /// 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. + /// + [Fact] + public void EveryOneUlpShrinkIsCaught() + { + var survivors = new List(); + 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))); + } + + /// + /// smoothstep's float32 evaluation steps down by one ulp in places (MeasureMonotonicity). Three + /// consecutive floats t0 < t1 < 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. + /// + [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); + } + + /// + /// 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. + /// + [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"); + } +} diff --git a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs new file mode 100644 index 00000000..45578ec7 --- /dev/null +++ b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs @@ -0,0 +1,427 @@ +using CultMath; +using Xunit; +using static CultMath.math; + +namespace CultMath.Tests; + +/// +/// 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. +/// +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}"); + } + } + } + + /// + /// 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. + /// + [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); + } + + /// + /// 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. + /// + [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); + + /// + /// 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. + /// + [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"); + } + + /// + /// 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. + /// + [Fact] + public void TightnessReport() + { + double Ratio(Func bound, Func 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 ---- + + /// + /// 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. + /// + 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); + } + } + + /// + /// 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. + /// + [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}"); + } +} From b8dddc90d99250a68bb6b4ff55494f3ead51d9df Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Fri, 2 Oct 2026 19:56:30 +0200 Subject: [PATCH 03/14] CultMath design.md and README: the Intervals section Representation, the op list, the one-ulp widening rule and its two cases, the Lipschitz provenance and pin, what the ball bounds guarantee and what they do not, tightness, the saving and its open fork, and the consumers. Every number is one a committed test prints. Co-Authored-By: Claude Opus 5.5 --- packages/cultmath/README.md | 10 +++- packages/cultmath/docs/design.md | 98 ++++++++++++++++++++++++++++++++ 2 files changed, 106 insertions(+), 2 deletions(-) diff --git a/packages/cultmath/README.md b/packages/cultmath/README.md index f70df445..6b493a31 100644 --- a/packages/cultmath/README.md +++ b/packages/cultmath/README.md @@ -33,8 +33,14 @@ shader semantics. - `shaders/CultMath.hlsl`, a canonical HLSL mirror include for shader-side parity. HLSL already owns `float2`, `float3`, and `float4`; the include exposes `cultmath_*` functions for shared semantics HLSL does not provide. - It includes `shaders/CultMath.Phacelle.hlsl` (MPL-2.0) from the same - directory, so ship both files. + It includes `shaders/CultMath.Phacelle.hlsl` (MPL-2.0) and + `shaders/CultMath.Interval.hlsl` from the same directory, so ship all three. +- interval arithmetic for culling empty space with a proof: `iv_*` functions + over `float2(lo, hi)` (`iv_add`, `iv_mul`, `iv_exp`, `iv_smoothstep`, ...) + and `iv_snoise_ball`/`iv_fbm_ball`, which bound `snoise` and `fbm_grad` over a + ball. Every one returns an interval containing the function's value at every + point of its input, so a raymarch can skip a segment whose `hi` is below its + cutoff; `cultmath_iv_*` in HLSL. See `docs/design.md`, "Intervals". - `Voronoi.SampleTones`, a C# batch surface that calls the Rust `cultmath-core` native kernel when `cultmath_core` is available and falls back to the managed parity path otherwise. diff --git a/packages/cultmath/docs/design.md b/packages/cultmath/docs/design.md index 0de8eb14..4b139234 100644 --- a/packages/cultmath/docs/design.md +++ b/packages/cultmath/docs/design.md @@ -341,6 +341,104 @@ Known C# friction: under `using static CultMath.math;` the constructor functions (`float3(...)`, `float3x3(...)`) hide the type names in member access, so static members need a qualified type (`CultMath.float3x3.Euler(...)`). +## Intervals + +`math.Interval.cs` and `shaders/CultMath.Interval.hlsl` (included from +`CultMath.hlsl` after the Phacelle include) carry interval arithmetic for +culling empty space with a proof. Their invariant, **intervals-enclose**: every +interval function returns a `[lo, hi]` that contains `f(x)` for every `x` in +its input region, where `f(x)` is the same function evaluated pointwise in +float32. A consumer that skips a region because `hi < cutoff` has skipped +nothing that was there. + +Representation: an interval is `float2(lo, hi)` with `lo <= hi`, both finite. +There is no empty interval; a caller that needs "provably empty" tests +`hi < cutoff`. `iv_point(x)` is `float2(x, x)`. + +Ops: `iv_point`, `iv_add`, `iv_sub`, `iv_neg`, `iv_mul` (the least and greatest +of the four corner products), `iv_scale(iv, float)`, `iv_abs` and `iv_sqr` +(an interval straddling zero starts at 0), `iv_min`, `iv_max`, `iv_sqrt` +(`lo` clamped at 0), `iv_exp`, `iv_clamp(iv, float, float)`, `iv_saturate`, +`iv_smoothstep(float, float, iv)` (its ordered endpoints) and +`iv_lerp(iv, iv, float)`. No others until a consumer names one. Each has an +HLSL mirror named `cultmath_iv_*`, compared bit for bit by +`HlslSourceCompatibilityTests`. + +The ulp-widening rule. Each bound is computed with the pointwise function's +own float32 operations. IEEE round-to-nearest is monotone, so an op that is +monotone on its region encloses its float32 evaluation exactly, with no +widening; that covers every op above except two, which widen by exactly one +ulp, in the implementation and never in the test: + +- `iv_exp`: `exp` is not required to be correctly rounded, so a libm need not + be monotone. .NET's `MathF.Exp` on Linux steps down 0 times over every float + in [-87, 88.7] (`IntervalTests.MeasureMonotonicity`), so on that CPU the + widening is a margin; it is there for other libms and GPUs. +- `iv_smoothstep`: `t * t * (3 - 2t)` in float32 steps down 337,095 times over + the floats in [0, 1], each time by at most 1 ulp (`MeasureMonotonicity`). + Since `t = saturate((x - minimum) / (maximum - minimum))` is itself monotone + in `x`, that covers every `smoothstep`. + +`iv_lerp` is the natural interval extension of `x + (y - x) * amount`, which is +how `lerp` evaluates; every step is monotone, so it encloses with no widening +for any `amount`, and it is tight only when `a` is a point, because `x` appears +twice. + +`IntervalTests` runs every op over 10,000 seeded intervals (centres in +[-20, 20], widths log-uniform in [1e-4, 10], every eighth a point) and 16 +points inside each, endpoints and zero included, with no tolerance; and for +every op it raises `lo` and lowers `hi` by one ulp and requires the harness to +catch each. A new op is one entry in that table. + +Noise bounds. `iv_snoise_ball(c, r)` is `[n - L r, n + L r]` intersected with +`[-1, 1]`, where `n = snoise(c)` and `L = SNOISE_LIPSCHITZ`: one `snoise`, value +only. No interval function reads `snoise_grad`, and none has a value-only or +gradient twin. `iv_fbm_ball(c, r, octaves, lacunarity, gain)` is the octave sum +of `iv_snoise_ball(c * f_i, r * |f_i|)` scaled by `a_i`, compounding frequency +and amplitude exactly as `fbm_grad` does, octaves clamped to [0, 16]. A domain +warp that moves a point by at most `D` is enclosed by the radius `r + D` +(`NoiseBoundTests.WarpedPointsStayEnclosed`); this is what the gamecult.org +ground's march rests on. + +`L` has provenance, not a proof. `NoiseBoundTests.MeasureLipschitz` (slow, +explicit) takes the largest `|snoise_grad|` over 1e6 seeded points in +[-256, 256]^3 (7.1933103), refines the 1e4 largest by gradient ascent +(9.181147), and multiplies by 1.10: `L = 10.099261`. It also reaches +`|snoise| = 0.9718335` by the same ascent on the value, which is what lets the +bound intersect with `[-1, 1]`. `LipschitzConstantPinsSampledGradients` +checks 1e5 more seeded points against `L` and climbs again from the start +`MeasureLipschitz` prints as its witness, failing if `L` is not 1.10 times that +maximum to within half a percent, so a 1% move of `L` fails. The enclosure tests +over 2,000 balls x 64 points are the defence. A change to the `snoise` kernel +must re-run `MeasureLipschitz` and re-pin `L` and the witness. + +What the ball bounds guarantee: enclosure of `snoise` and `fbm_grad(...).w` over +the ball, in float32. What they do not: anything about tapes, pruning, choice +tracking or affine forms; those are later steps of +`docs/cultmath-tape-target.md` (CultLib root). How loose they are is a number: +`NoiseBoundTests.TightnessReport` prints the mean interval width over the +sampled range: 5.08 for `iv_snoise_ball` at `r = 0.05`, 1.96 at 0.25 and 1.28 +at 1.0; 6.64 and 2.82 for 4-octave `iv_fbm_ball` at 0.05 and 0.25. + +The saving is a number too. `NoiseBoundTests.IntervalSkipHalvesEvaluations` +marches 1,000 seeded rays through the map's warped coverage field (`F0 = 0.35`, +cutoff 0.25, a slab of depth 12) once at the dense step and once with the +adaptive interval march, and counts `snoise` evaluations (three per dense +sample, one per probe; the flow is not counted, because the site reads it from +a texture). Both marches integrate the same cells, so their transmittance +agrees exactly. At the test's parameters (first step 0.25, dense step 0.125, +warp bound 0.05) the dense march spends 287,940 evaluations and the interval +march 309,550: 0.93x, not the 2x saving the proving-ground contract asks for. +Over the grid the test also prints, the best is 1.83x (first step 0.0625, no +warp), and at a first step of 1 it is 0.74x. With `L` near 10 a probe proves a +segment empty only where `snoise` is far below the cutoff, and a probe that +fails costs one evaluation on top of the three the dense sample costs. The test +skips with that open fork named instead of passing; whether the march, the +field or the bound changes is the campaign's decision. + +Consumers: the gamecult.org ground shader (through the GLSL lowering) and +Aetheria's nebula raymarch at its CultMath pin bump. + ## Rules - Keep type names HLSL-shaped: `float2`, `float3`, `float4`, `float3x3`, and `math`. From c9bd003bb37712755173f6c8773b12c1f1ad6552 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Fri, 2 Oct 2026 19:57:01 +0200 Subject: [PATCH 04/14] CultMath intervals: comments name the gradient without the snoise_grad/fbm_grad tokens The cut's negative check is a pattern over math.Interval.cs and CultMath.Interval.hlsl: the bounds use values only, so neither file names the gradient functions. Comments only. Co-Authored-By: Claude Opus 5.5 --- .../cultmath/shaders/CultMath.Interval.hlsl | 2 +- .../cultmath/src/CultMath/math.Interval.cs | 18 ++++++++++-------- .../Shaders/CultMath.Interval.hlsl | 2 +- 3 files changed, 12 insertions(+), 10 deletions(-) diff --git a/packages/cultmath/shaders/CultMath.Interval.hlsl b/packages/cultmath/shaders/CultMath.Interval.hlsl index ed6f83db..3fdac41d 100644 --- a/packages/cultmath/shaders/CultMath.Interval.hlsl +++ b/packages/cultmath/shaders/CultMath.Interval.hlsl @@ -87,7 +87,7 @@ float2 cultmath_iv_snoise_ball(float3 centre, float radius) return float2(max(n - e, -1.0), min(n + e, 1.0)); } -// Octave sum of cultmath_iv_snoise_ball, compounding frequency and amplitude as cultmath_fbm_grad does. +// Octave sum of cultmath_iv_snoise_ball, compounding frequency and amplitude as the value-and-gradient fBm does. float2 cultmath_iv_fbm_ball(float3 centre, float radius, int octaves, float lacunarity, float gain) { octaves = clamp(octaves, 0, 16); diff --git a/packages/cultmath/src/CultMath/math.Interval.cs b/packages/cultmath/src/CultMath/math.Interval.cs index 1430e385..7c2ba6c9 100644 --- a/packages/cultmath/src/CultMath/math.Interval.cs +++ b/packages/cultmath/src/CultMath/math.Interval.cs @@ -11,10 +11,11 @@ public static partial class math // the libm is not required to be, the bound is widened by exactly one ulp, and the doc comment // says so. There is no empty interval: a caller that needs "provably empty" tests hi < cutoff. - // Lipschitz constant of snoise(float3) in its input: an upper bound on |snoise_grad(p).xyz|. Provenance: - // NoiseBoundTests.MeasureLipschitz (slow, explicit) takes the largest |snoise_grad| over 1e6 seeded - // points, refines the 1e4 largest by gradient ascent, and multiplies the refined maximum by 1.10; - // NoiseBoundTests.LipschitzConstantPinsSampledGradients pins this value against a fast re-measure. + // Lipschitz constant of snoise(float3) in its input: an upper bound on the length of its analytic + // gradient. Provenance: NoiseBoundTests.MeasureLipschitz (slow, explicit) takes the largest gradient + // length over 1e6 seeded points, refines the 1e4 largest by gradient ascent, and multiplies the + // refined maximum by 1.10; NoiseBoundTests.LipschitzConstantPinsSampledGradients pins this value + // against a fast re-measure and the ascent from MeasureLipschitz's witness. // Empirical with a margin, not a proof; the enclosure tests are the defence. A change to the // snoise kernel must re-run MeasureLipschitz and re-pin this constant. public const float SNOISE_LIPSCHITZ = 10.099261f; @@ -142,10 +143,11 @@ public static partial class math } /// - /// Encloses fbm_grad(x, octaves, lacunarity, gain).w 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, with frequency f_i and - /// amplitude a_i compounded exactly as fbm_grad compounds them. octaves is clamped to [0, 16] like - /// fbm_grad. + /// 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, + /// with frequency f_i and amplitude a_i compounded exactly as the value-and-gradient fBm compounds + /// them (its .w is what NoiseBoundTests.FbmBallEnclosesPoints encloses). octaves is clamped to [0, 16] + /// as there. /// public static float2 iv_fbm_ball(float3 centre, float radius, int octaves, float lacunarity, float gain) { diff --git a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl index ed6f83db..3fdac41d 100644 --- a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl +++ b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl @@ -87,7 +87,7 @@ float2 cultmath_iv_snoise_ball(float3 centre, float radius) return float2(max(n - e, -1.0), min(n + e, 1.0)); } -// Octave sum of cultmath_iv_snoise_ball, compounding frequency and amplitude as cultmath_fbm_grad does. +// Octave sum of cultmath_iv_snoise_ball, compounding frequency and amplitude as the value-and-gradient fBm does. float2 cultmath_iv_fbm_ball(float3 centre, float radius, int octaves, float lacunarity, float gain) { octaves = clamp(octaves, 0, 16); From 9b9a8b3e050653530939cfcb8d5ac4429fd8513f Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Sat, 3 Oct 2026 00:04:44 +0200 Subject: [PATCH 05/14] CultMath intervals r2: iv_frustum_ball, the ball of a screen tile's rays over a depth segment One new function in C# and both tracked HLSL copies (one blob): the centre on the tile's central ray at mid-depth, radius ((z1 - z0) / 2) |(m_c, 1)| + z1 N / (sqrt(2) f) + warp. TileBallEnclosesEveryRaySegment checks 2,000 seeded tiles x 64 points (sub-pixel corners at both depth ends with full-length flows) with no tolerance, and pins the ball to the derivation so a grown ball fails too. The mirror's compared-function count moves from 52 to 53. Co-Authored-By: Claude Opus 5.5 --- .../cultmath/shaders/CultMath.Interval.hlsl | 11 ++++ .../cultmath/src/CultMath/math.Interval.cs | 23 +++++++ .../HlslSourceCompatibilityTests.cs | 2 +- .../tests/CultMath.Tests/NoiseBoundTests.cs | 60 ++++++++++++++++++- .../Shaders/CultMath.Interval.hlsl | 11 ++++ 5 files changed, 104 insertions(+), 3 deletions(-) diff --git a/packages/cultmath/shaders/CultMath.Interval.hlsl b/packages/cultmath/shaders/CultMath.Interval.hlsl index 3fdac41d..732ef3a0 100644 --- a/packages/cultmath/shaders/CultMath.Interval.hlsl +++ b/packages/cultmath/shaders/CultMath.Interval.hlsl @@ -104,4 +104,15 @@ float2 cultmath_iv_fbm_ball(float3 centre, float radius, int octaves, float lacu return sum; } +// A ball (xyz centre, w radius) enclosing every point of every ray of a screen tile over the depth segment +// [z0, z1], warped by at most warp. Camera frame: pinhole at the origin looking down +z, ray (m z, z). +// footprintPerDepth = N / (sqrt(2) f) for an N x N tile at focal length f. Derivation in math.Interval.cs. +float4 cultmath_iv_frustum_ball(float2 centreSlope, float z0, float z1, float footprintPerDepth, float warp) +{ + float zm = (z0 + z1) * 0.5; + float axis = sqrt(centreSlope.x * centreSlope.x + centreSlope.y * centreSlope.y + 1.0); + float radius = (z1 - z0) * 0.5 * axis + z1 * footprintPerDepth + warp; + return float4(centreSlope.x * zm, centreSlope.y * zm, zm, radius); +} + #endif diff --git a/packages/cultmath/src/CultMath/math.Interval.cs b/packages/cultmath/src/CultMath/math.Interval.cs index 7c2ba6c9..547c4acc 100644 --- a/packages/cultmath/src/CultMath/math.Interval.cs +++ b/packages/cultmath/src/CultMath/math.Interval.cs @@ -164,4 +164,27 @@ public static partial class math return sum; } + + /// + /// A ball (xyz centre, w radius) enclosing every point of every ray of a screen tile over the depth + /// segment [z0, z1], each point then moved by a warp of length at most warp: the ball to hand + /// iv_snoise_ball (scaled by the noise frequency) when one probe serves the whole tile. Camera frame: + /// a pinhole at the origin looking down +z, where the ray of slope m is the points (m z, z). + /// centreSlope is the tile's central slope m_c; footprintPerDepth is the largest distance of any + /// ray's slope from m_c, N / (sqrt(2) f) for an N x N pixel tile at focal length f pixels (the full + /// pixel footprint, so sub-pixel jitter is covered). The centre is (m_c z_m, z_m), z_m = (z0 + z1) / 2; + /// the radius is ((z1 - z0) / 2) |(m_c, 1)| + z1 footprintPerDepth + warp. Derivation: a point of slope + /// m_c + e (|e| <= footprintPerDepth) at depth z is c + (z - z_m)(m_c, 1) + z (e, 0); the first term is + /// at most the first radius term, the second at most the second, and the warp adds at most warp + /// (NoiseBoundTests.TileBallEnclosesEveryRaySegment). The two vectors are never parallel, so the + /// triangle inequality is strict. The caller rotates and translates the centre into world space (the + /// radius is unchanged) and scales both by its noise frequency. + /// + public static float4 iv_frustum_ball(float2 centreSlope, float z0, float z1, float footprintPerDepth, float warp) + { + var zm = (z0 + z1) * 0.5f; + var axis = sqrt(centreSlope.x * centreSlope.x + centreSlope.y * centreSlope.y + 1.0f); + var radius = (z1 - z0) * 0.5f * axis + z1 * footprintPerDepth + warp; + return new float4(centreSlope.x * zm, centreSlope.y * zm, zm, radius); + } } diff --git a/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs b/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs index 07f026ff..a5965a3d 100644 --- a/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs @@ -133,7 +133,7 @@ public sealed class HlslSourceCompatibilityTests } Assert.True(mismatches.Count == 0, $"{mismatches.Count} mismatches:{Environment.NewLine}{string.Join(Environment.NewLine, mismatches.Take(12))}"); - Assert.Equal(52, compared.Count); + Assert.Equal(53, compared.Count); } /// diff --git a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs index 45578ec7..0235b322 100644 --- a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs @@ -7,8 +7,12 @@ namespace CultMath.Tests; /// /// 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. +/// measured provenance, iv_frustum_ball encloses a screen tile's rays over a depth segment, 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. +/// +/// Hand mutations of iv_frustum_ball that TileBallEnclosesEveryRaySegment must kill: z0 for z1 in the +/// lateral term (z0 * footprintPerDepth); |m_c| for |(m_c, 1)| in the depth term; the warp dropped. /// public sealed class NoiseBoundTests { @@ -107,6 +111,57 @@ public sealed class NoiseBoundTests } } + // ---- The tile ball and the composed bound ---- + + /// + /// iv_frustum_ball encloses every point of every ray of an N x N tile over a depth segment, sub-pixel + /// jitter included, each point then moved by a flow of length at most warp, with no tolerance: 2,000 + /// seeded tiles (N in {1, 4, 8, 16}, f in [500, 3000], |m_c| up to 1.1, segments over the slab's, the + /// fog's and the void's depth ranges, a quarter starting at the camera) x 64 points, the first eight at + /// a corner of the footprint at either end of the segment with a full-length flow. It also pins the + /// ball to the derivation: centre and radius equal the formula evaluated in double to one part in + /// 1e5, so a ball that grows (sound but wasteful) fails too. Each tile's camera is drawn on its own, as + /// a moving camera's frames are: no ball, mask or probe result is carried from one to the next. + /// + [Fact] + public void TileBallEnclosesEveryRaySegment() + { + var random = new System.Random(0x711E); + int[] sizes = { 1, 4, 8, 16 }; + float[] depths = { 12.0f, 1000.0f, 2400.0f }; + for (var t = 0; t < 2000; t++) + { + var n = sizes[random.Next(sizes.Length)]; + var focal = Uniform(random, 500.0f, 3000.0f); + var slope = new float2(Uniform(random, -0.77f, 0.77f), Uniform(random, -0.77f, 0.77f)); + var depth = depths[random.Next(depths.Length)]; + var z0 = random.Next(4) == 0 ? 0.0f : Uniform(random, 0.0f, depth * 0.9f); + var z1 = z0 + LogUniform(random, 1.0e-3f * depth, depth - z0); + var warp = random.Next(4) == 0 ? 0.0f : LogUniform(random, 1.0e-3f, 60.0f); + var half = n / (2.0f * focal); + var footprint = n / (MathF.Sqrt(2.0f) * focal); + var ball = iv_frustum_ball(slope, z0, z1, footprint, warp); + var centre = new float3(ball.x, ball.y, ball.z); + + var zm = (z0 + (double)z1) / 2.0; + var radius = (z1 - (double)z0) / 2.0 * Math.Sqrt((double)slope.x * slope.x + (double)slope.y * slope.y + 1.0) + z1 * (double)footprint + warp; + Assert.True(Math.Abs(ball.w - radius) <= 1.0e-5 * radius, $"radius {ball.w:R}, derivation {radius:R}"); + Assert.True(Math.Abs(ball.x - slope.x * zm) + Math.Abs(ball.y - slope.y * zm) + Math.Abs(ball.z - zm) <= 1.0e-5 * zm, $"centre {centre}, derivation ({slope.x * zm:R}, {slope.y * zm:R}, {zm:R})"); + + for (var i = 0; i < 64; 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 flow = UnitVector(random) * (i < 8 ? warp : warp * MathF.Sqrt(random.NextSingle())); + var p = new float3(m.x * z, m.y * z, z) + flow; + Assert.True(length(p - centre) <= ball.w, $"iv_frustum_ball({slope}, {z0:R}, {z1:R}, {footprint:R}, {warp:R}) = {ball} misses {p} (slope {m}, depth {z:R}), distance {length(p - centre):R}"); + } + } + } + // ---- The Lipschitz constant ---- private static float GradientNorm(float3 p) @@ -281,6 +336,7 @@ public sealed class NoiseBoundTests Assert.All(snoiseRatios.Concat(fbmRatios), t => Assert.True(double.IsFinite(t.Item2))); } + // ---- The saving ---- /// diff --git a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl index 3fdac41d..732ef3a0 100644 --- a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl +++ b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl @@ -104,4 +104,15 @@ float2 cultmath_iv_fbm_ball(float3 centre, float radius, int octaves, float lacu return sum; } +// A ball (xyz centre, w radius) enclosing every point of every ray of a screen tile over the depth segment +// [z0, z1], warped by at most warp. Camera frame: pinhole at the origin looking down +z, ray (m z, z). +// footprintPerDepth = N / (sqrt(2) f) for an N x N tile at focal length f. Derivation in math.Interval.cs. +float4 cultmath_iv_frustum_ball(float2 centreSlope, float z0, float z1, float footprintPerDepth, float warp) +{ + float zm = (z0 + z1) * 0.5; + float axis = sqrt(centreSlope.x * centreSlope.x + centreSlope.y * centreSlope.y + 1.0); + float radius = (z1 - z0) * 0.5 * axis + z1 * footprintPerDepth + warp; + return float4(centreSlope.x * zm, centreSlope.y * zm, zm, radius); +} + #endif From 6f2f3a19a38e82a768654fc2c94391c4cbd72602 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Sat, 3 Oct 2026 00:04:54 +0200 Subject: [PATCH 06/14] CultMath intervals r2: cut r1's per-ray march; envelope-then-noise fields and their composed bound Deleted: SavingField's MarchIntervals (per-ray probes, the StepMax growth cap, probes attempted where the ball could never prove emptiness), its unconditional three-snoise sample, and the saving test that skipped naming r1's fork. Added, test-side only: the scenario model. FogField is Aetheria's nebula over its gravity-well height map (Volumetric.cginc, Zone.cs, Settings.asset cited by line): (a) height fog, (b) inside the fog, (c) the zone bowl and four wells. SlabField is r1's slab, kept as the worst case. VoidField is the shipped cavity with footprint octave weights. Each composes one bound from existing iv_* ops: the envelope over the tile slice's box, its fade, and the noise ball from iv_frustum_ball enlarged by the warp. The pointwise density is envelope-gated. EnvelopeBoundEnclosesDensity checks 2,000 slices x 64 warped points with no tolerance. Co-Authored-By: Claude Opus 5.5 --- .../tests/CultMath.Tests/NoiseBoundTests.cs | 753 +++++++++++++++--- 1 file changed, 634 insertions(+), 119 deletions(-) diff --git a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs index 0235b322..74edd778 100644 --- a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs @@ -7,9 +7,10 @@ namespace CultMath.Tests; /// /// 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, iv_frustum_ball encloses a screen tile's rays over a depth segment, 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. +/// measured provenance, iv_frustum_ball encloses a screen tile's rays over a depth segment, the bound +/// composed from an analytic envelope and the noise ball encloses the density, and a tile march over +/// fields with authored empty space is measured against a dense one for the same transmittance. Lines +/// starting IV-REPORT carry the numbers the cut report quotes. /// /// Hand mutations of iv_frustum_ball that TileBallEnclosesEveryRaySegment must kill: z0 for z1 in the /// lateral term (z0 * footprintPerDepth); |m_c| for |(m_c, 1)| in the depth term; the warp dropped. @@ -162,6 +163,49 @@ public sealed class NoiseBoundTests } } + /// + /// intervals-enclose for the composition the site and Aetheria use: over 2,000 seeded tile slices + /// (Aetheria's wells, inside the fog, the void at every grid point; N in {1, 4, 8, 16}; a camera drawn + /// per tile) the bound composed from the envelope over the slice's box, its fade, and the noise ball + /// enlarged by the warp encloses the density at 64 points of the slice, each sampled through the field's + /// own warp, with no tolerance. The void's footprint-truncated density (octave weights) is checked + /// against the same bound. + /// + [Fact] + public void EnvelopeBoundEnclosesDensity() + { + var random = new System.Random(0xE7B0); + int[] sizes = { 1, 4, 8, 16 }; + var fields = new List { FogField.Wells(FogField.AetheriaWarp), FogField.InsideFog(FogField.AetheriaWarp) }; + foreach (var rh in VoidField.HollowRadii) + foreach (var offset in VoidField.CameraOffsets) + foreach (var ramp in VoidField.Ramps) + fields.Add(new VoidField(rh, ramp, offset.Fraction, offset.Name)); + var counts = new Counts(); + for (var t = 0; t < 2000; t++) + { + var field = t % 2 == 0 ? fields[(t / 2) % 2] : fields[2 + random.Next(fields.Count - 2)]; + var tile = field.DrawTile(random, sizes[random.Next(sizes.Length)]); + var near = field.Grid[0]; + var far = field.Grid[^1]; + var z0 = Uniform(random, near, far * 0.9f); + var z1 = z0 + LogUniform(random, 1.0e-3f * far, far - z0); + var bound = field.Bound(tile, z0, z1, gated: false, counts); + for (var i = 0; i < 64; i++) + { + var m = tile.PixelSlope(random.Next(tile.N), random.Next(tile.N), random.NextSingle(), random.NextSingle()); + var z = i < 2 ? (i == 0 ? z0 : z1) : Uniform(random, z0, z1); + var density = field.Density(tile, m, z, lod: false, counts); + Assert.True(Inside(density, bound), $"{field.Name}: bound {bound} over [{z0:R}, {z1:R}] misses density {density:R} at slope {m}, depth {z:R}"); + if (field is VoidField) + { + var truncated = field.Density(tile, m, z, lod: true, counts); + Assert.True(Inside(truncated, bound), $"{field.Name}: bound {bound} misses the truncated density {truncated:R} at slope {m}, depth {z:R}"); + } + } + } + } + // ---- The Lipschitz constant ---- private static float GradientNorm(float3 p) @@ -337,147 +381,618 @@ public sealed class NoiseBoundTests } - // ---- The saving ---- + // ---- The saving: envelope-then-noise fields, marched over screen tiles ---- + + // Combined cost in snoise units: one envelope evaluation, pointwise or as an interval probe, costs this + // much. An estimate (a handful of pulses or one height-map fetch against ~150 flops of snoise); both + // counts are printed, so any other weight can be applied to them. + internal const double EnvelopeCost = 0.2; + + internal sealed class Counts + { + public long Snoise; + public long Envelope; + + public double Cost => Snoise + EnvelopeCost * Envelope; + } /// - /// 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. + /// One screen tile of one frame: a pinhole camera (origin and orthonormal basis; camera space looks + /// down +z and the ray of slope m is the points (m z, z)), the tile's central slope, its size N, the + /// focal length in pixels and the frame time. Every tile draws its own camera, so a probe result is + /// never carried to another frame. /// - internal sealed class SavingField + internal class Tile + { + public float3 Origin; + public float3 Right = new(1.0f, 0.0f, 0.0f); + public float3 Up = new(0.0f, 1.0f, 0.0f); + public float3 Forward = new(0.0f, 0.0f, 1.0f); + public float2 Slope; + public int N; + public float Focal; + public float Time; + + public float Half => N / (2.0f * Focal); + public float Footprint => N / (MathF.Sqrt(2.0f) * Focal); + public float2 SlopesX => new(Slope.x - Half, Slope.x + Half); + public float2 SlopesY => new(Slope.y - Half, Slope.y + Half); + + // Pixel (i, j) of the tile at sub-pixel offset (u, v) in [0, 1): its slope lies in SlopesX x SlopesY + // in float32, because (i + u) / f - N / (2f) is at most N / f - N / (2f) = N / (2f) exactly. + public float2 PixelSlope(int i, int j, float u, float v) => + Slope + new float2((i + u) / Focal - Half, (j + v) / Focal - Half); + + public float3 World(float3 c) => Origin + Right * c.x + Up * c.y + Forward * c.z; + + public float3 Camera(float3 w) + { + var d = w - Origin; + return new float3(dot(d, Right), dot(d, Up), dot(d, Forward)); + } + + public void Look(float3 forward) + { + Forward = normalize(forward); + Right = normalize(cross(new float3(0.0f, 1.0f, 0.0f), Forward)); + Up = cross(Forward, Right); + } + } + + /// + /// A field with an analytic envelope and the dense grid it is marched on. Density is the pointwise + /// field, envelope-gated (noise only where the envelope leaves the density undecided), counted. Bound + /// is the one composition of the envelope's interval over a tile slice's box and the noise ball over the + /// slice's iv_frustum_ball: an interval enclosing Density at every point of every ray of the tile in the + /// slice. A slice is provably empty iff Bound(...).y <= 0. gated is the march's form: it returns the + /// envelope's verdict without snoise when the envelope alone proves the slice empty, and attempts the + /// noise ball only when the best outcome any centre value could give would prove it; ungated, the ball + /// is always evaluated. + /// + internal abstract class Scenario + { + public required string Name { get; init; } + public float[] Grid { get; protected init; } = Array.Empty(); + public float Focal { get; init; } = 935.0f; // 1080 rows at a 60 degree vertical field of view + public float Extinction { get; init; } = 1.0f; + + public abstract Tile DrawTile(System.Random random, int n); + + public abstract float Density(Tile tile, float2 slope, float z, bool lod, Counts counts); + + public abstract float2 Bound(Tile tile, float z0, float z1, bool gated, Counts counts); + + // Cell boundaries near + (j / J)^2 (far - near): Aetheria's quadratic sample spacing + // (Assets/Shaders/Raymarching/CloudShader.shader:116). + protected static float[] Quadratic(float near, float far, int cells) => + Enumerable.Range(0, cells + 1).Select(j => near + (float)j / cells * ((float)j / cells) * (far - near)).ToArray(); + + protected static float[] Cells(float step, float end) => + Enumerable.Range(0, (int)MathF.Ceiling(end / step) + 1).Select(j => j * step).ToArray(); + } + + // iv_snoise_ball, with the snoise skipped where L r >= 2: there the bound is [-1, 1] whatever snoise(c) + // is (|snoise| < 1), so the result is the same and costs nothing. + private static float2 NoiseBall(float3 centre, float radius, Counts counts) + { + if (SNOISE_LIPSCHITZ * radius >= 2.0f) + return new float2(-1.0f, 1.0f); + counts.Snoise++; + return iv_snoise_ball(centre, radius); + } + + // The best interval any centre value could give a ball of this radius, toward low or high values; the + // gate on whether a probe can prove anything. snoise(c) > -1, so a real ball's hi is at least -1 + L r, + // and likewise its lo is at most 1 - L r. + private static float2 BestBall(float radius, bool low) + { + var e = SNOISE_LIPSCHITZ * radius; + if (e >= 2.0f) + return new float2(-1.0f, 1.0f); + return low ? new float2(-1.0f, min(-1.0f + e, 1.0f)) : new float2(max(1.0f - e, -1.0f), 1.0f); + } + + // Least and greatest distance along one axis from a box to a point, exact in float32: a coordinate in + // [lo, hi] is at least `near` and at most `far` from c, and squares, sums and sqrt are monotone. + private static (float Near, float Far) AxisDistance(float2 range, float c) => + (c < range.x ? range.x - c : c > range.y ? c - range.y : 0.0f, max(abs(range.x - c), abs(range.y - c))); + + /// + /// The noise of the fog and the void (r1's field shape): coverage is two snoise of the warped point, + /// cross-faded over two half-period phases; n = coverage + detail / 2, the detail one snoise at four + /// times the frequency. The flow is a unit-bounded 2D field times FlowAmplitude, so a phase moves a point + /// by at most Warp = FlowAmplitude * Period / 2. The flow is not counted: the site reads it from a + /// texture and Aetheria from its flow map. + /// + internal sealed class WarpedNoise + { + public required float F0 { get; init; } + public required float Warp { get; init; } + public float Period { get; init; } = 8.0f; + public float FlowFrequency { get; init; } = 1.0f / 512.0f; + public bool FlowInXY { get; init; } + + private float FlowAmplitude => Warp * 2.0f / Period; + + private float3 Flow(float3 p) + { + var q = (FlowInXY ? new float2(p.x, p.y) : new float2(p.x, p.z)) * FlowFrequency; + var v = new float2(snoise(q + new float2(17.0f, 3.0f)), snoise(q + new float2(-5.0f, 41.0f))); + v = v / max(1.0f, length(v)) * FlowAmplitude; + return FlowInXY ? new float3(v.x, v.y, 0.0f) : new float3(v.x, 0.0f, v.y); + } + + private (float W0, float Shift0, float Shift1) Phases(float time) + { + var phase0 = frac(time / Period); + var phase1 = frac(time / Period + 0.5f); + return (1.0f - abs(2.0f * phase0 - 1.0f), (phase0 - 0.5f) * Period, (phase1 - 0.5f) * Period); + } + + public float Coverage(float3 p, float time, Counts counts) + { + var (w0, shift0, shift1) = Phases(time); + var flow = Flow(p); + counts.Snoise += 2; + return w0 * snoise((p + flow * shift0) * F0) + (1.0f - w0) * snoise((p + flow * shift1) * F0); + } + + public float Detail(float3 p, Counts counts) + { + counts.Snoise++; + return snoise(p * (4.0f * F0) + 17.0f); + } + + // Both phases' warped points lie in a ball that carries the warp, so one ball encloses both snoise; + // the cross-fade is enclosed term by term with the frame's own weights. + public float2 CoverageBound(float2 ball, float time) + { + var w0 = Phases(time).W0; + return iv_add(iv_scale(ball, w0), iv_scale(ball, 1.0f - w0)); + } + + public float2 CoverageBall(float3 centre, float radius, Counts counts) => NoiseBall(centre * F0, radius * F0, counts); + + public float2 DetailBall(float3 centre, float radius, Counts counts) => NoiseBall(centre * (4.0f * F0) + 17.0f, radius * (4.0f * F0), counts); + + public float CoverageRadius(float radius) => radius * F0; + + public float DetailRadius(float radius) => radius * (4.0f * F0); + } + + /// + /// Aetheria's nebula (Assets/Shaders/Volumetric.cginc:165-221, cloudDensity :91-102) over its zone + /// height map (Assets/Scripts/ServerShared/Zone.cs:369-397 GetHeight, :426 PowerPulse), in its units + /// (Assets/Resources/Settings.asset: ZoneDepth :97-98, DefaultEnvironment :386-411). The pre-distortion + /// SDF is s(p) = p.y + h(p.xz), h >= 0 the well map: the zone bowl PowerPulse(|xz| / 2R, 2) * 64 plus, + /// per well, PowerPulse(|xz - c_b| / r_b, 16) * d_b, with r_b = 500 M^0.25, d_b = 30 M^0.175 and + /// PowerPulse(x, e) = (1 - 4x^2)^e on [0, 1/2], 0 beyond. Noise displaces only below the safety distance: + /// fade = 1 - smoothstep(0.75 S, S, s), s' = s + A fade n(warp(p)), density = max(0, (F - s') / B); the + /// patch band and the 1e-9 fill are left out. The warp moves only the noise argument; the height is read + /// at the unwarped xz, as Aetheria reads _NebulaSurfaceHeight at pos.xz. Envelope coordinates are the + /// tile's camera frame (the camera yaws about y, so horizontal distances are the world's). Extinction + /// 0.5 is the environment's (Settings.asset:397, applied at CloudShader.shader:91). + /// + internal sealed class FogField : Scenario + { + public const float Safety = 30.0f; // SafetyDistance, Settings.asset:399 + public const float FloorOffset = -20.0f; // :393 + public const float FloorBlend = 10.0f; // :394 + public const float Amplitude = 20.0f; + public const float NoiseScale = 414.2167f; // :409 + public const float AetheriaWarp = 60.0f; // Flow GlobalAmplitude 15 x Period 8 / 2, :402-405 + public const float ZoneRadius = 2000.0f; + private static readonly float2 FreeNoise = new(-1.5f, 1.5f); + + public required WarpedNoise Noise { get; init; } + public (float2 Centre, float Scale, float Exponent, float Depth)[] Bowls { get; init; } = Array.Empty<(float2, float, float, float)>(); + public required Func CameraAt { get; init; } + public required float2 SlopesY { get; init; } + + public FogField() + { + // Near 0.3 and far 1000: the ARPG scene's camera (Assets/Scenes/ARPG.unity:635-636). + Grid = Quadratic(0.3f, 1000.0f, 256); + Extinction = 0.5f; + } + + private static WarpedNoise AetheriaNoise(float warp) => new() { F0 = 1.0f / NoiseScale, Warp = warp }; + + /// (a) Height fog: no wells, the camera above the safety band at y = S + 30, rays level and up. + public static FogField HeightFog(float warp) => new() + { + Name = "(a) height fog", Noise = AetheriaNoise(warp), SlopesY = new float2(0.0f, 0.5f), + CameraAt = (random, _) => new float3(Uniform(random, -1000.0f, 1000.0f), Safety + 30.0f, Uniform(random, -1000.0f, 1000.0f)), + }; + + /// (b) Inside the fog: no wells, the camera at y = F - 10, rays level. + public static FogField InsideFog(float warp) => new() + { + Name = "(b) inside the fog", Noise = AetheriaNoise(warp), SlopesY = new float2(-0.05f, 0.05f), + CameraAt = (random, _) => new float3(Uniform(random, -1000.0f, 1000.0f), FloorOffset - 10.0f, Uniform(random, -1000.0f, 1000.0f)), + }; + + /// + /// (c) Aetheria-like: the zone bowl R = 2000 and four wells of mass 100, 1000, 10000 and 1000 at seeded + /// positions inside it; the camera at a seeded xz in the bowl, above the noise band at + /// y = max(0, S + 10 - h), gazing across with m_y in [-0.25, 0.1]. + /// + public static FogField Wells(float warp) + { + var random = new System.Random(0xA37E); + var bowls = new List<(float2, float, float, float)> { (new float2(0.0f, 0.0f), 2.0f * ZoneRadius, 2.0f, 64.0f) }; + foreach (var mass in new[] { 100.0f, 1000.0f, 10000.0f, 1000.0f }) + { + var angle = Uniform(random, 0.0f, 2.0f * MathF.PI); + var radius = ZoneRadius * MathF.Sqrt(random.NextSingle()); + bowls.Add((new float2(radius * MathF.Cos(angle), radius * MathF.Sin(angle)), 500.0f * MathF.Pow(mass, 0.25f), 16.0f, 30.0f * MathF.Pow(mass, 0.175f))); + } + + return new FogField + { + Name = "(c) Aetheria wells", Noise = AetheriaNoise(warp), Bowls = bowls.ToArray(), SlopesY = new float2(-0.25f, 0.1f), + CameraAt = (r, field) => + { + var angle = Uniform(r, 0.0f, 2.0f * MathF.PI); + var radius = ZoneRadius * MathF.Sqrt(r.NextSingle()); + var xz = new float2(radius * MathF.Cos(angle), radius * MathF.Sin(angle)); + var h = field.Height(field.Bowls.Select(b => b.Centre - xz).ToArray(), 0.0f, 0.0f); + return new float3(xz.x, max(0.0f, Safety + 10.0f - h), xz.y); + }, + }; + } + + private sealed class FogTile : Tile + { + public float2[] Centres = Array.Empty(); // bowl centres in the camera's (x, z) + } + + public override Tile DrawTile(System.Random random, int n) + { + var origin = CameraAt(random, this); + var yaw = Uniform(random, 0.0f, 2.0f * MathF.PI); + var tile = new FogTile + { + Origin = origin, N = n, Focal = Focal, Time = Uniform(random, 0.0f, 100.0f), + Slope = new float2(Uniform(random, -1.0f, 1.0f), Uniform(random, SlopesY.x, SlopesY.y)), + }; + tile.Look(new float3(MathF.Sin(yaw), 0.0f, MathF.Cos(yaw))); + tile.Centres = Bowls.Select(b => + { + var c = tile.Camera(new float3(b.Centre.x, origin.y, b.Centre.y)); + return new float2(c.x, c.z); + }).ToArray(); + return tile; + } + + private static float Pulse(float distance, float scale, float exponent) + { + var x = distance / scale; + var t = 1.0f - 4.0f * x * x; + return t <= 0.0f ? 0.0f : pow(t, exponent); + } + + // h at the camera-frame point (x, z), bowl centres in the same frame. + private float Height(float2[] centres, float x, float z) + { + var h = 0.0f; + for (var b = 0; b < Bowls.Length; b++) + { + var dx = x - centres[b].x; + var dz = z - centres[b].y; + h += Pulse(sqrt(dx * dx + dz * dz), Bowls[b].Scale, Bowls[b].Exponent) * Bowls[b].Depth; + } + + return h; + } + + // Each bowl's term over the box is [pulse(far), pulse(near)], pulse decreasing in distance; pow is + // not required to be correctly rounded, so both ends are widened by one ulp. + private float2 HeightBound(float2[] centres, float2 x, float2 z) + { + var h = iv_point(0.0f); + for (var b = 0; b < Bowls.Length; b++) + { + var (nx, fx) = AxisDistance(x, centres[b].x); + var (nz, fz) = AxisDistance(z, centres[b].y); + var near = Pulse(sqrt(nx * nx + nz * nz), Bowls[b].Scale, Bowls[b].Exponent); + var far = Pulse(sqrt(fx * fx + fz * fz), Bowls[b].Scale, Bowls[b].Exponent); + var pulse = new float2(far > 0.0f ? MathF.BitDecrement(far) : 0.0f, near > 0.0f ? MathF.BitIncrement(near) : 0.0f); + h = iv_add(h, iv_scale(pulse, Bowls[b].Depth)); + } + + return h; + } + + private static float2 Compose(float2 s, float2 fade, float2 n) + { + var displaced = iv_add(s, iv_scale(iv_mul(fade, n), Amplitude)); + return new float2(max(0.0f, (FloorOffset - displaced.y) / FloorBlend), max(0.0f, (FloorOffset - displaced.x) / FloorBlend)); + } + + public override float Density(Tile t, float2 slope, float z, bool lod, Counts counts) + { + var tile = (FogTile)t; + var x = slope.x * z; + var y = slope.y * z; + counts.Envelope++; + var s = tile.Origin.y + y + Height(tile.Centres, x, z); + var fade = 1.0f - smoothstep(0.75f * Safety, Safety, s); + if (Compose(iv_point(s), iv_point(fade), FreeNoise).y <= 0.0f) + return 0.0f; + var p = tile.World(new float3(x, y, z)); + var coverage = Noise.Coverage(p, tile.Time, counts); + if (Compose(iv_point(s), iv_point(fade), iv_add(iv_point(coverage), iv_scale(new float2(-1.0f, 1.0f), 0.5f))).y <= 0.0f) + return 0.0f; + var n = coverage + Noise.Detail(p, counts) * 0.5f; + return max(0.0f, (FloorOffset - (s + fade * n * Amplitude)) / FloorBlend); + } + + public override float2 Bound(Tile t, float z0, float z1, bool gated, Counts counts) + { + var tile = (FogTile)t; + var depth = new float2(z0, z1); + counts.Envelope++; + var height = HeightBound(tile.Centres, iv_mul(depth, tile.SlopesX), depth); + var s = iv_add(iv_add(iv_point(tile.Origin.y), iv_mul(depth, tile.SlopesY)), height); + var fade = iv_sub(iv_point(1.0f), iv_smoothstep(0.75f * Safety, Safety, s)); + var envelope = Compose(s, fade, FreeNoise); + if (gated && envelope.y <= 0.0f) + return envelope; + var ball = iv_frustum_ball(tile.Slope, z0, z1, tile.Footprint, Noise.Warp); + if (gated) + { + var provable = false; + foreach (var low in new[] { true, false }) + { + var best = iv_add(Noise.CoverageBound(BestBall(Noise.CoverageRadius(ball.w), low), tile.Time), iv_scale(BestBall(Noise.DetailRadius(ball.w), low), 0.5f)); + provable |= Compose(s, fade, best).y <= 0.0f; + } + + if (!provable) + return envelope; + } + + var centre = tile.World(new float3(ball.x, ball.y, ball.z)); + var noise = iv_add(Noise.CoverageBound(Noise.CoverageBall(centre, ball.w, counts), tile.Time), iv_scale(Noise.DetailBall(centre, ball.w, counts), 0.5f)); + return Compose(s, fade, noise); + } + } + + /// + /// (d) r1's uniform slab, kept as the worst case: the map's warped coverage field (F0 0.35, cutoff 0.25, + /// depth 12, period 4, extinction 2, flow in xy at 0.15, warp 0.05) with no envelope, so every cell is a + /// noise question. The dense sample evaluates the coverage first and the detail only where the coverage + /// passes the cutoff; ungated it would cost three snoise, as r1 counted. + /// + internal sealed class SlabField : Scenario { 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 const float SlabExtinction = 2.0f; + public const float SlabWarp = 0.05f; - public float Step0 { get; init; } = 0.25f; - public float Warp { get; init; } = 0.05f; + public required WarpedNoise Noise { get; init; } - public float StepDense => Step0 * 0.5f; - public float StepMax => Step0 * 8.0f; - public float FlowAmplitude => Warp * 2.0f / Period; + public SlabField() => Grid = Cells(0.125f, Depth); - public long Evaluations; - - private float2 Flow(float3 p) + public static SlabField Create(float warp) => new() { - 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; - } + Name = "(d) uniform slab", Noise = new WarpedNoise { F0 = F0, Warp = warp, Period = 4.0f, FlowFrequency = 0.15f, FlowInXY = true }, + }; - public float Density(float3 p, float time) + public override Tile DrawTile(System.Random random, int n) => new() { - 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); + Origin = new float3(Uniform(random, -50.0f, 50.0f), Uniform(random, -50.0f, 50.0f), 0.0f), + Slope = new float2(Uniform(random, -0.5f, 0.5f), Uniform(random, -0.5f, 0.5f)), + N = n, Focal = Focal, Time = Uniform(random, 0.0f, 100.0f), + }; + + public override float Density(Tile tile, float2 slope, float z, bool lod, Counts counts) + { + var p = tile.World(new float3(slope.x * z, slope.y * z, z)); + var coverage = Noise.Coverage(p, tile.Time, counts); + if (coverage <= Cutoff) + return 0.0f; + counts.Snoise++; 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); + return max(0.0f, coverage - Cutoff) * SlabExtinction * 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) + public override float2 Bound(Tile tile, float z0, float z1, bool gated, Counts counts) { - 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); + var maybe = new float2(0.0f, float.MaxValue); + var ball = iv_frustum_ball(tile.Slope, z0, z1, tile.Footprint, Noise.Warp); + if (gated && Noise.CoverageBound(BestBall(Noise.CoverageRadius(ball.w), low: true), tile.Time).y > Cutoff) + return maybe; + var coverage = Noise.CoverageBound(Noise.CoverageBall(tile.World(new float3(ball.x, ball.y, ball.z)), ball.w, counts), tile.Time); + return coverage.y <= Cutoff ? new float2(0.0f, 0.0f) : maybe; } } /// - /// 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. + /// (e) The void, the shipped scene: a cloud volume everywhere, carved by one negative pseudo-Gaussian + /// sphere. Authored knobs: hollow radius Rh (density exactly 0 inside) and ramp width S (the density + /// rises to full over [Rh, Rh + S]). The brush (1 - (d / Rc)^2)^e with CARVE = 1.5 reaches full at Rc, + /// so Rc = Rh + S, and Rh = Rc sqrt(1 - CARVE^(-1/e)) gives e = -ln CARVE / ln(1 - (Rh / Rc)^2): + /// density = K max(0, 1 - CARVE (1 - (d' / Rc)^2)^e), exactly 0 for d' <= Rh and K for d' >= Rc, + /// K = 1/30. The noise displaces the wall: fade = smoothstep(Rh - S, Rh, d), d' = d + A fade n(warp(p)), + /// A = 20, F0 = 0.01 (wavelengths 100 and 25), warp D = 10. The sun orbits the origin at radius 175 with + /// period 72 s and rests on the hollow's bottom (c_v = sun + (0, Rh - 12.5, 0)); the camera looks at + /// sun + (0, 60, 0) and, shipped, sits 120 behind the sun along the orbit tangent and 45 above it. The + /// footprint-truncated field (lod) weights octave i by w_i(z) = 1 - smoothstep(l_i / 2m, l_i / m, + /// z theta), m = 4, theta = 1 / f; the bound takes each weight over [w_i(z1), 1], which holds both the + /// truncated and the full field. The dense grid is the fixed-step reference's: cells of min(10, S / 4) + /// out to |camera - c_v| + Rc + 5 / K, past where every ray has gone opaque. /// - [Fact] - public void IntervalSkipHalvesEvaluations() + internal sealed class VoidField : Scenario { - 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 }) + public const float Carve = 1.5f; + public const float Amplitude = 20.0f; + public const float K = 1.0f / 30.0f; + public const float StepFootprints = 4.0f; // k + public const float OctaveFootprints = 4.0f; // m + private static readonly float2 FreeNoise = new(-1.5f, 1.5f); + + public static readonly float[] HollowRadii = { 100.0f, 200.0f, 400.0f, 1000.0f }; + public static readonly float[] Ramps = { 10.0f, 50.0f, 150.0f }; + public static readonly (string Name, float Fraction)[] CameraOffsets = { ("centred", 0.0f), ("mid", 0.5f), ("low", 0.94f) }; + + public readonly float HollowRadius; + public readonly float Ramp; + public readonly float CarveRadius; + public readonly float Exponent; + public readonly float RampStart; + public readonly float FixedStep; + public readonly float CameraFraction; // negative: the shipped camera + public readonly WarpedNoise Noise = new() { F0 = 0.01f, Warp = 10.0f }; + + private float CoarseWavelength => 1.0f / Noise.F0; + private float FineWavelength => 1.0f / (4.0f * Noise.F0); + private float Theta => 1.0f / Focal; + + [System.Diagnostics.CodeAnalysis.SetsRequiredMembers] + public VoidField(float hollowRadius, float ramp, float cameraFraction, string camera) { - 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}"); + HollowRadius = hollowRadius; + Ramp = ramp; + CarveRadius = hollowRadius + ramp; + Exponent = (float)(-Math.Log(Carve) / Math.Log(1.0 - (double)hollowRadius / CarveRadius * hollowRadius / CarveRadius)); + RampStart = hollowRadius - ramp; + FixedStep = min(10.0f, ramp / 4.0f); + CameraFraction = cameraFraction; + Name = $"(e) void Rh={hollowRadius:R} camera={camera} S={ramp:R}"; + var reach = cameraFraction < 0.0f ? length(new float2(120.0f, 45.0f - (hollowRadius - 12.5f))) : cameraFraction * hollowRadius; + Grid = Cells(FixedStep, reach + CarveRadius + 5.0f / K); } - 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}"); + public static VoidField Shipped() => new(198.0f, 50.0f, -1.0f, "shipped"); + + private sealed class VoidTile : Tile + { + public float3 Cavity; // c_v in the camera frame + } + + public override Tile DrawTile(System.Random random, int n) + { + var time = Uniform(random, 0.0f, 72.0f); + var angle = 2.0f * MathF.PI * time / 72.0f; + var sun = new float3(175.0f * MathF.Cos(angle), 0.0f, 175.0f * MathF.Sin(angle)); + var tangent = new float3(-MathF.Sin(angle), 0.0f, MathF.Cos(angle)); + var cavity = sun + new float3(0.0f, HollowRadius - 12.5f, 0.0f); + var target = sun + new float3(0.0f, 60.0f, 0.0f); + float3 origin; + float3 forward; + if (CameraFraction < 0.0f) + { + origin = sun - tangent * 120.0f + new float3(0.0f, 45.0f, 0.0f); + forward = target - origin; + } + else + { + // The shipped camera's direction from the centre, at a fraction of this hollow's radius; + // centred, it looks along the orbit, slightly down. + origin = cavity + normalize(tangent * -120.0f + new float3(0.0f, 45.0f - (198.0f - 12.5f), 0.0f)) * (CameraFraction * HollowRadius); + forward = CameraFraction == 0.0f ? tangent * -1.0f + new float3(0.0f, -0.2f, 0.0f) : target - origin; + } + + var tile = new VoidTile + { + Origin = origin, N = n, Focal = Focal, Time = time, + Slope = new float2((random.Next(1920 / n) * n + n * 0.5f - 960.0f) / Focal, (random.Next(1080 / n) * n + n * 0.5f - 540.0f) / Focal), + }; + tile.Look(forward); + tile.Cavity = tile.Camera(cavity); + return tile; + } + + private static float Weight(float wavelength, float footprint) => + 1.0f - smoothstep(wavelength / (2.0f * OctaveFootprints), wavelength / OctaveFootprints, footprint); + + private float2 WeightBound(float wavelength, float z0, float z1) => + new(iv_sub(iv_point(1.0f), iv_smoothstep(wavelength / (2.0f * OctaveFootprints), wavelength / OctaveFootprints, new float2(z0 * Theta, z1 * Theta))).x, 1.0f); + + /// The footprint step: max(k z theta, min(S / 4, half the finest live wavelength)). + public float LodStep(float z) + { + var footprint = z * Theta; + var fine = Weight(FineWavelength, footprint) > 0.0f ? FineWavelength : Weight(CoarseWavelength, footprint) > 0.0f ? CoarseWavelength : float.PositiveInfinity; + return max(StepFootprints * footprint, min(Ramp / 4.0f, fine / 2.0f)); + } + + private float Profile(float displaced) + { + if (displaced <= HollowRadius) + return 0.0f; + var x = displaced / CarveRadius; + var g = displaced < CarveRadius ? pow(max(0.0f, 1.0f - x * x), Exponent) : 0.0f; + return K * max(0.0f, 1.0f - Carve * g); + } + + private float2 Compose(float2 d, float2 fade, float2 n) + { + var displaced = iv_add(d, iv_scale(iv_mul(fade, n), Amplitude)); + return new float2(displaced.x >= CarveRadius ? K : 0.0f, displaced.y <= HollowRadius ? 0.0f : K); + } + + public override float Density(Tile t, float2 slope, float z, bool lod, Counts counts) + { + var tile = (VoidTile)t; + var x = slope.x * z; + var y = slope.y * z; + counts.Envelope++; + var dx = x - tile.Cavity.x; + var dy = y - tile.Cavity.y; + var dz = z - tile.Cavity.z; + var d = sqrt(dx * dx + dy * dy + dz * dz); + var fade = smoothstep(RampStart, HollowRadius, d); + var free = Compose(iv_point(d), iv_point(fade), FreeNoise); + if (free.x == free.y) + return free.x; + var p = tile.World(new float3(x, y, z)); + var wc = lod ? Weight(CoarseWavelength, z * Theta) : 1.0f; + var wd = lod ? Weight(FineWavelength, z * Theta) : 1.0f; + var coverage = wc > 0.0f ? Noise.Coverage(p, tile.Time, counts) * wc : 0.0f; + var partial = Compose(iv_point(d), iv_point(fade), iv_add(iv_point(coverage), iv_scale(new float2(-wd, wd), 0.5f))); + if (partial.x == partial.y) + return partial.x; + var detail = wd > 0.0f ? Noise.Detail(p, counts) * 0.5f * wd : 0.0f; + return Profile(d + fade * (coverage + detail) * Amplitude); + } + + public override float2 Bound(Tile t, float z0, float z1, bool gated, Counts counts) + { + var tile = (VoidTile)t; + var depth = new float2(z0, z1); + counts.Envelope++; + var (nx, fx) = AxisDistance(iv_mul(depth, tile.SlopesX), tile.Cavity.x); + var (ny, fy) = AxisDistance(iv_mul(depth, tile.SlopesY), tile.Cavity.y); + var (nz, fz) = AxisDistance(depth, tile.Cavity.z); + var d = new float2(sqrt(nx * nx + ny * ny + nz * nz), sqrt(fx * fx + fy * fy + fz * fz)); + var fade = iv_smoothstep(RampStart, HollowRadius, d); + var envelope = Compose(d, fade, FreeNoise); + if (gated && envelope.y <= 0.0f) + return envelope; + var ball = iv_frustum_ball(tile.Slope, z0, z1, tile.Footprint, Noise.Warp); + var wc = WeightBound(CoarseWavelength, z0, z1); + var wd = WeightBound(FineWavelength, z0, z1); + if (gated) + { + var provable = false; + foreach (var low in new[] { true, false }) + { + var best = iv_add(iv_mul(Noise.CoverageBound(BestBall(Noise.CoverageRadius(ball.w), low), tile.Time), wc), iv_mul(iv_scale(BestBall(Noise.DetailRadius(ball.w), low), 0.5f), wd)); + provable |= Compose(d, fade, best).y <= 0.0f; + } + + if (!provable) + return envelope; + } + + var centre = tile.World(new float3(ball.x, ball.y, ball.z)); + var noise = iv_add(iv_mul(Noise.CoverageBound(Noise.CoverageBall(centre, ball.w, counts), tile.Time), wc), iv_mul(iv_scale(Noise.DetailBall(centre, ball.w, counts), 0.5f), wd)); + return Compose(d, fade, noise); + } } } From 1111537b712a32646e55f04931008089f92847b4 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Sat, 3 Oct 2026 00:04:54 +0200 Subject: [PATCH 07/14] CultMath intervals r2: the tile march on cell ranges and the scenario measurement MarchTile runs the pre-pass once per tile on cell-index ranges. It starts from the whole grid, skips and doubles on proof, halves on maybe, and runs to the grid's end. Every ray of the tile integrates exactly the dense cells with the dense march's Integrate. MarchLod is the void's footprint-aware march inside the unmasked runs. IntervalSkipHalvesEvaluations prints (a)-(d) at N in {1, 4, 8, 16} and warp in {0, D}, the void grid (Rh x camera offset x S), and the shipped point with three marches. It asserts transmittance agreement and skips naming saving-2x-scenarios when (a), (c) or (e) misses its 2x contract. Co-Authored-By: Claude Opus 5.5 --- .../tests/CultMath.Tests/NoiseBoundTests.cs | 198 ++++++++++++++++++ 1 file changed, 198 insertions(+) diff --git a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs index 74edd778..b1c39d14 100644 --- a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs @@ -995,4 +995,202 @@ public sealed class NoiseBoundTests return Compose(d, fade, noise); } } + + // One cell's transmittance: the density at the cell's mid-depth, over the ray's length through it. + private static float Integrate(Scenario field, Tile tile, float2 slope, float z0, float z1, bool lod, Counts counts) => + exp(-(field.Density(tile, slope, (z0 + z1) * 0.5f, lod, counts) * field.Extinction) * ((z1 - z0) * sqrt(slope.x * slope.x + slope.y * slope.y + 1.0f))); + + /// + /// The tile pre-pass on cell-index ranges: probe the slice of the range; provably empty skips it and + /// doubles the next range; otherwise a range of more than one cell halves, and a single cell is dense. + /// It starts with the whole grid, runs to the grid's end (a pre-pass knows no ray's transmittance) and + /// yields the dense cells in order. Its probes are counted once per tile. + /// + private static List MarchTile(Scenario field, Tile tile, Counts probes) + { + var grid = field.Grid; + var cells = grid.Length - 1; + var dense = new List(); + var j = 0; + var span = cells; + while (j < cells) + { + var end = Math.Min(j + span, cells); + if (field.Bound(tile, grid[j], grid[end], gated: true, probes).y <= 0.0f) + { + j = end; + span = Math.Min(span * 2, cells); + } + else if (end - j > 1) + { + span = (end - j + 1) / 2; + } + else + { + dense.Add(j); + j++; + } + } + + return dense; + } + + // The footprint-aware march, the void's shipped one: steps of LodStep(z) inside the unmasked runs only, + // octave weights applied. + private static float MarchLod(VoidField field, Tile tile, float2 slope, List cells, Counts counts, ref long steps) + { + var grid = field.Grid; + var transmittance = 1.0f; + var k = 0; + while (k < cells.Count && transmittance >= 0.02f) + { + var z = grid[cells[k]]; + var end = cells[k] + 1; + for (k++; k < cells.Count && cells[k] == end; k++) + end++; + while (z < grid[end] && transmittance >= 0.02f) + { + var next = min(z + field.LodStep(z), grid[end]); + transmittance *= Integrate(field, tile, slope, z, next, lod: true, counts); + steps++; + z = next; + } + } + + return transmittance; + } + + internal sealed class RunStats + { + public readonly Counts Dense = new(); + public readonly Counts Masked = new(); + public readonly Counts Probes = new(); + public readonly Counts Lod = new(); + public long Rays; + public long Tiles; + public long DenseSteps; + public long MaskedSteps; + public long LodSteps; + public long Unfinished; + public float MaxDifference; + public float MaxLodDifference; + + public double PerRay(double value) => value / Rays; + + // Dense cost over the tile march's: each ray's own dense cells plus the tile's probes, shared. + public double Ratio => Dense.Cost / (Masked.Cost + Probes.Cost); + + // What no amortization can beat: the dense cost over the cost of the dense cells alone. + public double Ceiling => Dense.Cost / Masked.Cost; + + public double LodRatio => Dense.Cost / (Lod.Cost + Probes.Cost); + } + + /// + /// Draws `tiles` tiles, runs the pre-pass once per tile, and marches each of the N^2 rays (one per pixel, + /// jittered inside it) densely over every cell and masked over the tile's dense cells, both with the + /// early-out at 0.02 and the same Integrate; with lod, the footprint-aware march as well. + /// + private static RunStats Run(Scenario field, int n, int tiles, int seed, bool lod = false) + { + var random = new System.Random(seed); + var stats = new RunStats(); + var grid = field.Grid; + for (var t = 0; t < tiles; t++) + { + var tile = field.DrawTile(random, n); + var cells = MarchTile(field, tile, stats.Probes); + stats.Tiles++; + for (var j = 0; j < n; j++) + for (var i = 0; i < n; i++) + { + var slope = tile.PixelSlope(i, j, random.NextSingle(), random.NextSingle()); + var dense = 1.0f; + for (var c = 0; c < grid.Length - 1 && dense >= 0.02f; c++, stats.DenseSteps++) + dense *= Integrate(field, tile, slope, grid[c], grid[c + 1], lod: false, stats.Dense); + var masked = 1.0f; + for (var c = 0; c < cells.Count && masked >= 0.02f; c++, stats.MaskedSteps++) + masked *= Integrate(field, tile, slope, grid[cells[c]], grid[cells[c] + 1], lod: false, stats.Masked); + stats.Rays++; + stats.MaxDifference = MathF.Max(stats.MaxDifference, MathF.Abs(dense - masked)); + if (dense >= 0.02f) + stats.Unfinished++; + if (lod) + { + var marched = MarchLod((VoidField)field, tile, slope, cells, stats.Lod, ref stats.LodSteps); + stats.MaxLodDifference = MathF.Max(stats.MaxLodDifference, MathF.Abs(dense - marched)); + } + } + } + + return stats; + } + + private void Report(string label, RunStats s) => output.WriteLine( + $"IV-REPORT {label}: dense snoise/ray {s.PerRay(s.Dense.Snoise):F2} env/ray {s.PerRay(s.Dense.Envelope):F2} (samples/ray {s.PerRay(s.DenseSteps):F2}); " + + $"tile snoise/ray {s.PerRay(s.Masked.Snoise + (double)s.Probes.Snoise):F2} env/ray {s.PerRay(s.Masked.Envelope + (double)s.Probes.Envelope):F2} " + + $"(probes/tile env {(double)s.Probes.Envelope / s.Tiles:F1} snoise {(double)s.Probes.Snoise / s.Tiles:F2}; dense cells/ray {s.PerRay(s.MaskedSteps):F2}); " + + $"cost ratio {s.Ratio:F2}x, ceiling {s.Ceiling:F2}x; max |dT| {s.MaxDifference:R}; unfinished rays {s.Unfinished}"); + + private void ReportVoid(string label, RunStats s) => output.WriteLine( + $"IV-REPORT {label}: steps/px ref {s.PerRay(s.DenseSteps):F2}, masked {s.PerRay(s.MaskedSteps):F2}, LOD {s.PerRay(s.LodSteps):F2} ({(double)s.DenseSteps / s.LodSteps:F2}x fewer); " + + $"snoise/px ref {s.PerRay(s.Dense.Snoise):F2}, masked {s.PerRay(s.Masked.Snoise):F2}, LOD {s.PerRay(s.Lod.Snoise):F2}; " + + $"cost/px ref {s.PerRay(s.Dense.Cost):F2}, masked {s.PerRay(s.Masked.Cost + s.Probes.Cost):F2}, LOD {s.PerRay(s.Lod.Cost + s.Probes.Cost):F2} incl. probes (ref/LOD {s.LodRatio:F2}x, ref/masked {s.Ratio:F2}x); " + + $"probes/tile env {(double)s.Probes.Envelope / s.Tiles:F1} snoise {(double)s.Probes.Snoise / s.Tiles:F2}; max |dT| masked {s.MaxDifference:R}, LOD {s.MaxLodDifference:R}; unfinished rays {s.Unfinished}"); + + /// + /// The proving ground (docs/cultmath-interval-ground-cut.md, "Pass 3"). Scenarios (a) height fog, + /// (b) inside the fog, (c) Aetheria's wells and (d) r1's uniform slab, each over 64 tiles at + /// N in {1, 4, 8, 16} and warp in {0, D}: the tile march integrates exactly the dense march's cells, so + /// the transmittance agrees within 1e-3 (a difference would be a skipped cell that was not empty), and + /// its combined cost must be at most half the dense march's for (a) and (c) at N = 8. (e) the void: the + /// grid over hollow radius, camera offset and ramp width (32 tiles at N = 8) and the shipped point + /// (200 tiles), three marches each; at the shipped point the masked fixed-step march agrees with the + /// reference exactly, the footprint-aware march within 0.02, and it must take at most half the + /// reference's steps per pixel. A shortfall in (a), (c) or (e) skips naming saving-2x-scenarios; the + /// fields, the grids and the envelope cost are not tuned toward it. + /// + [Fact] + public void IntervalSkipHalvesEvaluations() + { + var ratios = new Dictionary(); + var scenarios = new (string Key, float Warp, Func Make)[] + { + ("a", FogField.AetheriaWarp, FogField.HeightFog), + ("b", FogField.AetheriaWarp, FogField.InsideFog), + ("c", FogField.AetheriaWarp, FogField.Wells), + ("d", SlabField.SlabWarp, SlabField.Create), + }; + foreach (var (key, warp, make) in scenarios) + foreach (var n in new[] { 1, 4, 8, 16 }) + foreach (var w in new[] { 0.0f, warp }) + { + var field = make(w); + var stats = Run(field, n, 64, 0x5A7E + n); + Report($"{field.Name} N={n} warp={w:R}", stats); + Assert.True(stats.MaxDifference <= 1.0e-3f, $"{field.Name} N={n} warp={w:R}: transmittance differs by {stats.MaxDifference:R}"); + if (key == "d") + output.WriteLine($"IV-REPORT (d) N={n} warp={w:R}: ungated dense snoise/ray {3.0 * stats.PerRay(stats.DenseSteps):F2} (three per dense sample, as r1 counted)"); + if (n == 8) + ratios[key] = Math.Min(ratios.GetValueOrDefault(key, double.MaxValue), stats.Ratio); + } + + foreach (var rh in VoidField.HollowRadii) + foreach (var offset in VoidField.CameraOffsets) + foreach (var ramp in VoidField.Ramps) + { + var field = new VoidField(rh, ramp, offset.Fraction, offset.Name); + ReportVoid($"{field.Name} N=8", Run(field, 8, 32, 0x701D, lod: true)); + } + + var shipped = VoidField.Shipped(); + var headline = Run(shipped, 8, 200, 0x5417, lod: true); + ReportVoid($"{shipped.Name} N=8 (headline, 200 tiles)", headline); + Assert.True(headline.MaxDifference == 0.0f, $"the masked fixed-step march differs from the reference by {headline.MaxDifference:R}"); + Assert.True(headline.MaxLodDifference <= 0.02f, $"the footprint-aware march differs from the reference by {headline.MaxLodDifference:R}"); + var steps = (double)headline.DenseSteps / headline.LodSteps; + + if (ratios["a"] < 2.0 || ratios["c"] < 2.0 || steps < 2.0) + Assert.Skip($"saving-2x-scenarios: at N = 8, (a) saves {ratios["a"]:F2}x and (c) {ratios["c"]:F2}x combined cost (the lower of warp 0 and D); (e) takes {steps:F2}x fewer steps per pixel at the shipped point; the contract is 2x each"); + } } From 43183276d377041e478fd368114aed5d6df80e68 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Sat, 3 Oct 2026 00:08:35 +0200 Subject: [PATCH 08/14] CultMath design.md and README: iv_frustum_ball, the envelope composition, the scenario saving design.md "Intervals" gains the frustum-ball paragraph and the envelope-then-noise paragraph: how a consumer composes the bound, the short-circuit, and that the warp lives on the noise only. r1's slab saving paragraph is replaced by r2's printed numbers per scenario. README names iv_frustum_ball and the envelope composition. Co-Authored-By: Claude Opus 5.5 --- packages/cultmath/README.md | 11 ++-- packages/cultmath/docs/design.md | 90 ++++++++++++++++++++++++++------ 2 files changed, 83 insertions(+), 18 deletions(-) diff --git a/packages/cultmath/README.md b/packages/cultmath/README.md index 6b493a31..aba9c1a5 100644 --- a/packages/cultmath/README.md +++ b/packages/cultmath/README.md @@ -38,9 +38,14 @@ shader semantics. - interval arithmetic for culling empty space with a proof: `iv_*` functions over `float2(lo, hi)` (`iv_add`, `iv_mul`, `iv_exp`, `iv_smoothstep`, ...) and `iv_snoise_ball`/`iv_fbm_ball`, which bound `snoise` and `fbm_grad` over a - ball. Every one returns an interval containing the function's value at every - point of its input, so a raymarch can skip a segment whose `hi` is below its - cutoff; `cultmath_iv_*` in HLSL. See `docs/design.md`, "Intervals". + ball, and `iv_frustum_ball`, the ball around a screen tile's rays over a depth + segment, so one probe serves the whole tile. Every one returns an interval + containing the function's value at every point of its input, so a raymarch + can skip a segment whose `hi` is below its cutoff; `cultmath_iv_*` in HLSL. + The bounds compose with a consumer's own analytic envelope (a height fog, a + carved sphere): its interval over the segment proves empty space in one probe, + and the noise ball is consulted only near the surface. See `docs/design.md`, + "Intervals". - `Voronoi.SampleTones`, a C# batch surface that calls the Rust `cultmath-core` native kernel when `cultmath_core` is available and falls back to the managed parity path otherwise. diff --git a/packages/cultmath/docs/design.md b/packages/cultmath/docs/design.md index 4b139234..aa348744 100644 --- a/packages/cultmath/docs/design.md +++ b/packages/cultmath/docs/design.md @@ -400,6 +400,44 @@ warp that moves a point by at most `D` is enclosed by the radius `r + D` (`NoiseBoundTests.WarpedPointsStayEnclosed`); this is what the gamecult.org ground's march rests on. +`iv_frustum_ball(m_c, z0, z1, footprintPerDepth, warp)` is the ball one probe +serves a whole screen tile with. In a pinhole camera frame looking down `+z`, +where the ray of slope `m` is the points `(m z, z)`, it encloses every point of +every ray whose slope lies within `footprintPerDepth` of `m_c` (`N / (sqrt(2) f)` +for an `N x N` tile at focal length `f` pixels: the full pixel footprint, so +sub-pixel jitter is covered), over depths `[z0, z1]`, each point moved by a warp +of length at most `warp`. The centre is `(m_c z_m, z_m)`, `z_m` the mid-depth; +the radius is `((z1 - z0) / 2) |(m_c, 1)| + z1 footprintPerDepth + warp`. A point +of slope `m_c + e` at depth `z` is the centre plus `(z - z_m)(m_c, 1)` plus +`z (e, 0)`, and the two vectors are never parallel, so the bound is strict. The +caller rotates and translates the centre into world space and scales centre and +radius by its noise frequency. `NoiseBoundTests.TileBallEnclosesEveryRaySegment` +checks 2,000 seeded tiles x 64 points, tile corners at both depth ends with +full-length flows included, with no tolerance, and pins the ball to the +derivation, so a ball grown past it fails too. A moving camera needs nothing +more: each frame's probes use that frame's camera, and no ball, mask or probe +result is carried from one frame to the next. + +The noise ball alone proves little: at `L` near 10 it bounds `snoise` to a width +of `2 L r`, which covers all of `[-1, 1]` once `L r` reaches 1. Empty space +becomes provable through an analytic envelope. The consumer's density is a +pre-distortion SDF (a height fog, a well map, a carved sphere) that noise +displaces only near its surface, and the consumer composes the bound from `iv_*` +ops. The envelope is taken over the slice's exact box (`iv_mul` of the depth +range by each slope range, then the SDF's own monotone pieces), together with its +fade. The noise is `iv_snoise_ball` over `iv_frustum_ball`, scaled by the noise +frequency. The warp moves only the noise argument, so it enlarges the ball and +the envelope's box takes none of it. Where the fade's upper end is 0 the noise +term is identically 0 and no `snoise` is evaluated. Where the envelope already +proves the slice empty with the noise at its full range, the ball is not +evaluated either. A probe is worth attempting only when the best result any +centre value could give would prove the slice, and a ball with `L r >= 2` is +`[-1, 1]` whatever `snoise(c)` is, so it costs nothing. CultMath owns none of the +envelope: the interval files have no envelope function, and each consumer composes +its own. `NoiseBoundTests.EnvelopeBoundEnclosesDensity` checks that composition, +the one the site and Aetheria use, over 2,000 seeded slices x 64 warped points +with no tolerance. + `L` has provenance, not a proof. `NoiseBoundTests.MeasureLipschitz` (slow, explicit) takes the largest `|snoise_grad|` over 1e6 seeded points in [-256, 256]^3 (7.1933103), refines the 1e4 largest by gradient ascent @@ -420,21 +458,43 @@ tracking or affine forms; those are later steps of sampled range: 5.08 for `iv_snoise_ball` at `r = 0.05`, 1.96 at 0.25 and 1.28 at 1.0; 6.64 and 2.82 for 4-octave `iv_fbm_ball` at 0.05 and 0.25. -The saving is a number too. `NoiseBoundTests.IntervalSkipHalvesEvaluations` -marches 1,000 seeded rays through the map's warped coverage field (`F0 = 0.35`, -cutoff 0.25, a slab of depth 12) once at the dense step and once with the -adaptive interval march, and counts `snoise` evaluations (three per dense -sample, one per probe; the flow is not counted, because the site reads it from -a texture). Both marches integrate the same cells, so their transmittance -agrees exactly. At the test's parameters (first step 0.25, dense step 0.125, -warp bound 0.05) the dense march spends 287,940 evaluations and the interval -march 309,550: 0.93x, not the 2x saving the proving-ground contract asks for. -Over the grid the test also prints, the best is 1.83x (first step 0.0625, no -warp), and at a first step of 1 it is 0.74x. With `L` near 10 a probe proves a -segment empty only where `snoise` is far below the cutoff, and a probe that -fails costs one evaluation on top of the three the dense sample costs. The test -skips with that open fork named instead of passing; whether the march, the -field or the bound changes is the campaign's decision. +The saving is measured where there is empty space to find. +`NoiseBoundTests.IntervalSkipHalvesEvaluations` draws screen tiles, each with +its own camera, and runs one pre-pass per tile over the dense grid's cell ranges: +it starts with the whole grid, skips and doubles a range it proves empty, halves +one it cannot, and runs to the grid's end. Each of the tile's `N^2` rays is then +marched twice with the same integrator and the same early-out at transmittance +0.02: once over every cell, once over the tile's unproven cells only. The two +integrate the same cells, so their transmittance agrees exactly in every +configuration. Cost is counted in `snoise` evaluations, with an envelope +evaluation or probe weighted 0.2. That weight is an estimate, and both counts are +printed. The flow is not counted. At `N = 8` with the warp on, in Aetheria's units +on its 256-cell quadratic grid: + +- (a) Height fog, camera above the safety band, rays level and up: the whole + grid is one envelope probe per tile, against 256 envelope evaluations per ray, + 16384x. +- (b) Inside the fog: 0.99x. Every cell is fog, and the pre-pass probes each one. +- (c) The zone bowl and four wells, the camera gazing across: 1.41x, and 3.32x + with no warp. The pre-pass cuts envelope evaluations from 212.93 to 37.12 per + ray. Almost all of the dense march's 78.67 `snoise` per ray fall in the noise + band, where the rays need them (ceiling 1.42x). The warp bound `D = 60` is what + leaves the noise ball unable to prove anything at this noise scale. +- (d) r1's uniform slab, which has no envelope: 1.23x, and 1.48x with no warp. +- (e) The void at the shipped point (`Rh = 198`, `S = 50`, the low camera, + 200 tiles), with three marches. The fixed-step reference takes 36.34 steps per + pixel. The same cells masked by the pre-pass take 20.13, with identical + transmittance. The footprint-aware march takes 16.45 inside the unmasked + cells, within 0.008339431 of the reference's transmittance: 2.21x fewer steps. + `snoise` per pixel goes from 32.98 to 26.56, and combined cost improves 1.34x. + Over the grid of hollow radius, camera offset and ramp width the steps ratio + runs from 1.57x (`Rh = 100`, low, `S = 150`) to 12.40x (`Rh = 1000`, centred, + `S = 10`), and the cost ratio from 1.12x to 3.47x. + +The contract is 2x on combined cost for (a) and (c) and 2x fewer steps per pixel +for (e). (a) and (e) meet it. (c) does not, so the test skips naming +`saving-2x-scenarios`. The fields, the grids and the 0.2 weight are not tuned +toward the contract. Consumers: the gamecult.org ground shader (through the GLSL lowering) and Aetheria's nebula raymarch at its CultMath pin bump. From 5f6bbfd4b51a84af0854b1845fd03d10fa7cdbc3 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Sat, 3 Oct 2026 11:37:30 +0200 Subject: [PATCH 09/14] CultMath intervals r3: iv_frustum_ball widened by its derived float32 rounding bound; extreme-family enclosure The radius gains (|c|_1 + radius) 2^-20 in C# and both HLSL copies; the doc's strictness claim leaves. TileBallEnclosesEveryRaySegment measures in double over the r2 domain and seven extreme families with outward corner flows, pins the radius two-sided, and counts misses per family. The noise ball enclosure tests draw radii to 1e3. Co-Authored-By: Claude Opus 5.5 --- packages/cultmath/README.md | 3 +- packages/cultmath/docs/design.md | 26 ++- .../cultmath/shaders/CultMath.Interval.hlsl | 8 +- .../cultmath/src/CultMath/math.Interval.cs | 19 +- .../tests/CultMath.Tests/NoiseBoundTests.cs | 165 +++++++++++++----- .../Shaders/CultMath.Interval.hlsl | 8 +- 6 files changed, 173 insertions(+), 56 deletions(-) diff --git a/packages/cultmath/README.md b/packages/cultmath/README.md index aba9c1a5..81d15859 100644 --- a/packages/cultmath/README.md +++ b/packages/cultmath/README.md @@ -39,7 +39,8 @@ shader semantics. over `float2(lo, hi)` (`iv_add`, `iv_mul`, `iv_exp`, `iv_smoothstep`, ...) and `iv_snoise_ball`/`iv_fbm_ball`, which bound `snoise` and `fbm_grad` over a ball, and `iv_frustum_ball`, the ball around a screen tile's rays over a depth - segment, so one probe serves the whole tile. Every one returns an interval + segment, so one probe serves the whole tile (its radius carries its own + float32 rounding bound). Every one returns an interval containing the function's value at every point of its input, so a raymarch can skip a segment whose `hi` is below its cutoff; `cultmath_iv_*` in HLSL. The bounds compose with a consumer's own analytic envelope (a height fog, a diff --git a/packages/cultmath/docs/design.md b/packages/cultmath/docs/design.md index aa348744..77688dc5 100644 --- a/packages/cultmath/docs/design.md +++ b/packages/cultmath/docs/design.md @@ -409,12 +409,26 @@ sub-pixel jitter is covered), over depths `[z0, z1]`, each point moved by a warp of length at most `warp`. The centre is `(m_c z_m, z_m)`, `z_m` the mid-depth; the radius is `((z1 - z0) / 2) |(m_c, 1)| + z1 footprintPerDepth + warp`. A point of slope `m_c + e` at depth `z` is the centre plus `(z - z_m)(m_c, 1)` plus -`z (e, 0)`, and the two vectors are never parallel, so the bound is strict. The -caller rotates and translates the centre into world space and scales centre and -radius by its noise frequency. `NoiseBoundTests.TileBallEnclosesEveryRaySegment` -checks 2,000 seeded tiles x 64 points, tile corners at both depth ends with -full-length flows included, with no tolerance, and pins the ball to the -derivation, so a ball grown past it fails too. A moving camera needs nothing +`z (e, 0)`. The returned ball is float32, and on a degenerate (`z0 == z1`), +thin or far segment the rounding of the centre alone exceeds the gap the +triangle inequality leaves, so the radius is widened by +`(|c|_1 + radius) 2^-20`. With `u = 2^-24`, each centre component carries at most +two roundings (`z_m`, then the product), so the float centre is within +`2u |c|_1` of the exact one; the radius carries at most seven (`z1 - z0`, the +axis length at `3u`, the product, `z1 footprintPerDepth`, the sum, `+ warp`), so +the exact radius is at most `(1 + 7u)` times the float one; the widening rounds +once more. `2^-20 = 16u` covers `8u (|c|_1 + radius)` twice over, and fma +contraction only removes roundings, so HLSL and GLSL compilers are covered by +the same bound. At the shipped scales (`|c|_1` near 600, `r` near 5) it adds +`6e-4`. The caller rotates and translates the centre into world space and scales +centre and radius by its noise frequency; the rounding of that transform is the +caller's. `NoiseBoundTests.TileBallEnclosesEveryRaySegment` checks 2,000 seeded +tiles x 64 points of the r2 domain and 2,000 tiles x 8 corner points of each of +seven extreme families (`z0 == z1`, `z1 - z0 = 1e-6 z1`, `|m_c|` to 1000, `z` to +`1e7`, `z0 = 0` with `z1` to `1e-3`, `f` down to 1, warp to `1e6`), corner flows +full-length and outward, distances in double, with no tolerance. It pins the ball +to the derivation both ways (`r <= ball.w <= r + 2^-19 (r + |c|_1)`), so a ball +that forgets the widening fails and one grown past it fails too. A moving camera needs nothing more: each frame's probes use that frame's camera, and no ball, mask or probe result is carried from one frame to the next. diff --git a/packages/cultmath/shaders/CultMath.Interval.hlsl b/packages/cultmath/shaders/CultMath.Interval.hlsl index 732ef3a0..b6b16665 100644 --- a/packages/cultmath/shaders/CultMath.Interval.hlsl +++ b/packages/cultmath/shaders/CultMath.Interval.hlsl @@ -106,13 +106,17 @@ float2 cultmath_iv_fbm_ball(float3 centre, float radius, int octaves, float lacu // A ball (xyz centre, w radius) enclosing every point of every ray of a screen tile over the depth segment // [z0, z1], warped by at most warp. Camera frame: pinhole at the origin looking down +z, ray (m z, z). -// footprintPerDepth = N / (sqrt(2) f) for an N x N tile at focal length f. Derivation in math.Interval.cs. +// footprintPerDepth = N / (sqrt(2) f) for an N x N tile at focal length f. The radius is widened by +// (|c|_1 + radius) 2^-20, which covers the ball's own float32 rounding. Derivation in math.Interval.cs. float4 cultmath_iv_frustum_ball(float2 centreSlope, float z0, float z1, float footprintPerDepth, float warp) { float zm = (z0 + z1) * 0.5; float axis = sqrt(centreSlope.x * centreSlope.x + centreSlope.y * centreSlope.y + 1.0); + float cx = centreSlope.x * zm; + float cy = centreSlope.y * zm; float radius = (z1 - z0) * 0.5 * axis + z1 * footprintPerDepth + warp; - return float4(centreSlope.x * zm, centreSlope.y * zm, zm, radius); + radius += (abs(cx) + abs(cy) + abs(zm) + radius) * 9.5367431640625e-7; + return float4(cx, cy, zm, radius); } #endif diff --git a/packages/cultmath/src/CultMath/math.Interval.cs b/packages/cultmath/src/CultMath/math.Interval.cs index 547c4acc..91e99f87 100644 --- a/packages/cultmath/src/CultMath/math.Interval.cs +++ b/packages/cultmath/src/CultMath/math.Interval.cs @@ -175,16 +175,25 @@ public static partial class math /// pixel footprint, so sub-pixel jitter is covered). The centre is (m_c z_m, z_m), z_m = (z0 + z1) / 2; /// the radius is ((z1 - z0) / 2) |(m_c, 1)| + z1 footprintPerDepth + warp. Derivation: a point of slope /// m_c + e (|e| <= footprintPerDepth) at depth z is c + (z - z_m)(m_c, 1) + z (e, 0); the first term is - /// at most the first radius term, the second at most the second, and the warp adds at most warp - /// (NoiseBoundTests.TileBallEnclosesEveryRaySegment). The two vectors are never parallel, so the - /// triangle inequality is strict. The caller rotates and translates the centre into world space (the - /// radius is unchanged) and scales both by its noise frequency. + /// at most the first radius term, the second at most the second, and the warp adds at most warp. + /// The ball also carries its own float32 rounding: the radius is widened by (|c|_1 + radius) 2^-20. + /// With u = 2^-24, each centre component carries at most two roundings (z_m, then the product), so the + /// float centre is within 2u |c|_1 of the exact one; the radius carries at most seven, so the exact + /// radius is at most (1 + 7u) times the float one; the widening rounds once more. 2^-20 = 16u covers + /// 8u (|c|_1 + radius) twice over, and fma contraction only removes roundings. So every exact point of + /// every ray, at any depth and slope, lies in the returned float ball + /// (NoiseBoundTests.TileBallEnclosesEveryRaySegment, degenerate and far segments included). The caller + /// rotates and translates the centre into world space (the radius is unchanged) and scales both by its + /// noise frequency; the rounding of that transform is the caller's. /// public static float4 iv_frustum_ball(float2 centreSlope, float z0, float z1, float footprintPerDepth, float warp) { var zm = (z0 + z1) * 0.5f; var axis = sqrt(centreSlope.x * centreSlope.x + centreSlope.y * centreSlope.y + 1.0f); + var cx = centreSlope.x * zm; + var cy = centreSlope.y * zm; var radius = (z1 - z0) * 0.5f * axis + z1 * footprintPerDepth + warp; - return new float4(centreSlope.x * zm, centreSlope.y * zm, zm, radius); + radius += (abs(cx) + abs(cy) + abs(zm) + radius) * 9.5367431640625e-7f; + return new float4(cx, cy, zm, radius); } } diff --git a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs index b1c39d14..d895828e 100644 --- a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs @@ -13,7 +13,11 @@ namespace CultMath.Tests; /// starting IV-REPORT carry the numbers the cut report quotes. /// /// Hand mutations of iv_frustum_ball that TileBallEnclosesEveryRaySegment must kill: z0 for z1 in the -/// lateral term (z0 * footprintPerDepth); |m_c| for |(m_c, 1)| in the depth term; the warp dropped. +/// lateral term (z0 * footprintPerDepth); |m_c| for |(m_c, 1)| in the depth term; the warp dropped; the +/// rounding widening dropped (which the degenerate family's enclosure check alone must catch). Hand +/// mutations of the test's own fields: w(z0) for w(z1) in VoidField.WeightBound must fail the truncated +/// check of EnvelopeBoundEnclosesDensity on the LOD-far configuration; CarveRadius for HollowRadius in +/// the lo test of VoidField.Compose must fail AnalyticBodyMatchesFineMarch. /// public sealed class NoiseBoundTests { @@ -55,7 +59,7 @@ public sealed class NoiseBoundTests for (var b = 0; b < 2000; b++) { var centre = RandomCentre(random); - var radius = LogUniform(random, 1.0e-3f, 2.0f); + var radius = LogUniform(random, 1.0e-3f, 1.0e3f); var bound = iv_snoise_ball(centre, radius); for (var i = 0; i < 64; i++) { @@ -72,7 +76,7 @@ public sealed class NoiseBoundTests for (var b = 0; b < 2000; b++) { var centre = RandomCentre(random); - var radius = LogUniform(random, 1.0e-3f, 2.0f); + var radius = LogUniform(random, 1.0e-3f, 1.0e3f); var octaves = random.Next(0, 7); var lacunarity = Uniform(random, 1.5f, 2.5f); var gain = Uniform(random, 0.3f, 0.7f); @@ -114,53 +118,134 @@ public sealed class NoiseBoundTests // ---- The tile ball and the composed bound ---- + // 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) + { + int[] sizes = { 1, 4, 8, 16 }; + float[] depths = { 12.0f, 1000.0f, 2400.0f }; + var n = sizes[random.Next(sizes.Length)]; + var focal = Uniform(random, 500.0f, 3000.0f); + var slope = new float2(Uniform(random, -0.77f, 0.77f), Uniform(random, -0.77f, 0.77f)); + var depth = depths[random.Next(depths.Length)]; + var z0 = random.Next(4) == 0 ? 0.0f : Uniform(random, 0.0f, depth * 0.9f); + var z1 = z0 + LogUniform(random, 1.0e-3f * depth, depth - z0); + var warp = random.Next(4) == 0 ? 0.0f : LogUniform(random, 1.0e-3f, 60.0f); + float Signed(float magnitude) => random.Next(2) == 0 ? -magnitude : magnitude; + switch (family) + { + case "degenerate": + z0 = z1 = LogUniform(random, 1.0e-3f * depth, depth); + break; + case "thin": + z0 = z1 - 1.0e-6f * z1; + break; + case "wide slope": + slope = new float2(Signed(LogUniform(random, 1.0e-3f, 707.0f)), Signed(LogUniform(random, 1.0e-3f, 707.0f))); + break; + case "far": + z1 = LogUniform(random, 1.0e4f, 1.0e7f); + z0 = z1 - z1 * LogUniform(random, 1.0e-6f, 1.0f); + break; + case "from the camera": + z0 = 0.0f; + z1 = LogUniform(random, 1.0e-3f, 1.0f); + break; + case "tiny focal": + focal = LogUniform(random, 1.0f, 500.0f); + break; + case "huge warp": + warp = LogUniform(random, 60.0f, 1.0e6f); + break; + } + + return (n, focal, slope, z0, z1, warp); + } + /// - /// iv_frustum_ball encloses every point of every ray of an N x N tile over a depth segment, sub-pixel - /// jitter included, each point then moved by a flow of length at most warp, with no tolerance: 2,000 - /// seeded tiles (N in {1, 4, 8, 16}, f in [500, 3000], |m_c| up to 1.1, segments over the slab's, the - /// fog's and the void's depth ranges, a quarter starting at the camera) x 64 points, the first eight at - /// a corner of the footprint at either end of the segment with a full-length flow. It also pins the - /// ball to the derivation: centre and radius equal the formula evaluated in double to one part in - /// 1e5, so a ball that grows (sound but wasteful) fails too. Each tile's camera is drawn on its own, as - /// a moving camera's frames are: no ball, mask or probe result is carried from one to the next. + /// iv_frustum_ball encloses every exact point of every ray of an N x N tile over a depth segment, + /// sub-pixel jitter included, each point then moved by a flow of length at most warp, with no + /// tolerance: the distance is taken in double from the float32 ball to the point (m z, z) formed in + /// double from the float32 slope and depth. 2,000 seeded tiles of the r2 domain (N in {1, 4, 8, 16}, + /// f in [500, 3000], |m_c| up to 1.1, segments over the slab's, the fog's and the void's depth ranges, + /// a quarter starting at the camera) x 64 points, and 2,000 tiles of each extreme family (z0 == z1; + /// z1 - z0 = 1e-6 z1; |m_c| up to 1000; z up to 1e7; z0 = 0, z1 down to 1e-3; f down to 1; warp up to + /// 1e6) x the 8 points at the footprint's corners at both depth ends. Every corner point's flow is + /// full-length and points outward from the ball's centre. It also pins the ball to the derivation in + /// double, both ways: the centre within 2^-22 |c|_1, and r <= ball.w <= r + 2^-19 (r + |c|_1), so a + /// ball that forgets its rounding widening fails and one that grows past it fails. Misses and pin + /// failures are counted per family and reported together. Each tile's camera is drawn on its own, as a + /// moving camera's frames are: no ball, mask or probe result is carried from one to the next. /// [Fact] public void TileBallEnclosesEveryRaySegment() { var random = new System.Random(0x711E); - int[] sizes = { 1, 4, 8, 16 }; - float[] depths = { 12.0f, 1000.0f, 2400.0f }; - for (var t = 0; t < 2000; t++) + string[] families = { "r2 domain", "degenerate", "thin", "wide slope", "far", "from the camera", "tiny focal", "huge warp" }; + var misses = new Dictionary(); + var pins = new Dictionary(); + var witness = new List(); + foreach (var family in families) { - var n = sizes[random.Next(sizes.Length)]; - var focal = Uniform(random, 500.0f, 3000.0f); - var slope = new float2(Uniform(random, -0.77f, 0.77f), Uniform(random, -0.77f, 0.77f)); - var depth = depths[random.Next(depths.Length)]; - var z0 = random.Next(4) == 0 ? 0.0f : Uniform(random, 0.0f, depth * 0.9f); - var z1 = z0 + LogUniform(random, 1.0e-3f * depth, depth - z0); - var warp = random.Next(4) == 0 ? 0.0f : LogUniform(random, 1.0e-3f, 60.0f); - var half = n / (2.0f * focal); - var footprint = n / (MathF.Sqrt(2.0f) * focal); - var ball = iv_frustum_ball(slope, z0, z1, footprint, warp); - var centre = new float3(ball.x, ball.y, ball.z); - - var zm = (z0 + (double)z1) / 2.0; - var radius = (z1 - (double)z0) / 2.0 * Math.Sqrt((double)slope.x * slope.x + (double)slope.y * slope.y + 1.0) + z1 * (double)footprint + warp; - Assert.True(Math.Abs(ball.w - radius) <= 1.0e-5 * radius, $"radius {ball.w:R}, derivation {radius:R}"); - Assert.True(Math.Abs(ball.x - slope.x * zm) + Math.Abs(ball.y - slope.y * zm) + Math.Abs(ball.z - zm) <= 1.0e-5 * zm, $"centre {centre}, derivation ({slope.x * zm:R}, {slope.y * zm:R}, {zm:R})"); - - for (var i = 0; i < 64; i++) + misses[family] = 0; + pins[family] = 0; + for (var t = 0; t < 2000; t++) { - 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 flow = UnitVector(random) * (i < 8 ? warp : warp * MathF.Sqrt(random.NextSingle())); - var p = new float3(m.x * z, m.y * z, z) + flow; - Assert.True(length(p - centre) <= ball.w, $"iv_frustum_ball({slope}, {z0:R}, {z1:R}, {footprint:R}, {warp:R}) = {ball} misses {p} (slope {m}, depth {z:R}), distance {length(p - centre):R}"); + var (n, focal, slope, z0, z1, warp) = DrawBallCase(random, family); + var half = n / (2.0f * focal); + var footprint = n / (MathF.Sqrt(2.0f) * focal); + var ball = iv_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 zm = (z0 + (double)z1) / 2.0; + var (ex, ey) = (slope.x * zm, slope.y * zm); + var c1 = Math.Abs(ex) + Math.Abs(ey) + Math.Abs(zm); + var radius = (z1 - (double)z0) / 2.0 * Math.Sqrt((double)slope.x * slope.x + (double)slope.y * slope.y + 1.0) + z1 * (double)footprint + warp; + var centreOff = Math.Abs(cx - ex) + Math.Abs(cy - ey) + Math.Abs(cz - zm); + if (!(centreOff <= Math.ScaleB(c1, -22) && radius <= w && w <= radius + Math.ScaleB(radius + c1, -19))) + { + pins[family]++; + if (witness.Count < 8) + witness.Add($"{family} pin: iv_frustum_ball({slope}, {z0:R}, {z1:R}, {footprint:R}, {warp:R}) = {ball}; derivation centre ({ex:R}, {ey:R}, {zm:R}), radius {radius: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 (px, py, pz) = (m.x * (double)z, m.y * (double)z, (double)z); + double fx, fy, fz; + if (i < 8) + { + // Outward: along the corner's direction from the centre, full length. + var (dx, dy, dz) = (px - cx, py - cy, pz - cz); + var away = Math.Sqrt(dx * dx + dy * dy + dz * dz); + (fx, fy, fz) = away > 0.0 ? (dx / away * warp, dy / away * warp, dz / away * warp) : (warp, 0.0, 0.0); + } + else + { + var flow = UnitVector(random) * (warp * MathF.Sqrt(random.NextSingle())); + (fx, fy, fz) = (flow.x, flow.y, flow.z); + } + + var (qx, qy, qz) = (px + fx - cx, py + fy - cy, pz + fz - cz); + var distance = Math.Sqrt(qx * qx + qy * qy + qz * qz); + if (!(distance <= w)) + { + misses[family]++; + if (witness.Count < 8) + witness.Add($"{family}: iv_frustum_ball({slope}, {z0:R}, {z1:R}, {footprint:R}, {warp:R}) = {ball} misses slope {m}, depth {z:R} by {(distance - w) / w:R} r"); + } + } } } + + Assert.True(misses.Values.Sum() + pins.Values.Sum() == 0, + "misses " + string.Join(", ", families.Select(f => $"{f} {misses[f]}")) + "; pin failures " + string.Join(", ", families.Select(f => $"{f} {pins[f]}")) + + Environment.NewLine + string.Join(Environment.NewLine, witness)); } /// diff --git a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl index 732ef3a0..b6b16665 100644 --- a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl +++ b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Interval.hlsl @@ -106,13 +106,17 @@ float2 cultmath_iv_fbm_ball(float3 centre, float radius, int octaves, float lacu // A ball (xyz centre, w radius) enclosing every point of every ray of a screen tile over the depth segment // [z0, z1], warped by at most warp. Camera frame: pinhole at the origin looking down +z, ray (m z, z). -// footprintPerDepth = N / (sqrt(2) f) for an N x N tile at focal length f. Derivation in math.Interval.cs. +// footprintPerDepth = N / (sqrt(2) f) for an N x N tile at focal length f. The radius is widened by +// (|c|_1 + radius) 2^-20, which covers the ball's own float32 rounding. Derivation in math.Interval.cs. float4 cultmath_iv_frustum_ball(float2 centreSlope, float z0, float z1, float footprintPerDepth, float warp) { float zm = (z0 + z1) * 0.5; float axis = sqrt(centreSlope.x * centreSlope.x + centreSlope.y * centreSlope.y + 1.0); + float cx = centreSlope.x * zm; + float cy = centreSlope.y * zm; float radius = (z1 - z0) * 0.5 * axis + z1 * footprintPerDepth + warp; - return float4(centreSlope.x * zm, centreSlope.y * zm, zm, radius); + radius += (abs(cx) + abs(cy) + abs(zm) + radius) * 9.5367431640625e-7; + return float4(cx, cy, zm, radius); } #endif From ebb12cbc938fed69945bc059e9e1e003739bfc42 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Sat, 3 Oct 2026 11:38:12 +0200 Subject: [PATCH 10/14] CultMath intervals r3: the fog grid's far plane is Aetheria's Main Camera's (2048), cited by line The Tint camera's 1000 was the wrong camera. Co-Authored-By: Claude Opus 5.5 --- packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs index d895828e..30ec1e95 100644 --- a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs @@ -676,8 +676,10 @@ public sealed class NoiseBoundTests public FogField() { - // Near 0.3 and far 1000: the ARPG scene's camera (Assets/Scenes/ARPG.unity:635-636). - Grid = Quadratic(0.3f, 1000.0f, 256); + // Near 0.3 and far 2048: the ARPG scene's Main Camera, which carries the cloud march + // (Aetheria d3075eb0 Assets/Scenes/ARPG.unity:31777 m_Name, :31830 far clip plane; + // Assets/Scripts/Zone Display/VolumeCloudRenderer.cs:36 RequireComponent(Camera)). + Grid = Quadratic(0.3f, 2048.0f, 256); Extinction = 0.5f; } From 3ac7250aa92569caefee45869c0c21174578936d Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Sat, 3 Oct 2026 11:42:57 +0200 Subject: [PATCH 11/14] CultMath intervals r3: the analytic body, LOD-far, one shipped void, (c) against its ceiling MarchTile reports a body start where the gated bound over the rest of the grid is a positive point and appends the rest unprobed; MarchLod takes one closed-form Beer-Lambert step there. AnalyticBodyMatchesFineMarch pins it against the ungated field and the fixed-step product, and pins that no (a)-(d) tile reports a body. VoidField takes a focal length; LOD-far (Rh 4000, S 50, f 467) joins the enclosure test, which now asserts it sampled weights in (0, 1) and at 0. The void grid ends at Rc + 6/K. The doc derives the one shipped void (198, 50) and records the [w(z1), 1] deviation. (c) is held to 0.90 of its probe-free ceiling. Co-Authored-By: Claude Opus 5.5 --- .../tests/CultMath.Tests/NoiseBoundTests.cs | 276 +++++++++++++++--- 1 file changed, 231 insertions(+), 45 deletions(-) diff --git a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs index 30ec1e95..301979f5 100644 --- a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs @@ -250,11 +250,12 @@ public sealed class NoiseBoundTests /// /// intervals-enclose for the composition the site and Aetheria use: over 2,000 seeded tile slices - /// (Aetheria's wells, inside the fog, the void at every grid point; N in {1, 4, 8, 16}; a camera drawn - /// per tile) the bound composed from the envelope over the slice's box, its fade, and the noise ball - /// enlarged by the warp encloses the density at 64 points of the slice, each sampled through the field's - /// own warp, with no tolerance. The void's footprint-truncated density (octave weights) is checked - /// against the same bound. + /// (Aetheria's wells, inside the fog, the void at every grid point, the shipped void and LOD-far; + /// N in {1, 4, 8, 16}; a camera drawn per tile) the bound composed from the envelope over the slice's + /// box, its fade, and the noise ball enlarged by the warp encloses the density at 64 points of the + /// slice, each sampled through the field's own warp, with no tolerance. The void's footprint-truncated density (octave weights) is checked + /// against the same bound; a quarter of the slices are LOD-far's, where the weights are strictly + /// between 0 and 1 and at 0, and the test asserts it sampled both. /// [Fact] public void EnvelopeBoundEnclosesDensity() @@ -266,10 +267,13 @@ public sealed class NoiseBoundTests foreach (var offset in VoidField.CameraOffsets) foreach (var ramp in VoidField.Ramps) fields.Add(new VoidField(rh, ramp, offset.Fraction, offset.Name)); + fields.Add(VoidField.Shipped()); + fields.Add(VoidField.LodFar()); var counts = new Counts(); + var (partial, zero) = (0, 0); for (var t = 0; t < 2000; t++) { - var field = t % 2 == 0 ? fields[(t / 2) % 2] : fields[2 + random.Next(fields.Count - 2)]; + var field = t % 2 == 0 ? fields[(t / 2) % 2] : t % 4 == 1 ? fields[^1] : fields[2 + random.Next(fields.Count - 2)]; var tile = field.DrawTile(random, sizes[random.Next(sizes.Length)]); var near = field.Grid[0]; var far = field.Grid[^1]; @@ -282,13 +286,20 @@ public sealed class NoiseBoundTests var z = i < 2 ? (i == 0 ? z0 : z1) : Uniform(random, z0, z1); var density = field.Density(tile, m, z, lod: false, counts); Assert.True(Inside(density, bound), $"{field.Name}: bound {bound} over [{z0:R}, {z1:R}] misses density {density:R} at slope {m}, depth {z:R}"); - if (field is VoidField) + if (field is VoidField voidField) { + var (wc, wd) = voidField.Weights(z); + partial += (wc > 0.0f && wc < 1.0f) || (wd > 0.0f && wd < 1.0f) ? 1 : 0; + zero += wd == 0.0f ? 1 : 0; var truncated = field.Density(tile, m, z, lod: true, counts); Assert.True(Inside(truncated, bound), $"{field.Name}: bound {bound} misses the truncated density {truncated:R} at slope {m}, depth {z:R}"); } } } + + // The truncated check is not vacuous: some samples carry an octave weight strictly inside (0, 1) + // and some a fine weight of 0 (LOD-far reaches both). + Assert.True(partial > 0 && zero > 0, $"octave weights in (0, 1) at {partial} samples and fine weight 0 at {zero}"); } // ---- The Lipschitz constant ---- @@ -899,19 +910,32 @@ public sealed class NoiseBoundTests } /// - /// (e) The void, the shipped scene: a cloud volume everywhere, carved by one negative pseudo-Gaussian - /// sphere. Authored knobs: hollow radius Rh (density exactly 0 inside) and ramp width S (the density - /// rises to full over [Rh, Rh + S]). The brush (1 - (d / Rc)^2)^e with CARVE = 1.5 reaches full at Rc, - /// so Rc = Rh + S, and Rh = Rc sqrt(1 - CARVE^(-1/e)) gives e = -ln CARVE / ln(1 - (Rh / Rc)^2): - /// density = K max(0, 1 - CARVE (1 - (d' / Rc)^2)^e), exactly 0 for d' <= Rh and K for d' >= Rc, - /// K = 1/30. The noise displaces the wall: fade = smoothstep(Rh - S, Rh, d), d' = d + A fade n(warp(p)), - /// A = 20, F0 = 0.01 (wavelengths 100 and 25), warp D = 10. The sun orbits the origin at radius 175 with - /// period 72 s and rests on the hollow's bottom (c_v = sun + (0, Rh - 12.5, 0)); the camera looks at - /// sun + (0, 60, 0) and, shipped, sits 120 behind the sun along the orbit tangent and 45 above it. The - /// footprint-truncated field (lod) weights octave i by w_i(z) = 1 - smoothstep(l_i / 2m, l_i / m, - /// z theta), m = 4, theta = 1 / f; the bound takes each weight over [w_i(z1), 1], which holds both the - /// truncated and the full field. The dense grid is the fixed-step reference's: cells of min(10, S / 4) - /// out to |camera - c_v| + Rc + 5 / K, past where every ray has gone opaque. + /// (e) The void: a cloud volume everywhere, carved by one negative pseudo-Gaussian sphere. Authored + /// knobs: hollow radius Rh (density exactly 0 inside) and ramp width S (the density rises to full over + /// [Rh, Rh + S]); everything else is derived from them. The brush (1 - (d / Rc)^2)^e with CARVE = 1.5 + /// reaches full at Rc, so Rc = Rh + S, and Rh = Rc sqrt(1 - CARVE^(-1/e)) gives + /// e = -ln CARVE / ln(1 - (Rh / Rc)^2): density = K max(0, 1 - CARVE (1 - (d' / Rc)^2)^e), exactly 0 for + /// d' <= Rh and K for d' >= Rc, K = 1/30. The noise displaces the wall: fade = + /// smoothstep(Rh - S, Rh, d), d' = d + A fade n(warp(p)), A = 20, F0 = 0.01 (wavelengths 100 and 25), + /// warp D = 10. The sun orbits the origin at radius 175 with period 72 s and rests on the hollow's + /// bottom (c_v = sun + (0, Rh - 12.5, 0)); the camera looks at sun + (0, 60, 0) and, shipped, sits 120 + /// behind the sun along the orbit tangent and 45 above it. + /// + /// The one shipped void (ruling operator-shipped-void) is Shipped(): (Rh, S) = (198, 50), so Rc = 248 + /// and e = -ln 1.5 / ln(1 - (198 / 248)^2), which the constructor derives; no other void is shipped. + /// + /// The footprint-truncated field (lod) weights octave i by w_i(z) = 1 - smoothstep(l_i / 2m, l_i / m, + /// z theta), m = 4, theta = 1 / f. Deviation from the spec's [w_i(z1), w_i(z0)]: the bound takes each + /// weight over [w_i(z1), 1], because one Bound serves the full field (w = 1) and the truncated one + /// (w in [w_i(z1), w_i(z0)], w decreasing in z); [w_i(z1), 1] contains both and iv_mul is monotone in + /// its interval argument, so the composed interval encloses both fields. The slack is only in far + /// tiles, whose step the footprint already sets. At f = 935 every weight is 1 inside every grid; + /// LodFar() (f = 467) reaches weights in (0, 1) and at 0. + /// + /// The dense grid is the fixed-step reference's: cells of min(10, S / 4) out to + /// |camera - c_v| + Rc + 6 / K. The density is provably K only from d >= Rc + A 1.5 (the free noise + /// interval), so 6 / K leaves at least 4.9 of optical depth past it and every ray goes opaque before + /// the grid ends. /// internal sealed class VoidField : Scenario { @@ -940,8 +964,9 @@ public sealed class NoiseBoundTests private float Theta => 1.0f / Focal; [System.Diagnostics.CodeAnalysis.SetsRequiredMembers] - public VoidField(float hollowRadius, float ramp, float cameraFraction, string camera) + public VoidField(float hollowRadius, float ramp, float cameraFraction, string camera, float focal = 935.0f) { + Focal = focal; HollowRadius = hollowRadius; Ramp = ramp; CarveRadius = hollowRadius + ramp; @@ -949,13 +974,23 @@ public sealed class NoiseBoundTests RampStart = hollowRadius - ramp; FixedStep = min(10.0f, ramp / 4.0f); CameraFraction = cameraFraction; - Name = $"(e) void Rh={hollowRadius:R} camera={camera} S={ramp:R}"; + Name = $"(e) void Rh={hollowRadius:R} camera={camera} S={ramp:R}" + (focal == 935.0f ? "" : $" f={focal:R}"); var reach = cameraFraction < 0.0f ? length(new float2(120.0f, 45.0f - (hollowRadius - 12.5f))) : cameraFraction * hollowRadius; - Grid = Cells(FixedStep, reach + CarveRadius + 5.0f / K); + Grid = Cells(FixedStep, reach + CarveRadius + 6.0f / K); } public static VoidField Shipped() => new(198.0f, 50.0f, -1.0f, "shipped"); + /// + /// The configuration that reaches the octave weights: Rh = 4000, S = 50, the low camera, f = 467 + /// (1080 rows at a 98 degree vertical field of view). The fine weight is below 1 past z = 1460 and 0 + /// past 2920, the coarse below 1 past 5840, and the grid reaches about 7990. Printed, not asserted. + /// + public static VoidField LodFar() => new(4000.0f, 50.0f, 0.94f, "low", 467.0f); + + /// The coarse and fine octave weights at depth z (1 for the full field). + public (float Coarse, float Fine) Weights(float z) => (Weight(CoarseWavelength, z * Theta), Weight(FineWavelength, z * Theta)); + private sealed class VoidTile : Tile { public float3 Cavity; // c_v in the camera frame @@ -1048,6 +1083,26 @@ public sealed class NoiseBoundTests return Profile(d + fade * (coverage + detail) * Amplitude); } + /// + /// The field with no envelope gate: both noise terms always evaluated, the profile applied to the + /// displaced distance. Density's shortcuts go through Compose; this does not, so it is the value a + /// wrong Compose cannot also corrupt. Not counted. + /// + public float Exact(Tile t, float2 slope, float z, bool lod) + { + var tile = (VoidTile)t; + var counts = new Counts(); + var (x, y) = (slope.x * z, slope.y * z); + var (dx, dy, dz) = (x - tile.Cavity.x, y - tile.Cavity.y, z - tile.Cavity.z); + var d = sqrt(dx * dx + dy * dy + dz * dz); + var fade = smoothstep(RampStart, HollowRadius, d); + var p = tile.World(new float3(x, y, z)); + var wc = lod ? Weight(CoarseWavelength, z * Theta) : 1.0f; + var wd = lod ? Weight(FineWavelength, z * Theta) : 1.0f; + var n = Noise.Coverage(p, tile.Time, counts) * wc + Noise.Detail(p, counts) * 0.5f * wd; + return Profile(d + fade * n * Amplitude); + } + public override float2 Bound(Tile t, float z0, float z1, bool gated, Counts counts) { var tile = (VoidTile)t; @@ -1092,8 +1147,15 @@ public sealed class NoiseBoundTests /// doubles the next range; otherwise a range of more than one cell halves, and a single cell is dense. /// It starts with the whole grid, runs to the grid's end (a pre-pass knows no ray's transmittance) and /// yields the dense cells in order. Its probes are counted once per tile. + /// + /// The body (ruling operator-analytic-wall-body): when a single dense cell's own bound is a positive + /// point (the cell is provably constant), one more probe over [grid[j], grid end]; if that is a + /// positive point too, every point of every ray of the tile over the rest of the grid has exactly that + /// density, so BodyStart = j, BodyDensity = its value, and every remaining cell is appended dense + /// without probing (adding cells is always sound, so the mask keeps its meaning). BodyStart is -1 when + /// no body is proven. /// - private static List MarchTile(Scenario field, Tile tile, Counts probes) + private static (List Dense, int BodyStart, float BodyDensity) MarchTile(Scenario field, Tile tile, Counts probes) { var grid = field.Grid; var cells = grid.Length - 1; @@ -1103,7 +1165,8 @@ public sealed class NoiseBoundTests while (j < cells) { var end = Math.Min(j + span, cells); - if (field.Bound(tile, grid[j], grid[end], gated: true, probes).y <= 0.0f) + var bound = field.Bound(tile, grid[j], grid[end], gated: true, probes); + if (bound.y <= 0.0f) { j = end; span = Math.Min(span * 2, cells); @@ -1114,19 +1177,35 @@ public sealed class NoiseBoundTests } else { + if (bound.x == bound.y && bound.x > 0.0f) + { + var rest = field.Bound(tile, grid[j], grid[cells], gated: true, probes); + if (rest.x == rest.y && rest.x > 0.0f) + { + dense.AddRange(Enumerable.Range(j, cells - j)); + return (dense, j, rest.x); + } + } + dense.Add(j); j++; } } - return dense; + return (dense, -1, 0.0f); } + // The body's one closed-form Beer-Lambert step from z to the grid's end, at the density the field's + // pointwise envelope gives at z (one envelope evaluation; no probe in the main pass). + private static float BodyFactor(Scenario field, Tile tile, float2 slope, float z, bool lod, Counts counts) => + exp(-(field.Density(tile, slope, z, lod, counts) * field.Extinction) * ((field.Grid[^1] - z) * sqrt(slope.x * slope.x + slope.y * slope.y + 1.0f))); + // The footprint-aware march, the void's shipped one: steps of LodStep(z) inside the unmasked runs only, - // octave weights applied. - private static float MarchLod(VoidField field, Tile tile, float2 slope, List cells, Counts counts, ref long steps) + // octave weights applied; at the body's start, one analytic step to the grid's end, and stop. + private static float MarchLod(VoidField field, Tile tile, float2 slope, List cells, int bodyStart, Counts counts, ref long steps, ref long bodySteps) { var grid = field.Grid; + var body = bodyStart >= 0 ? grid[bodyStart] : float.PositiveInfinity; var transmittance = 1.0f; var k = 0; while (k < cells.Count && transmittance >= 0.02f) @@ -1137,7 +1216,15 @@ public sealed class NoiseBoundTests end++; while (z < grid[end] && transmittance >= 0.02f) { - var next = min(z + field.LodStep(z), grid[end]); + if (z >= body) + { + transmittance *= BodyFactor(field, tile, slope, z, lod: true, counts); + steps++; + bodySteps++; + return transmittance; + } + + var next = min(min(z + field.LodStep(z), grid[end]), z < body ? body : grid[end]); transmittance *= Integrate(field, tile, slope, z, next, lod: true, counts); steps++; z = next; @@ -1158,6 +1245,8 @@ public sealed class NoiseBoundTests public long DenseSteps; public long MaskedSteps; public long LodSteps; + public long BodySteps; + public long BodyTiles; public long Unfinished; public float MaxDifference; public float MaxLodDifference; @@ -1186,8 +1275,10 @@ public sealed class NoiseBoundTests for (var t = 0; t < tiles; t++) { var tile = field.DrawTile(random, n); - var cells = MarchTile(field, tile, stats.Probes); + var (cells, bodyStart, _) = MarchTile(field, tile, stats.Probes); + Assert.True(field is VoidField || bodyStart < 0, $"{field.Name}: a body reported at cell {bodyStart}; only the void's density can be provably constant"); stats.Tiles++; + stats.BodyTiles += bodyStart >= 0 ? 1 : 0; for (var j = 0; j < n; j++) for (var i = 0; i < n; i++) { @@ -1204,7 +1295,7 @@ public sealed class NoiseBoundTests stats.Unfinished++; if (lod) { - var marched = MarchLod((VoidField)field, tile, slope, cells, stats.Lod, ref stats.LodSteps); + var marched = MarchLod((VoidField)field, tile, slope, cells, bodyStart, stats.Lod, ref stats.LodSteps, ref stats.BodySteps); stats.MaxLodDifference = MathF.Max(stats.MaxLodDifference, MathF.Abs(dense - marched)); } } @@ -1223,19 +1314,103 @@ public sealed class NoiseBoundTests $"IV-REPORT {label}: steps/px ref {s.PerRay(s.DenseSteps):F2}, masked {s.PerRay(s.MaskedSteps):F2}, LOD {s.PerRay(s.LodSteps):F2} ({(double)s.DenseSteps / s.LodSteps:F2}x fewer); " + $"snoise/px ref {s.PerRay(s.Dense.Snoise):F2}, masked {s.PerRay(s.Masked.Snoise):F2}, LOD {s.PerRay(s.Lod.Snoise):F2}; " + $"cost/px ref {s.PerRay(s.Dense.Cost):F2}, masked {s.PerRay(s.Masked.Cost + s.Probes.Cost):F2}, LOD {s.PerRay(s.Lod.Cost + s.Probes.Cost):F2} incl. probes (ref/LOD {s.LodRatio:F2}x, ref/masked {s.Ratio:F2}x); " - + $"probes/tile env {(double)s.Probes.Envelope / s.Tiles:F1} snoise {(double)s.Probes.Snoise / s.Tiles:F2}; max |dT| masked {s.MaxDifference:R}, LOD {s.MaxLodDifference:R}; unfinished rays {s.Unfinished}"); + + $"probes/tile env {(double)s.Probes.Envelope / s.Tiles:F1} snoise {(double)s.Probes.Snoise / s.Tiles:F2}; body tiles {s.BodyTiles}/{s.Tiles}, analytic steps/px {s.PerRay(s.BodySteps):F2}; " + + $"max |dT| masked {s.MaxDifference:R}, LOD {s.MaxLodDifference:R}; unfinished rays {s.Unfinished}"); /// - /// The proving ground (docs/cultmath-interval-ground-cut.md, "Pass 3"). Scenarios (a) height fog, - /// (b) inside the fog, (c) Aetheria's wells and (d) r1's uniform slab, each over 64 tiles at - /// N in {1, 4, 8, 16} and warp in {0, D}: the tile march integrates exactly the dense march's cells, so - /// the transmittance agrees within 1e-3 (a difference would be a skipped cell that was not empty), and - /// its combined cost must be at most half the dense march's for (a) and (c) at N = 8. (e) the void: the - /// grid over hollow radius, camera offset and ramp width (32 tiles at N = 8) and the shipped point - /// (200 tiles), three marches each; at the shipped point the masked fixed-step march agrees with the - /// reference exactly, the footprint-aware march within 0.02, and it must take at most half the - /// reference's steps per pixel. A shortfall in (a), (c) or (e) skips naming saving-2x-scenarios; the - /// fields, the grids and the envelope cost are not tuned toward it. + /// The analytic body's soundness pin (ruling operator-analytic-wall-body). Over 2,000 seeded tiles of + /// the (e) grid, the shipped void and LOD-far at N in {1, 4, 8, 16}, wherever the pre-pass reports a + /// body start b: (1) the gated bound over [grid[b], grid end] is a positive point equal to the reported + /// density, every cell from b to the end is in the mask, and at 64 random (pixel, sub-pixel, cell >= b, + /// depth in cell) points the full and the truncated Density and the ungated Exact field all equal it + /// exactly, no tolerance; (2) the one closed-form factor exp(-density ext L) over the body equals the + /// fixed-step product of the body cells' Integrate factors (no early-out, same cells) within 1e-4 + /// relative. (3) No tile of (a)-(d), at either warp, reports a body: their bounds are never a positive + /// point. + /// + [Fact] + public void AnalyticBodyMatchesFineMarch() + { + var random = new System.Random(0xB0D1); + int[] sizes = { 1, 4, 8, 16 }; + var voids = new List(); + foreach (var rh in VoidField.HollowRadii) + foreach (var offset in VoidField.CameraOffsets) + foreach (var ramp in VoidField.Ramps) + voids.Add(new VoidField(rh, ramp, offset.Fraction, offset.Name)); + voids.Add(VoidField.Shipped()); + voids.Add(VoidField.LodFar()); + var counts = new Counts(); + var bodies = new Dictionary(); + for (var t = 0; t < 2000; t++) + { + var field = t % 4 == 0 ? voids[^2] : t % 4 == 1 ? voids[^1] : voids[random.Next(voids.Count)]; + var tile = field.DrawTile(random, sizes[random.Next(sizes.Length)]); + var (cells, b, density) = MarchTile(field, tile, counts); + if (b < 0) + continue; + bodies[field.Name] = bodies.GetValueOrDefault(field.Name) + 1; + var grid = field.Grid; + var last = grid.Length - 1; + var rest = field.Bound(tile, grid[b], grid[last], gated: true, counts); + Assert.True(rest.x == rest.y && rest.x > 0.0f && rest.x == density, $"{field.Name}: body at cell {b} with density {density:R}, but the bound over the rest is {rest}"); + Assert.Equal(Enumerable.Range(b, last - b), cells.Skip(cells.Count - (last - b))); + for (var i = 0; i < 64; i++) + { + var slope = tile.PixelSlope(random.Next(tile.N), random.Next(tile.N), random.NextSingle(), random.NextSingle()); + var c = random.Next(b, last); + var z = Uniform(random, grid[c], grid[c + 1]); + foreach (var lod in new[] { false, true }) + { + var pointwise = field.Density(tile, slope, z, lod, counts); + var exact = field.Exact(tile, slope, z, lod); + Assert.True(pointwise == density && exact == density, $"{field.Name}: body from cell {b} claims {density:R}; at slope {slope}, depth {z:R} (cell {c}, lod {lod}) Density is {pointwise:R} and the ungated field {exact:R}"); + } + } + + var ray = tile.PixelSlope(random.Next(tile.N), random.Next(tile.N), random.NextSingle(), random.NextSingle()); + var product = 1.0; + for (var c = b; c < last; c++) + product *= Integrate(field, tile, ray, grid[c], grid[c + 1], lod: true, counts); + var tau = density * field.Extinction * ((grid[last] - grid[b]) * sqrt(ray.x * ray.x + ray.y * ray.y + 1.0f)); + var analytic = Math.Exp(-(double)tau); + Assert.True(Math.Abs(analytic - product) <= 1.0e-4 * product, $"{field.Name}: the analytic body factor {analytic:R} differs from the fixed-step product {product:R} over cells {b}..{last}"); + } + + output.WriteLine("IV-REPORT analytic body: tiles with a body " + string.Join(", ", bodies.OrderBy(kv => kv.Key).Select(kv => $"{kv.Key} {kv.Value}"))); + Assert.True(bodies.ContainsKey(VoidField.Shipped().Name) && bodies.ContainsKey(VoidField.LodFar().Name), "no body proven at the shipped void or at LOD-far"); + + var scenarios = new (float Warp, Func Make)[] + { + (FogField.AetheriaWarp, FogField.HeightFog), (FogField.AetheriaWarp, FogField.InsideFog), + (FogField.AetheriaWarp, FogField.Wells), (SlabField.SlabWarp, SlabField.Create), + }; + foreach (var (warp, make) in scenarios) + foreach (var w in new[] { 0.0f, warp }) + { + var field = make(w); + for (var t = 0; t < 250; t++) + { + var (_, b, _) = MarchTile(field, field.DrawTile(random, sizes[random.Next(sizes.Length)]), counts); + Assert.True(b < 0, $"{field.Name} warp {w:R}: a body reported at cell {b}"); + } + } + } + + /// + /// The proving ground (docs/cultmath-interval-ground-cut.md, "Pass 3" and "Pass 4"). Scenarios + /// (a) height fog, (b) inside the fog, (c) Aetheria's wells and (d) r1's uniform slab, each over 64 + /// tiles at N in {1, 4, 8, 16} and warp in {0, D}: the tile march integrates exactly the dense march's + /// cells, so the transmittance agrees within 1e-3 (a difference would be a skipped cell that was not + /// empty), and no (a)-(d) tile reports a body. (e) the void: the grid over hollow radius, camera offset + /// and ramp width (32 tiles at N = 8), LOD-far (32 tiles) and the one shipped void (200 tiles), three + /// marches each, the footprint-aware one taking the analytic body step. Contracts at N = 8, the lower + /// of warp 0 and D: (a) at least 2x combined cost; (c) at least 0.90 of its probe-free ceiling + /// (ruling operator-wells-contract-vs-ceiling: the ceiling is what intervals control, the fog band's + /// cost is the field's); (e) at the shipped void the masked fixed-step march agrees with the reference + /// exactly, the footprint-aware march within 0.02, and it takes at most half the reference's steps per + /// pixel. A shortfall skips naming saving-2x-scenarios and every shortfall; the fields, the grids and + /// the envelope cost are not tuned toward it. /// [Fact] public void IntervalSkipHalvesEvaluations() @@ -1259,7 +1434,7 @@ public sealed class NoiseBoundTests if (key == "d") output.WriteLine($"IV-REPORT (d) N={n} warp={w:R}: ungated dense snoise/ray {3.0 * stats.PerRay(stats.DenseSteps):F2} (three per dense sample, as r1 counted)"); if (n == 8) - ratios[key] = Math.Min(ratios.GetValueOrDefault(key, double.MaxValue), stats.Ratio); + ratios[key] = Math.Min(ratios.GetValueOrDefault(key, double.MaxValue), key == "c" ? stats.Ratio / stats.Ceiling : stats.Ratio); } foreach (var rh in VoidField.HollowRadii) @@ -1270,6 +1445,9 @@ public sealed class NoiseBoundTests ReportVoid($"{field.Name} N=8", Run(field, 8, 32, 0x701D, lod: true)); } + var lodFar = VoidField.LodFar(); + ReportVoid($"{lodFar.Name} N=8 (LOD far)", Run(lodFar, 8, 32, 0x701D, lod: true)); + var shipped = VoidField.Shipped(); var headline = Run(shipped, 8, 200, 0x5417, lod: true); ReportVoid($"{shipped.Name} N=8 (headline, 200 tiles)", headline); @@ -1277,7 +1455,15 @@ public sealed class NoiseBoundTests Assert.True(headline.MaxLodDifference <= 0.02f, $"the footprint-aware march differs from the reference by {headline.MaxLodDifference:R}"); var steps = (double)headline.DenseSteps / headline.LodSteps; - if (ratios["a"] < 2.0 || ratios["c"] < 2.0 || steps < 2.0) - Assert.Skip($"saving-2x-scenarios: at N = 8, (a) saves {ratios["a"]:F2}x and (c) {ratios["c"]:F2}x combined cost (the lower of warp 0 and D); (e) takes {steps:F2}x fewer steps per pixel at the shipped point; the contract is 2x each"); + var shortfalls = new List(); + if (ratios["a"] < 2.0) + shortfalls.Add($"(a) saves {ratios["a"]:F2}x combined cost, under 2x"); + if (ratios["c"] < 0.90) + shortfalls.Add($"(c) reaches {ratios["c"]:F2} of its probe-free ceiling, under 0.90"); + if (steps < 2.0) + shortfalls.Add($"(e) takes {steps:F2}x fewer steps per pixel at the shipped void, under 2x"); + output.WriteLine($"IV-REPORT contracts at N=8 (lower of warp 0 and D): (a) {ratios["a"]:F2}x cost, (c) {ratios["c"]:F2} of ceiling, (e) {steps:F2}x fewer steps/px"); + if (shortfalls.Count > 0) + Assert.Skip("saving-2x-scenarios: " + string.Join("; ", shortfalls)); } } From 3b1b5d358cfc212c51fb044a9977ba7ba4fb79f8 Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Sat, 3 Oct 2026 13:54:36 +0200 Subject: [PATCH 12/14] CultMath intervals r4: LOD agreement against a converged march at the body start, BodyDepth, (c) against an oracle The footprint march's optical depth at the body start is measured against a converged midpoint march (h = 0.5, every cell) and may be at most twice as wrong as the fixed-step march; the end-of-ray assertion goes. BodyDepth owns the body's optical depth and its pin compares depths, not transmittances. (c) asserts probe overhead >= 0.90 and a hard cull floor of 0.5 (ruling wells-overhead-plus-floor); the oracle ceiling, efficiency and cull fraction are printed. WeightBoundEnclosesOneSignedNoise pins VoidField's weight bound. Co-Authored-By: Claude Opus 5.5 --- .../tests/CultMath.Tests/NoiseBoundTests.cs | 217 ++++++++++++++---- 1 file changed, 172 insertions(+), 45 deletions(-) diff --git a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs index 301979f5..2b00ca88 100644 --- a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs @@ -15,9 +15,10 @@ namespace CultMath.Tests; /// Hand mutations of iv_frustum_ball that TileBallEnclosesEveryRaySegment must kill: z0 for z1 in the /// lateral term (z0 * footprintPerDepth); |m_c| for |(m_c, 1)| in the depth term; the warp dropped; the /// rounding widening dropped (which the degenerate family's enclosure check alone must catch). Hand -/// mutations of the test's own fields: w(z0) for w(z1) in VoidField.WeightBound must fail the truncated -/// check of EnvelopeBoundEnclosesDensity on the LOD-far configuration; CarveRadius for HollowRadius in -/// the lo test of VoidField.Compose must fail AnalyticBodyMatchesFineMarch. +/// mutations of the test's own fields and marches: w(z0) for w(z1) in VoidField.WeightBound must fail +/// WeightBoundEnclosesOneSignedNoise; CarveRadius for HollowRadius in the lo test of VoidField.Compose +/// and BodyDepth without |dir| must fail AnalyticBodyMatchesFineMarch; the footprint march's pre-body +/// integration dropped, or LodStep doubled, must fail IntervalSkipHalvesEvaluations's depth agreement. /// public sealed class NoiseBoundTests { @@ -302,6 +303,40 @@ public sealed class NoiseBoundTests Assert.True(partial > 0 && zero > 0, $"octave weights in (0, 1) at {partial} samples and fine weight 0 at {zero}"); } + /// + /// VoidField's octave weight bound over [z0, z1] encloses each octave's weight at every depth of the + /// segment, and, for a one-signed noise interval n = [0.25, 1], the weighted noise w(z) n_s lies in + /// iv_mul(bound, n). 2,000 LOD-far segments with z0 ~ U(1460, 0.95 grid end), z1 ~ U(z0, grid end), + /// 16 depths each including both ends. A bound whose lower end is w(z0) misses w(z1) 0.25 wherever the + /// weight falls across the segment. + /// + [Fact] + public void WeightBoundEnclosesOneSignedNoise() + { + var field = VoidField.LodFar(); + var random = new System.Random(0x3B0D); + var end = field.Grid[^1]; + var n = new float2(0.25f, 1.0f); + for (var i = 0; i < 2000; i++) + { + var z0 = Uniform(random, 1460.0f, 0.95f * end); + var z1 = Uniform(random, z0, end); + var (coarseBound, fineBound) = field.WeightBounds(z0, z1); + for (var s = 0; s < 16; s++) + { + var z = s == 0 ? z0 : s == 15 ? z1 : Uniform(random, z0, z1); + var (coarse, fine) = field.Weights(z); + foreach (var (octave, w, bound) in new[] { ("coarse", coarse, coarseBound), ("fine", fine, fineBound) }) + { + Assert.True(Inside(w, bound), $"{octave} weight {w:R} at depth {z:R} lies outside its bound {bound} over [{z0:R}, {z1:R}]"); + var weighted = iv_mul(bound, n); + foreach (var ns in new[] { 0.25f, Uniform(random, 0.25f, 1.0f), 1.0f }) + Assert.True(Inside(w * ns, weighted), $"{octave} weight {w:R} x noise {ns:R} at depth {z:R} lies outside iv_mul({bound}, {n}) = {weighted} over [{z0:R}, {z1:R}]"); + } + } + } + } + // ---- The Lipschitz constant ---- private static float GradientNorm(float3 p) @@ -991,6 +1026,9 @@ public sealed class NoiseBoundTests /// The coarse and fine octave weights at depth z (1 for the full field). public (float Coarse, float Fine) Weights(float z) => (Weight(CoarseWavelength, z * Theta), Weight(FineWavelength, z * Theta)); + /// The coarse and fine octave weight bounds Bound uses over [z0, z1]. + public (float2 Coarse, float2 Fine) WeightBounds(float z0, float z1) => (WeightBound(CoarseWavelength, z0, z1), WeightBound(FineWavelength, z0, z1)); + private sealed class VoidTile : Tile { public float3 Cavity; // c_v in the camera frame @@ -1195,33 +1233,36 @@ public sealed class NoiseBoundTests return (dense, -1, 0.0f); } - // The body's one closed-form Beer-Lambert step from z to the grid's end, at the density the field's - // pointwise envelope gives at z (one envelope evaluation; no probe in the main pass). - private static float BodyFactor(Scenario field, Tile tile, float2 slope, float z, bool lod, Counts counts) => - exp(-(field.Density(tile, slope, z, lod, counts) * field.Extinction) * ((field.Grid[^1] - z) * sqrt(slope.x * slope.x + slope.y * slope.y + 1.0f))); + // The body's optical depth from z to the grid's end, at the density the field's pointwise envelope gives + // at z (one envelope evaluation; no probe in the main pass). MarchLod takes exp(-BodyDepth) as one step. + private static float BodyDepth(Scenario field, Tile tile, float2 slope, float z, bool lod, Counts counts) => + field.Density(tile, slope, z, lod, counts) * field.Extinction * ((field.Grid[^1] - z) * sqrt(slope.x * slope.x + slope.y * slope.y + 1.0f)); // The footprint-aware march, the void's shipped one: steps of LodStep(z) inside the unmasked runs only, - // octave weights applied; at the body's start, one analytic step to the grid's end, and stop. - private static float MarchLod(VoidField field, Tile tile, float2 slope, List cells, int bodyStart, Counts counts, ref long steps, ref long bodySteps) + // octave weights applied, early-out below cutoff; at the body's start, one analytic step to the grid's + // end, and stop. AtBody is the transmittance when the march reaches the body's start, before that step, + // or the final transmittance when there is no body. + private static (float Transmittance, float AtBody) MarchLod(VoidField field, Tile tile, float2 slope, List cells, int bodyStart, float cutoff, Counts counts, ref long steps, ref long bodySteps) { var grid = field.Grid; var body = bodyStart >= 0 ? grid[bodyStart] : float.PositiveInfinity; var transmittance = 1.0f; var k = 0; - while (k < cells.Count && transmittance >= 0.02f) + while (k < cells.Count && transmittance >= cutoff) { var z = grid[cells[k]]; var end = cells[k] + 1; for (k++; k < cells.Count && cells[k] == end; k++) end++; - while (z < grid[end] && transmittance >= 0.02f) + while (z < grid[end] && transmittance >= cutoff) { if (z >= body) { - transmittance *= BodyFactor(field, tile, slope, z, lod: true, counts); + var atBody = transmittance; + transmittance *= exp(-BodyDepth(field, tile, slope, z, lod: true, counts)); steps++; bodySteps++; - return transmittance; + return (transmittance, atBody); } var next = min(min(z + field.LodStep(z), grid[end]), z < body ? body : grid[end]); @@ -1231,7 +1272,26 @@ public sealed class NoiseBoundTests } } - return transmittance; + return (transmittance, transmittance); + } + + // The accuracy oracle before the body: the optical depth over every grid cell below `stop` (never the + // mask's, so it trusts no pre-pass), by midpoint quadrature at step h, accumulated in double. Not counted. + private static double ConvergedDepth(Scenario field, Tile tile, float2 slope, int stop, float h, bool lod) + { + var grid = field.Grid; + var counts = new Counts(); + var dir = Math.Sqrt(slope.x * (double)slope.x + slope.y * (double)slope.y + 1.0); + var tau = 0.0; + for (var c = 0; c < stop; c++) + { + var samples = Math.Max(1, (int)Math.Ceiling((grid[c + 1] - grid[c]) / h)); + var dz = (grid[c + 1] - grid[c]) / (double)samples; + for (var s = 0; s < samples; s++) + tau += field.Density(tile, slope, (float)(grid[c] + (s + 0.5) * dz), lod, counts) * (double)field.Extinction * dz * dir; + } + + return tau; } internal sealed class RunStats @@ -1240,6 +1300,7 @@ public sealed class NoiseBoundTests public readonly Counts Masked = new(); public readonly Counts Probes = new(); public readonly Counts Lod = new(); + public readonly Counts Oracle = new(); public long Rays; public long Tiles; public long DenseSteps; @@ -1250,27 +1311,46 @@ public sealed class NoiseBoundTests public long Unfinished; public float MaxDifference; public float MaxLodDifference; + public float MaxOracleDifference; + public long OracleEmpty; + public long OracleEmptySkipped; + public double MaxLodDepthError; + public double MaxRefDepthError; public double PerRay(double value) => value / Rays; // Dense cost over the tile march's: each ray's own dense cells plus the tile's probes, shared. public double Ratio => Dense.Cost / (Masked.Cost + Probes.Cost); - // What no amortization can beat: the dense cost over the cost of the dense cells alone. - public double Ceiling => Dense.Cost / Masked.Cost; - public double LodRatio => Dense.Cost / (Lod.Cost + Probes.Cost); + + // The oracle mask (the cells some ray of the tile samples nonzero) is a ceiling no pre-pass can move. + public double OracleCeiling => Dense.Cost / Oracle.Cost; + + public double Efficiency => Oracle.Cost / (Masked.Cost + Probes.Cost); + + // The tile march over the tile march plus its probes: what the probes cost, not what they cull. + public double ProbeOverhead => Masked.Cost / (Masked.Cost + Probes.Cost); + + // The oracle-empty cells the pre-pass proves empty; a pre-pass that skips nothing scores 0. + public double CullFraction => (double)OracleEmptySkipped / OracleEmpty; } /// - /// Draws `tiles` tiles, runs the pre-pass once per tile, and marches each of the N^2 rays (one per pixel, - /// jittered inside it) densely over every cell and masked over the tile's dense cells, both with the - /// early-out at 0.02 and the same Integrate; with lod, the footprint-aware march as well. + /// Draws `tiles` tiles, runs the pre-pass once per tile, draws the tile's N^2 ray slopes (one per pixel, + /// jittered inside it), builds the oracle mask (the cells whose mid-depth density is nonzero for at least + /// one of those slopes), and marches each ray densely over every cell, masked over the tile's dense cells + /// and over the oracle cells, all with the early-out at 0.02 and the same Integrate. With lod, the + /// footprint-aware march as well, and the agreement at the body start (or the grid's end): the optical + /// depth of a cutoff-0 footprint march and of the fixed-step cells, each against ConvergedDepth at + /// h = 0.5. The oracle mask and the agreement marches count into throwaway Counts, the oracle march into + /// its own, so the cost rows do not move. /// private static RunStats Run(Scenario field, int n, int tiles, int seed, bool lod = false) { var random = new System.Random(seed); var stats = new RunStats(); + var scratch = new Counts(); var grid = field.Grid; for (var t = 0; t < tiles; t++) { @@ -1279,24 +1359,58 @@ public sealed class NoiseBoundTests Assert.True(field is VoidField || bodyStart < 0, $"{field.Name}: a body reported at cell {bodyStart}; only the void's density can be provably constant"); stats.Tiles++; stats.BodyTiles += bodyStart >= 0 ? 1 : 0; + var slopes = new List(); for (var j = 0; j < n; j++) for (var i = 0; i < n; i++) + slopes.Add(tile.PixelSlope(i, j, random.NextSingle(), random.NextSingle())); + var inMask = new HashSet(cells); + var oracle = new List(); + for (var c = 0; c < grid.Length - 1; c++) + { + var mid = (grid[c] + grid[c + 1]) * 0.5f; + if (slopes.Any(s => field.Density(tile, s, mid, false, scratch) > 0.0f)) + { + oracle.Add(c); + } + else + { + stats.OracleEmpty++; + stats.OracleEmptySkipped += inMask.Contains(c) ? 0 : 1; + } + } + + var stop = bodyStart >= 0 ? bodyStart : grid.Length - 1; + foreach (var slope in slopes) { - var slope = tile.PixelSlope(i, j, random.NextSingle(), random.NextSingle()); var dense = 1.0f; for (var c = 0; c < grid.Length - 1 && dense >= 0.02f; c++, stats.DenseSteps++) dense *= Integrate(field, tile, slope, grid[c], grid[c + 1], lod: false, stats.Dense); var masked = 1.0f; for (var c = 0; c < cells.Count && masked >= 0.02f; c++, stats.MaskedSteps++) masked *= Integrate(field, tile, slope, grid[cells[c]], grid[cells[c] + 1], lod: false, stats.Masked); + var exact = 1.0f; + for (var c = 0; c < oracle.Count && exact >= 0.02f; c++) + exact *= Integrate(field, tile, slope, grid[oracle[c]], grid[oracle[c] + 1], lod: false, stats.Oracle); stats.Rays++; stats.MaxDifference = MathF.Max(stats.MaxDifference, MathF.Abs(dense - masked)); + stats.MaxOracleDifference = MathF.Max(stats.MaxOracleDifference, MathF.Abs(dense - exact)); if (dense >= 0.02f) stats.Unfinished++; if (lod) { - var marched = MarchLod((VoidField)field, tile, slope, cells, bodyStart, stats.Lod, ref stats.LodSteps, ref stats.BodySteps); + var (marched, _) = MarchLod((VoidField)field, tile, slope, cells, bodyStart, 0.02f, stats.Lod, ref stats.LodSteps, ref stats.BodySteps); stats.MaxLodDifference = MathF.Max(stats.MaxLodDifference, MathF.Abs(dense - marched)); + long ignoredSteps = 0, ignoredBody = 0; + var (_, atBody) = MarchLod((VoidField)field, tile, slope, cells, bodyStart, 0.0f, scratch, ref ignoredSteps, ref ignoredBody); + var reference = 1.0; + for (var c = 0; c < stop; c++) + reference *= Integrate(field, tile, slope, grid[c], grid[c + 1], lod: false, scratch); + var tauConv = ConvergedDepth(field, tile, slope, stop, 0.5f, lod: true); + var tauLod = -Math.Log(Math.Max(atBody, 1.0e-30f)); + var tauRef = -Math.Log(Math.Max(reference, 1.0e-300)); + var scale = Math.Max(tauConv, 1.0e-3); + stats.MaxLodDepthError = Math.Max(stats.MaxLodDepthError, Math.Abs(tauLod - tauConv) / scale); + stats.MaxRefDepthError = Math.Max(stats.MaxRefDepthError, Math.Abs(tauRef - tauConv) / scale); } } } @@ -1308,14 +1422,16 @@ public sealed class NoiseBoundTests $"IV-REPORT {label}: dense snoise/ray {s.PerRay(s.Dense.Snoise):F2} env/ray {s.PerRay(s.Dense.Envelope):F2} (samples/ray {s.PerRay(s.DenseSteps):F2}); " + $"tile snoise/ray {s.PerRay(s.Masked.Snoise + (double)s.Probes.Snoise):F2} env/ray {s.PerRay(s.Masked.Envelope + (double)s.Probes.Envelope):F2} " + $"(probes/tile env {(double)s.Probes.Envelope / s.Tiles:F1} snoise {(double)s.Probes.Snoise / s.Tiles:F2}; dense cells/ray {s.PerRay(s.MaskedSteps):F2}); " - + $"cost ratio {s.Ratio:F2}x, ceiling {s.Ceiling:F2}x; max |dT| {s.MaxDifference:R}; unfinished rays {s.Unfinished}"); + + $"cost ratio {s.Ratio:F2}x, oracle ceiling {s.OracleCeiling:F2}x, efficiency {s.Efficiency:F3}, cull fraction {s.CullFraction:F3} ({s.OracleEmptySkipped}/{s.OracleEmpty}), " + + $"probe overhead {s.ProbeOverhead:F3}; max |dT| {s.MaxDifference:R}, oracle {s.MaxOracleDifference:R}; unfinished rays {s.Unfinished}"); private void ReportVoid(string label, RunStats s) => output.WriteLine( $"IV-REPORT {label}: steps/px ref {s.PerRay(s.DenseSteps):F2}, masked {s.PerRay(s.MaskedSteps):F2}, LOD {s.PerRay(s.LodSteps):F2} ({(double)s.DenseSteps / s.LodSteps:F2}x fewer); " + $"snoise/px ref {s.PerRay(s.Dense.Snoise):F2}, masked {s.PerRay(s.Masked.Snoise):F2}, LOD {s.PerRay(s.Lod.Snoise):F2}; " + $"cost/px ref {s.PerRay(s.Dense.Cost):F2}, masked {s.PerRay(s.Masked.Cost + s.Probes.Cost):F2}, LOD {s.PerRay(s.Lod.Cost + s.Probes.Cost):F2} incl. probes (ref/LOD {s.LodRatio:F2}x, ref/masked {s.Ratio:F2}x); " + $"probes/tile env {(double)s.Probes.Envelope / s.Tiles:F1} snoise {(double)s.Probes.Snoise / s.Tiles:F2}; body tiles {s.BodyTiles}/{s.Tiles}, analytic steps/px {s.PerRay(s.BodySteps):F2}; " - + $"max |dT| masked {s.MaxDifference:R}, LOD {s.MaxLodDifference:R}; unfinished rays {s.Unfinished}"); + + $"depth error at the body start vs converged: LOD {s.MaxLodDepthError:G4}, fixed-step {s.MaxRefDepthError:G4} (ratio {s.MaxLodDepthError / s.MaxRefDepthError:F2}); " + + $"end-of-ray max |dT| masked {s.MaxDifference:R}, LOD {s.MaxLodDifference:R}; unfinished rays {s.Unfinished}"); /// /// The analytic body's soundness pin (ruling operator-analytic-wall-body). Over 2,000 seeded tiles of @@ -1323,9 +1439,9 @@ public sealed class NoiseBoundTests /// body start b: (1) the gated bound over [grid[b], grid end] is a positive point equal to the reported /// density, every cell from b to the end is in the mask, and at 64 random (pixel, sub-pixel, cell >= b, /// depth in cell) points the full and the truncated Density and the ungated Exact field all equal it - /// exactly, no tolerance; (2) the one closed-form factor exp(-density ext L) over the body equals the - /// fixed-step product of the body cells' Integrate factors (no early-out, same cells) within 1e-4 - /// relative. (3) No tile of (a)-(d), at either warp, reports a body: their bounds are never a positive + /// exactly, no tolerance; (2) BodyDepth, the optical depth MarchLod's one closed-form step uses, equals + /// the summed fixed-step depth -ln Integrate of the body cells (no early-out, same cells) within 1e-4 + /// relative, compared in the depth domain so LOD-far's deep bodies do not underflow. (3) No tile of (a)-(d), at either warp, reports a body: their bounds are never a positive /// point. /// [Fact] @@ -1369,12 +1485,11 @@ public sealed class NoiseBoundTests } var ray = tile.PixelSlope(random.Next(tile.N), random.Next(tile.N), random.NextSingle(), random.NextSingle()); - var product = 1.0; + var depth = BodyDepth(field, tile, ray, grid[b], lod: true, counts); + var sum = 0.0; for (var c = b; c < last; c++) - product *= Integrate(field, tile, ray, grid[c], grid[c + 1], lod: true, counts); - var tau = density * field.Extinction * ((grid[last] - grid[b]) * sqrt(ray.x * ray.x + ray.y * ray.y + 1.0f)); - var analytic = Math.Exp(-(double)tau); - Assert.True(Math.Abs(analytic - product) <= 1.0e-4 * product, $"{field.Name}: the analytic body factor {analytic:R} differs from the fixed-step product {product:R} over cells {b}..{last}"); + sum += -Math.Log(Integrate(field, tile, ray, grid[c], grid[c + 1], lod: true, counts)); + Assert.True(Math.Abs(depth - sum) <= 1.0e-4 * sum, $"{field.Name}: the body's depth from cell {b}, BodyDepth {depth:R}, differs from the fixed-step cells' summed depth {sum:R} over cells {b}..{last}"); } output.WriteLine("IV-REPORT analytic body: tiles with a body " + string.Join(", ", bodies.OrderBy(kv => kv.Key).Select(kv => $"{kv.Key} {kv.Value}"))); @@ -1402,20 +1517,25 @@ public sealed class NoiseBoundTests /// (a) height fog, (b) inside the fog, (c) Aetheria's wells and (d) r1's uniform slab, each over 64 /// tiles at N in {1, 4, 8, 16} and warp in {0, D}: the tile march integrates exactly the dense march's /// cells, so the transmittance agrees within 1e-3 (a difference would be a skipped cell that was not - /// empty), and no (a)-(d) tile reports a body. (e) the void: the grid over hollow radius, camera offset - /// and ramp width (32 tiles at N = 8), LOD-far (32 tiles) and the one shipped void (200 tiles), three - /// marches each, the footprint-aware one taking the analytic body step. Contracts at N = 8, the lower - /// of warp 0 and D: (a) at least 2x combined cost; (c) at least 0.90 of its probe-free ceiling - /// (ruling operator-wells-contract-vs-ceiling: the ceiling is what intervals control, the fog band's - /// cost is the field's); (e) at the shipped void the masked fixed-step march agrees with the reference - /// exactly, the footprint-aware march within 0.02, and it takes at most half the reference's steps per - /// pixel. A shortfall skips naming saving-2x-scenarios and every shortfall; the fields, the grids and - /// the envelope cost are not tuned toward it. + /// empty), the oracle march (each ray over the oracle mask: the cells whose mid-depth density is nonzero + /// for at least one of the tile's rays) agrees with the dense march exactly, and no (a)-(d) tile reports + /// a body. (e) the void: the grid over hollow radius, camera offset and ramp width (32 tiles at N = 8), + /// LOD-far (32 tiles) and the one shipped void (200 tiles), three marches each, the footprint-aware one + /// taking the analytic body step. Contracts at N = 8, the lower of warp 0 and D: (a) at least 2x + /// combined cost; (c) (ruling wells-overhead-plus-floor) probe overhead, the tile march over the tile + /// march plus its probes, at least 0.90, and, hard, the pre-pass proves at least half of the + /// oracle-empty cells empty; the oracle ceiling, the efficiency against it and the cull fraction are + /// printed; (e) at the shipped void the masked fixed-step march agrees with the reference exactly, the + /// footprint-aware march's optical depth at the body start is off a converged march (h = 0.5, every + /// cell) by at most twice the fixed-step march's error, and it takes at most half the reference's + /// steps per pixel. A shortfall skips naming saving-2x-scenarios and every shortfall; the fields, the + /// grids and the envelope cost are not tuned toward it. /// [Fact] public void IntervalSkipHalvesEvaluations() { var ratios = new Dictionary(); + var cullFloor = double.MaxValue; var scenarios = new (string Key, float Warp, Func Make)[] { ("a", FogField.AetheriaWarp, FogField.HeightFog), @@ -1431,10 +1551,16 @@ public sealed class NoiseBoundTests var stats = Run(field, n, 64, 0x5A7E + n); Report($"{field.Name} N={n} warp={w:R}", stats); Assert.True(stats.MaxDifference <= 1.0e-3f, $"{field.Name} N={n} warp={w:R}: transmittance differs by {stats.MaxDifference:R}"); + Assert.True(stats.MaxOracleDifference == 0.0f, $"{field.Name} N={n} warp={w:R}: the oracle march differs from the dense march by {stats.MaxOracleDifference:R}"); if (key == "d") output.WriteLine($"IV-REPORT (d) N={n} warp={w:R}: ungated dense snoise/ray {3.0 * stats.PerRay(stats.DenseSteps):F2} (three per dense sample, as r1 counted)"); if (n == 8) - ratios[key] = Math.Min(ratios.GetValueOrDefault(key, double.MaxValue), key == "c" ? stats.Ratio / stats.Ceiling : stats.Ratio); + ratios[key] = Math.Min(ratios.GetValueOrDefault(key, double.MaxValue), key == "c" ? stats.ProbeOverhead : stats.Ratio); + if (n == 8 && key == "c") + { + Assert.True(stats.CullFraction >= 0.5, $"(c) warp={w:R}: the pre-pass proves {stats.CullFraction:F3} of the oracle-empty cells empty, under the 0.5 floor (ruling wells-overhead-plus-floor)"); + cullFloor = Math.Min(cullFloor, stats.CullFraction); + } } foreach (var rh in VoidField.HollowRadii) @@ -1452,17 +1578,18 @@ public sealed class NoiseBoundTests var headline = Run(shipped, 8, 200, 0x5417, lod: true); ReportVoid($"{shipped.Name} N=8 (headline, 200 tiles)", headline); Assert.True(headline.MaxDifference == 0.0f, $"the masked fixed-step march differs from the reference by {headline.MaxDifference:R}"); - Assert.True(headline.MaxLodDifference <= 0.02f, $"the footprint-aware march differs from the reference by {headline.MaxLodDifference:R}"); + Assert.True(headline.MaxLodDepthError <= 2.0 * headline.MaxRefDepthError, $"at the body start the footprint-aware march's optical depth is off the converged march's by {headline.MaxLodDepthError:G6} relative, more than twice the fixed-step march's {headline.MaxRefDepthError:G6}"); var steps = (double)headline.DenseSteps / headline.LodSteps; var shortfalls = new List(); if (ratios["a"] < 2.0) shortfalls.Add($"(a) saves {ratios["a"]:F2}x combined cost, under 2x"); if (ratios["c"] < 0.90) - shortfalls.Add($"(c) reaches {ratios["c"]:F2} of its probe-free ceiling, under 0.90"); + shortfalls.Add($"(c) probe overhead {ratios["c"]:F3}, under 0.90"); if (steps < 2.0) shortfalls.Add($"(e) takes {steps:F2}x fewer steps per pixel at the shipped void, under 2x"); - output.WriteLine($"IV-REPORT contracts at N=8 (lower of warp 0 and D): (a) {ratios["a"]:F2}x cost, (c) {ratios["c"]:F2} of ceiling, (e) {steps:F2}x fewer steps/px"); + output.WriteLine($"IV-REPORT contracts at N=8 (lower of warp 0 and D): (a) {ratios["a"]:F2}x cost; (c) probe overhead {ratios["c"]:F3} (at least 0.90), cull fraction {cullFloor:F3} (at least 0.5, hard); " + + $"(e) {steps:F2}x fewer steps/px, depth error at the body start LOD {headline.MaxLodDepthError:G4} vs fixed-step {headline.MaxRefDepthError:G4} (at most 2x)"); if (shortfalls.Count > 0) Assert.Skip("saving-2x-scenarios: " + string.Join("; ", shortfalls)); } From eb8bcc55952ce1e8ef04e84afae62c82760627fc Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Sat, 3 Oct 2026 13:58:10 +0200 Subject: [PATCH 13/14] CultMath design.md: the saving section from r4's run; the unprinted 6e-4 sentence goes Every figure is from r4's IV-REPORT lines: (a)-(d) with (c)'s oracle ceiling, efficiency, cull fraction and probe overhead; (e) at the shipped void with the two depth errors against the converged march; the analytic body and its pin; the [w(z1), 1] weight-bound deviation; the contracts as asserted, (c) in the words of ruling wells-overhead-plus-floor. Co-Authored-By: Claude Opus 5.5 --- packages/cultmath/docs/design.md | 108 +++++++++++++++++++++++-------- 1 file changed, 80 insertions(+), 28 deletions(-) diff --git a/packages/cultmath/docs/design.md b/packages/cultmath/docs/design.md index 77688dc5..a8adb795 100644 --- a/packages/cultmath/docs/design.md +++ b/packages/cultmath/docs/design.md @@ -419,8 +419,7 @@ axis length at `3u`, the product, `z1 footprintPerDepth`, the sum, `+ warp`), so the exact radius is at most `(1 + 7u)` times the float one; the widening rounds once more. `2^-20 = 16u` covers `8u (|c|_1 + radius)` twice over, and fma contraction only removes roundings, so HLSL and GLSL compilers are covered by -the same bound. At the shipped scales (`|c|_1` near 600, `r` near 5) it adds -`6e-4`. The caller rotates and translates the centre into world space and scales +the same bound. The caller rotates and translates the centre into world space and scales centre and radius by its noise frequency; the rounding of that transform is the caller's. `NoiseBoundTests.TileBallEnclosesEveryRaySegment` checks 2,000 seeded tiles x 64 points of the r2 domain and 2,000 tiles x 8 corner points of each of @@ -476,39 +475,92 @@ The saving is measured where there is empty space to find. `NoiseBoundTests.IntervalSkipHalvesEvaluations` draws screen tiles, each with its own camera, and runs one pre-pass per tile over the dense grid's cell ranges: it starts with the whole grid, skips and doubles a range it proves empty, halves -one it cannot, and runs to the grid's end. Each of the tile's `N^2` rays is then -marched twice with the same integrator and the same early-out at transmittance -0.02: once over every cell, once over the tile's unproven cells only. The two -integrate the same cells, so their transmittance agrees exactly in every -configuration. Cost is counted in `snoise` evaluations, with an envelope +one it cannot, and runs to the grid's end. It then draws the tile's `N^2` ray +slopes, one per pixel, jittered inside it. Each ray is marched three times with +the same integrator and the same early-out at transmittance 0.02: over every +cell (dense), over the tile's unproven cells (the tile march), and over the +oracle mask, the cells whose mid-depth density is nonzero for at least one of +the tile's rays. The tile march and the oracle march integrate the dense march's +nonzero cells, so their transmittance agrees with it exactly in every +configuration. The oracle is a ceiling the pre-pass cannot move: dense cost over +oracle cost. Against it the test prints the efficiency (oracle cost over the +tile march's cost with its probes) and the cull fraction (the oracle-empty cells +the pre-pass proves empty). The probe overhead is the tile march's cost over that +cost plus the probes'. Cost is counted in `snoise` evaluations, with an envelope evaluation or probe weighted 0.2. That weight is an estimate, and both counts are -printed. The flow is not counted. At `N = 8` with the warp on, in Aetheria's units -on its 256-cell quadratic grid: +printed. The flow is not counted. At `N = 8`, in Aetheria's units on its +256-cell quadratic grid out to the Main Camera's far plane, 2048: - (a) Height fog, camera above the safety band, rays level and up: the whole grid is one envelope probe per tile, against 256 envelope evaluations per ray, - 16384x. -- (b) Inside the fog: 0.99x. Every cell is fog, and the pre-pass probes each one. -- (c) The zone bowl and four wells, the camera gazing across: 1.41x, and 3.32x - with no warp. The pre-pass cuts envelope evaluations from 212.93 to 37.12 per - ray. Almost all of the dense march's 78.67 `snoise` per ray fall in the noise - band, where the rays need them (ceiling 1.42x). The warp bound `D = 60` is what - leaves the noise ball unable to prove anything at this noise scale. + 16384x. Every cell is oracle-empty and the pre-pass proves them all. +- (b) Inside the fog: 1.15x, and 0.99x with the warp. Most cells are fog, and + the pre-pass probes each one. +- (c) The zone bowl and four wells, the camera gazing across. With no warp: + 2.74x, oracle ceiling 7.26x, efficiency 0.378, cull fraction 0.952, probe + overhead 0.961. With the warp `D = 60`: 1.39x, oracle ceiling 7.17x, + efficiency 0.194, cull fraction 0.842, probe overhead 0.993. With the warp the + pre-pass cuts envelope evaluations from 177.46 to 33.49 per ray, but the + oracle march costs 0.194 of the tile march, so about four fifths of the tile + march's cost goes to cells no ray samples nonzero: the warp bound swamps the + noise ball, so the ball cannot prove those cells empty. That gap belongs to the tape target's affine forms, not to a + better pre-pass. - (d) r1's uniform slab, which has no envelope: 1.23x, and 1.48x with no warp. -- (e) The void at the shipped point (`Rh = 198`, `S = 50`, the low camera, - 200 tiles), with three marches. The fixed-step reference takes 36.34 steps per +- (e) The void at the shipped point (ruling `operator-shipped-void`): `Rh = 198`, + `S = 50`, `Rc = Rh + S = 248`, `e = -ln 1.5 / ln(1 - (Rh / Rc)^2)`, the + shipped camera, 200 tiles. The fixed-step reference takes 36.34 steps per pixel. The same cells masked by the pre-pass take 20.13, with identical - transmittance. The footprint-aware march takes 16.45 inside the unmasked - cells, within 0.008339431 of the reference's transmittance: 2.21x fewer steps. - `snoise` per pixel goes from 32.98 to 26.56, and combined cost improves 1.34x. - Over the grid of hollow radius, camera offset and ramp width the steps ratio - runs from 1.57x (`Rh = 100`, low, `S = 150`) to 12.40x (`Rh = 1000`, centred, - `S = 10`), and the cost ratio from 1.12x to 3.47x. + transmittance. The footprint-aware march takes 13.54: 2.68x fewer. `snoise` + per pixel goes from 32.98 to 26.69, and combined cost from 40.25 to 29.52 with + the probes, 1.36x. Every one of the 200 tiles proves a body. Before the body + both marches are measured against a converged one (midpoint quadrature at + `h = 0.5` over every cell, trusting no pre-pass): the relative optical-depth + error at the body start is at most 0.0608 for the fixed-step march and 0.08639 + for the footprint march, 1.42x. Over the grid of hollow radius, camera offset + and ramp width the steps ratio runs from 1.57x (`Rh = 100`, low, `S = 150`) to + 27.45x (`Rh = 1000`, centred, `S = 10`), the cost ratio from 1.26x to 3.89x, + and the depth-error ratio from 0.85 to 4.72 (`Rh = 200`, centred, `S = 150`, + where the errors are 0.01504 and 0.003188). LOD-far takes 14.19x fewer steps; + its depth errors, 0.2849 and 0.205, are not comparable, because its + fixed-step reference keeps the octaves the footprint march truncates. -The contract is 2x on combined cost for (a) and (c) and 2x fewer steps per pixel -for (e). (a) and (e) meet it. (c) does not, so the test skips naming -`saving-2x-scenarios`. The fields, the grids and the 0.2 weight are not tuned -toward the contract. +The analytic body (ruling `operator-analytic-wall-body`). When a single dense +cell's own bound is a positive point, the pre-pass spends one more probe over +the rest of the grid. If that bound is a positive point too, every point of +every ray of the tile from there on has exactly that density, so the footprint +march takes one closed-form step, `exp(-BodyDepth)`, from the body start to the +grid's end and stops. `BodyDepth` is density x extinction x length x `|dir|`. +`NoiseBoundTests.AnalyticBodyMatchesFineMarch` pins it: the bound, the mask and +the pointwise densities agree exactly at the body, and `BodyDepth` equals the +summed fixed-step depth of the body cells within 1e-4 relative. The comparison +is in the depth domain, so LOD-far's deep bodies do not underflow. + +The footprint-truncated void weights its octaves by `w_i(z)`, decreasing in `z`. +Its bound takes each weight over `[w_i(z1), 1]` rather than +`[w_i(z1), w_i(z0)]`, because one bound serves both the full field (`w = 1`) and +the truncated one. `[w_i(z1), 1]` contains both, and `iv_mul` is monotone in its +interval argument, so the composed interval encloses both fields. +`NoiseBoundTests.WeightBoundEnclosesOneSignedNoise` checks it on LOD-far over +2,000 segments: each weight lies in its bound, and `w(z) n` lies in +`iv_mul(bound, n)` for a one-signed `n`. A lower end of `w(z0)` fails it. + +The contracts, at `N = 8`, the lower of warp 0 and `D`: + +- (a) at least 2x combined cost; +- (c) (ruling `wells-overhead-plus-floor`): probe overhead, the tile march over + the tile march plus its probes, at least 0.90; and, hard, the pre-pass proves + at least half of the oracle-empty cells empty, so a pre-pass that skips + nothing fails. The oracle ceiling, the efficiency against it and the cull + fraction are printed, not asserted; +- (e) at the shipped void: the masked fixed-step march agrees with the reference + exactly; the footprint march's optical depth at the body start is off the + converged march's by at most twice the fixed-step march's error; and it takes + at most half the reference's steps per pixel. + +This run meets every contract: (a) 16384x, (c) probe overhead 0.961 and cull +fraction 0.842, (e) 2.68x fewer steps and a depth error 1.42x the fixed-step +march's. The fields, the grids and the 0.2 weight are not tuned toward the +contracts. Consumers: the gamecult.org ground shader (through the GLSL lowering) and Aetheria's nebula raymarch at its CultMath pin bump. From 2c8585bdb1686cd8e9e997c98a359d66c2a40c7b Mon Sep 17 00:00:00 2001 From: Meta Vens Date: Sat, 3 Oct 2026 14:12:00 +0200 Subject: [PATCH 14/14] CultMath design.md: rewrap the (c) paragraph Co-Authored-By: Claude Opus 5.5 --- packages/cultmath/docs/design.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/packages/cultmath/docs/design.md b/packages/cultmath/docs/design.md index a8adb795..d102353b 100644 --- a/packages/cultmath/docs/design.md +++ b/packages/cultmath/docs/design.md @@ -503,8 +503,8 @@ printed. The flow is not counted. At `N = 8`, in Aetheria's units on its pre-pass cuts envelope evaluations from 177.46 to 33.49 per ray, but the oracle march costs 0.194 of the tile march, so about four fifths of the tile march's cost goes to cells no ray samples nonzero: the warp bound swamps the - noise ball, so the ball cannot prove those cells empty. That gap belongs to the tape target's affine forms, not to a - better pre-pass. + noise ball, so the ball cannot prove those cells empty. That gap belongs to + the tape target's affine forms, not to a better pre-pass. - (d) r1's uniform slab, which has no envelope: 1.23x, and 1.48x with no warp. - (e) The void at the shipped point (ruling `operator-shipped-void`): `Rh = 198`, `S = 50`, `Rc = Rh + S = 248`, `e = -ln 1.5 / ln(1 - (Rh / Rc)^2)`, the