diff --git a/packages/cultmath/shaders/CultMath.Affine.hlsl b/packages/cultmath/shaders/CultMath.Affine.hlsl index f165148b..8727081e 100644 --- a/packages/cultmath/shaders/CultMath.Affine.hlsl +++ b/packages/cultmath/shaders/CultMath.Affine.hlsl @@ -12,6 +12,10 @@ // math.Affine.cs (SNOISE_HESSIAN). static const float CULTMATH_SNOISE_HESSIAN = 56.050385; +// An upper bound on the operator norm of the third-derivative tensor of cultmath_snoise(float3), the +// next order after CULTMATH_SNOISE_HESSIAN; provenance in math.Span.cs (SNOISE_THIRD). +static const float CULTMATH_SNOISE_THIRD = 1.0; + float3 cultmath_af_point(float value) { return float3(value, 0.0, 0.0); } float3 cultmath_af_symbol(float x0, float a) { return float3(x0, a, 0.0); } diff --git a/packages/cultmath/shaders/CultMath.hlsl b/packages/cultmath/shaders/CultMath.hlsl index c29fd603..18cd1fc9 100644 --- a/packages/cultmath/shaders/CultMath.hlsl +++ b/packages/cultmath/shaders/CultMath.hlsl @@ -228,6 +228,96 @@ float4 cultmath_snoise_grad(float3 value) return float4(grad, value2); } +// Analytic Hessian of cultmath_snoise(float3); see the C# math.snoise_hessian comment for the +// derivation. Every local up through p0..p3 and x0..x3 is cultmath_snoise_grad's own. A float3x3 built +// from nine scalars reads as rows here and as columns in GLSL's mat3: the matrix is symmetric. +float3x3 cultmath_snoise_hessian(float3 value) +{ + float2 c = float2(1.0 / 6.0, 1.0 / 3.0); + float3 i = floor(value + dot(value, c.yyy)); + float3 x0 = value - i + dot(i, c.xxx); + float3 g = float3(step(x0.y, x0.x), step(x0.z, x0.y), 1.0 - step(x0.z, x0.x)); + float3 l = 1.0 - g; + float3 i1 = min(g.xyz, l.zxy); + float3 i2 = max(g.xyz, l.zxy); + float3 x1 = x0 - i1 + c.xxx; + float3 x2 = x0 - i2 + c.yyy; + float3 x3 = x0 - 0.5; + + i = cultmath_snoise_mod289(i); + float4 p = cultmath_snoise_permute(cultmath_snoise_permute(cultmath_snoise_permute( + i.z + float4(0.0, i1.z, i2.z, 1.0)) + i.y + float4(0.0, i1.y, i2.y, 1.0)) + i.x + float4(0.0, i1.x, i2.x, 1.0)); + + const float n_ = 0.142857142857; + float3 ns = float3(2.0 * n_, 0.5 * n_ - 1.0, n_); + float4 j = p - 49.0 * floor(p * ns.z * ns.z); + float4 x_ = floor(j * ns.z); + float4 y_ = floor(j - 7.0 * x_); + float4 x = x_ * ns.x + ns.yyyy; + float4 y = y_ * ns.x + ns.yyyy; + float4 h = 1.0 - abs(x) - abs(y); + + float4 b0 = float4(x.xy, y.xy); + float4 b1 = float4(x.zw, y.zw); + float4 s0 = floor(b0) * 2.0 + 1.0; + float4 s1 = floor(b1) * 2.0 + 1.0; + float4 sh = -step(h, float4(0.0, 0.0, 0.0, 0.0)); + float4 a0 = b0.xzyw + s0.xzyw * sh.xxyy; + float4 a1 = b1.xzyw + s1.xzyw * sh.zzww; + + float3 p0 = float3(a0.xy, h.x); + float3 p1 = float3(a0.zw, h.y); + float3 p2 = float3(a1.xy, h.z); + float3 p3 = float3(a1.zw, h.w); + float4 norm = 1.79284291400159 - 0.85373472095314 * float4(dot(p0, p0), dot(p1, p1), dot(p2, p2), dot(p3, p3)); + p0 *= norm.x; + p1 *= norm.y; + p2 *= norm.z; + p3 *= norm.w; + + float4 m0 = max(0.5 - float4(dot(x0, x0), dot(x1, x1), dot(x2, x2), dot(x3, x3)), 0.0); + float4 m2 = m0 * m0; + float4 m3 = m2 * m0; + float4 px = float4(dot(p0, x0), dot(p1, x1), dot(p2, x2), dot(p3, x3)); + float4 wide = 48.0 * m2 * px; + float4 narrow = 8.0 * m3; + + float4 cxx = float4(x0.x, x1.x, x2.x, x3.x); + float4 cxy = float4(x0.y, x1.y, x2.y, x3.y); + float4 cxz = float4(x0.z, x1.z, x2.z, x3.z); + float4 cpx = float4(p0.x, p1.x, p2.x, p3.x); + float4 cpy = float4(p0.y, p1.y, p2.y, p3.y); + float4 cpz = float4(p0.z, p1.z, p2.z, p3.z); + + float hxx = 105.0 * (dot(wide * cxx, cxx) - dot(narrow, 2.0 * cpx * cxx + px)); + float hyy = 105.0 * (dot(wide * cxy, cxy) - dot(narrow, 2.0 * cpy * cxy + px)); + float hzz = 105.0 * (dot(wide * cxz, cxz) - dot(narrow, 2.0 * cpz * cxz + px)); + float hxy = 105.0 * (dot(wide * cxx, cxy) - dot(narrow, cpx * cxy + cpy * cxx)); + float hxz = 105.0 * (dot(wide * cxx, cxz) - dot(narrow, cpx * cxz + cpz * cxx)); + float hyz = 105.0 * (dot(wide * cxy, cxz) - dot(narrow, cpy * cxz + cpz * cxy)); + return float3x3(hxx, hxy, hxz, hxy, hyy, hyz, hxz, hyz, hzz); +} + +// The integral over [t0, t1] of clamp(a + b t, lo, hi); see the C# math.integral_clamp_linear comment. +float cultmath_integral_clamp_linear(float a, float b, float lo, float hi, float t0, float t1) +{ + if (b == 0.0) + return clamp(a, lo, hi) * (t1 - t0); + + lo = min(lo, hi); + float direction = t1 < t0 ? -1.0 : 1.0; + float start = min(t0, t1); + float stop = max(t0, t1); + float atLo = (lo - a) / b; + float atHi = (hi - a) / b; + float enter = clamp(min(atLo, atHi), start, stop); + float leave = clamp(max(atLo, atHi), start, stop); + float before = b > 0.0 ? lo : hi; + float after = b > 0.0 ? hi : lo; + float total = before * (enter - start) + (a + b * (0.5 * (enter + leave))) * (leave - enter) + after * (stop - leave); + return direction * total; +} + // Octave sum of cultmath_snoise_grad (invariant 8); see the C# math.fbm_grad comment for the // per-octave frequency/amplitude derivation and why octaves is clamped. float4 cultmath_fbm_grad(float3 p, int octaves, float lacunarity, float gain) diff --git a/packages/cultmath/src/CultMath/math.Affine.cs b/packages/cultmath/src/CultMath/math.Affine.cs index 3c00e397..8c4ce94a 100644 --- a/packages/cultmath/src/CultMath/math.Affine.cs +++ b/packages/cultmath/src/CultMath/math.Affine.cs @@ -31,7 +31,7 @@ public static partial class math // curvature is realised if the constant is under it. The margin over the function: a scan of the // whole 289^3 hash period found a realised maximum of 50.58, and bounding each corner's gradient // independently over every hash configuration gives 52.22 (Soul pass s1). A change to the snoise - // kernel must re-run MeasureHessian and re-pin it. + // kernel must re-run MeasureHessian and re-pin it. SNOISE_THIRD (math.Span.cs) is the next-order pin. public const float SNOISE_HESSIAN = 56.050385f; /// The point form (x, 0, 0): the value x at every eps. Exact. diff --git a/packages/cultmath/src/CultMath/math.Span.cs b/packages/cultmath/src/CultMath/math.Span.cs new file mode 100644 index 00000000..4adc1139 --- /dev/null +++ b/packages/cultmath/src/CultMath/math.Span.cs @@ -0,0 +1,122 @@ +namespace CultMath; + +public static partial class math +{ + // An upper bound on the operator norm of the third-derivative tensor of snoise(float3): + // max over unit u of |D^3 snoise(p)[u, u, u]|, the next order after SNOISE_LIPSCHITZ (gradient) + // and SNOISE_HESSIAN (curvature). Provenance: NoiseBoundTests.MeasureThird (slow, explicit) takes + // the largest such norm over 1e6 seeded points of [-256, 256]^3 (central differences of + // snoise_hessian at h = 1e-3, then an ascent over unit vectors for the tensor's operator norm), + // refines the 1e4 largest by ascent over the point, and multiplies the maximum by 1.10; it prints + // the witness start. NoiseBoundTests.ThirdConstantPinsSampled pins the value against a fast + // re-measure and the ascent from that witness. Measured on the tied kernel (snoise_grad's order), + // not on any pre-tie splice. Empirical with a margin, not a proof. A change to the snoise kernel + // must re-run MeasureThird and re-pin it. + public const float SNOISE_THIRD = 1.0f; + + // Analytic Hessian of snoise(float3): the symmetric matrix of second partial derivatives, same + // lattice, hash and tie order as snoise_grad, whose kernel it differentiates once more. Every local + // up through p0..p3 and x0..x3 is snoise_grad's own derivation verbatim. Per corner, with + // m = max(0.5 - dot(x,x), 0) and s = dot(p,x), the corner's gradient is -8 m^3 s x + m^4 p, so + // d/dx_j of its i-th component is 48 m^2 s x_i x_j - 8 m^3 (p_j x_i + p_i x_j + s delta_ij) (the + // derivative of m^3 is -6 m^2 x_j, of s is p_j); summed over the four corners and scaled by the + // same 105. A clamped corner (m = 0) contributes nothing, and the m^2 and m^3 factors make the + // contribution vanish continuously as it clamps. Each off-diagonal entry is computed once and + // placed twice, so the matrix is its own transpose bit for bit. A float3x3 is built from its nine + // scalars, which the GLSL lowering reads as columns and HLSL as rows: the matrix is symmetric. + public static float3x3 snoise_hessian(float3 value) + { + const float cx = 1.0f / 6.0f; + const float cy = 1.0f / 3.0f; + var i = floor(value + dot(value, new float3(cy, cy, cy))); + var x0 = value - i + dot(i, new float3(cx, cx, cx)); + var g = new float3(step(x0.y, x0.x), step(x0.z, x0.y), 1.0f - step(x0.z, x0.x)); + var l = 1.0f - g; + var lzxy = new float3(l.z, l.x, l.y); + var i1 = min(g, lzxy); + var i2 = max(g, lzxy); + var x1 = x0 - i1 + cx; + var x2 = x0 - i2 + cy; + var x3 = x0 - 0.5f; + + i = snoise_mod289(i); + var p = snoise_permute(snoise_permute(snoise_permute( + i.z + new float4(0.0f, i1.z, i2.z, 1.0f)) + i.y + new float4(0.0f, i1.y, i2.y, 1.0f)) + i.x + new float4(0.0f, i1.x, i2.x, 1.0f)); + + const float n = 0.142857142857f; + var ns = new float3(2.0f * n, 0.5f * n - 1.0f, n); + var j = p - 49.0f * floor(p * ns.z * ns.z); + var xs = floor(j * ns.z); + var ys = floor(j - 7.0f * xs); + var x = xs * ns.x + ns.y; + var y = ys * ns.x + ns.y; + var h = 1.0f - abs(x) - abs(y); + + var b0 = new float4(x.x, x.y, y.x, y.y); + var b1 = new float4(x.z, x.w, y.z, y.w); + var s0 = floor(b0) * 2.0f + 1.0f; + var s1 = floor(b1) * 2.0f + 1.0f; + var sh = -step(h, 0.0f); + var a0 = new float4(b0.x, b0.z, b0.y, b0.w) + new float4(s0.x, s0.z, s0.y, s0.w) * new float4(sh.x, sh.x, sh.y, sh.y); + var a1 = new float4(b1.x, b1.z, b1.y, b1.w) + new float4(s1.x, s1.z, s1.y, s1.w) * new float4(sh.z, sh.z, sh.w, sh.w); + + var p0 = new float3(a0.x, a0.y, h.x); + var p1 = new float3(a0.z, a0.w, h.y); + var p2 = new float3(a1.x, a1.y, h.z); + var p3 = new float3(a1.z, a1.w, h.w); + var norm = 1.79284291400159f - 0.85373472095314f * new float4(dot(p0, p0), dot(p1, p1), dot(p2, p2), dot(p3, p3)); + p0 *= norm.x; + p1 *= norm.y; + p2 *= norm.z; + p3 *= norm.w; + + var m0 = max(0.5f - new float4(dot(x0, x0), dot(x1, x1), dot(x2, x2), dot(x3, x3)), 0.0f); + var m2 = m0 * m0; + var m3 = m2 * m0; + var px = new float4(dot(p0, x0), dot(p1, x1), dot(p2, x2), dot(p3, x3)); + var wide = 48.0f * m2 * px; + var narrow = 8.0f * m3; + + // One component of every corner's x and p per register. + var cxx = new float4(x0.x, x1.x, x2.x, x3.x); + var cxy = new float4(x0.y, x1.y, x2.y, x3.y); + var cxz = new float4(x0.z, x1.z, x2.z, x3.z); + var cpx = new float4(p0.x, p1.x, p2.x, p3.x); + var cpy = new float4(p0.y, p1.y, p2.y, p3.y); + var cpz = new float4(p0.z, p1.z, p2.z, p3.z); + + var hxx = 105.0f * (dot(wide * cxx, cxx) - dot(narrow, 2.0f * cpx * cxx + px)); + var hyy = 105.0f * (dot(wide * cxy, cxy) - dot(narrow, 2.0f * cpy * cxy + px)); + var hzz = 105.0f * (dot(wide * cxz, cxz) - dot(narrow, 2.0f * cpz * cxz + px)); + var hxy = 105.0f * (dot(wide * cxx, cxy) - dot(narrow, cpx * cxy + cpy * cxx)); + var hxz = 105.0f * (dot(wide * cxx, cxz) - dot(narrow, cpx * cxz + cpz * cxx)); + var hyz = 105.0f * (dot(wide * cxy, cxz) - dot(narrow, cpy * cxz + cpz * cxy)); + return new float3x3(hxx, hxy, hxz, hxy, hyy, hyz, hxz, hyz, hzz); + } + + // The integral over [t0, t1] of clamp(a + b t, lo, hi) dt, in closed form. The line meets lo at + // (lo - a) / b and hi at (hi - a) / b; the span is cut at those two crossings into the stretch held + // at the bound the line starts beyond, the stretch on the ramp (a trapezoid: its length times the + // line's value at the stretch's midpoint), and the stretch held at the far bound. b = 0 is the + // constant clamp(a, lo, hi); t1 < t0 is the negative of the integral over [t1, t0]; lo above hi + // clamps to hi everywhere, as clamp does. A crossing that overflows (b tiny) lies past the span and + // cuts nothing. + public static float integral_clamp_linear(float a, float b, float lo, float hi, float t0, float t1) + { + if (b == 0.0f) + return clamp(a, lo, hi) * (t1 - t0); + + lo = min(lo, hi); + var direction = t1 < t0 ? -1.0f : 1.0f; + var start = min(t0, t1); + var stop = max(t0, t1); + var atLo = (lo - a) / b; + var atHi = (hi - a) / b; + var enter = clamp(min(atLo, atHi), start, stop); + var leave = clamp(max(atLo, atHi), start, stop); + var before = b > 0.0f ? lo : hi; + var after = b > 0.0f ? hi : lo; + var total = before * (enter - start) + (a + b * (0.5f * (enter + leave))) * (leave - enter) + after * (stop - leave); + return direction * total; + } +} diff --git a/packages/cultmath/src/CultMath/math.Span.cs.meta b/packages/cultmath/src/CultMath/math.Span.cs.meta new file mode 100644 index 00000000..31f6c409 --- /dev/null +++ b/packages/cultmath/src/CultMath/math.Span.cs.meta @@ -0,0 +1,2 @@ +fileFormatVersion: 2 +guid: 1466848dfbfc420ab4f3f67090944a67 diff --git a/packages/cultmath/tests/CultMath.Tests/GlslLowering.cs b/packages/cultmath/tests/CultMath.Tests/GlslLowering.cs index 5d119d86..5684535e 100644 --- a/packages/cultmath/tests/CultMath.Tests/GlslLowering.cs +++ b/packages/cultmath/tests/CultMath.Tests/GlslLowering.cs @@ -74,7 +74,9 @@ internal static class GlslLowering source = Regex.Replace(source, @"(?m)^[ \t]*#(ifndef|define|endif)\b.*$", string.Empty); source = Regex.Replace(source, @"(?ms)^\w+ \w+\([^)]*\b(Texture2D|SamplerState)\b[^)]*\)\s*\{.*?^\}[ \t]*$", string.Empty); - // 3. Type names: floatN, intN, uintN and boolN become vecN, ivecN, uvecN and bvecN. + // 3. Type names: floatN, intN, uintN and boolN become vecN, ivecN, uvecN and bvecN, and float3x3 + // becomes mat3 (the one matrix type the library returns: a float3x3 is built from nine scalars, + // which HLSL reads as rows and GLSL as columns, so only symmetric matrices are exported). source = LowerTypeNames(source); // 4. File-scope `static const` becomes `const`. @@ -110,7 +112,8 @@ internal static class GlslLowering source = Regex.Replace(source, @"\bfloat([234])\b", "vec$1"); source = Regex.Replace(source, @"\bint([234])\b", "ivec$1"); source = Regex.Replace(source, @"\buint([234])\b", "uvec$1"); - return Regex.Replace(source, @"\bbool([234])\b", "bvec$1"); + source = Regex.Replace(source, @"\bbool([234])\b", "bvec$1"); + return Regex.Replace(source, @"\bfloat3x3\b", "mat3"); } private static readonly Regex Cast = new(@"\((float|int|uint)\)"); diff --git a/packages/cultmath/tests/CultMath.Tests/GlslMirrorTests.cs b/packages/cultmath/tests/CultMath.Tests/GlslMirrorTests.cs index f7699379..14ea07c0 100644 --- a/packages/cultmath/tests/CultMath.Tests/GlslMirrorTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/GlslMirrorTests.cs @@ -47,7 +47,7 @@ public sealed class GlslMirrorTests // Identifier boundaries: cultmath_lerp is a CultMath function, lerp( is HLSL. private static readonly Regex HlslOnly = new( - @"\b(float|u?int|bool)[234]\b|\blerp\(|\bfrac\(|\bsaturate\(|\brsqrt\(|\basuint\(|\basfloat\(|\batan2\(|static const|\((float|int|uint)\)|Texture2D|SamplerState|#include|precision (high|medium|low)p"); + @"\b(float|u?int|bool)[234]\b|\bfloat3x3\b|\blerp\(|\bfrac\(|\bsaturate\(|\brsqrt\(|\basuint\(|\basfloat\(|\batan2\(|static const|\((float|int|uint)\)|Texture2D|SamplerState|#include|precision (high|medium|low)p"); [Fact] public void NoHlslOnlyTokensSurvive() @@ -630,6 +630,14 @@ public sealed class GlslMirrorTests new Family("float cultmath_cycles_arcsine(float)", "ulp-bounded, platform asin", r => new object[] { Uniform(r, -1.5f, 1.5f) }, Platform: Platform.Ulp), new Family("float cultmath_cycles_arccosine(float)", "ulp-bounded, platform acos", r => new object[] { Uniform(r, -1.5f, 1.5f) }, Platform: Platform.Ulp), new Family("float cultmath_cycles_arctan2(float, float)", "ulp-bounded, platform atan", r => new object[] { Uniform(r, -20.0f, 20.0f), Uniform(r, -20.0f, 20.0f) }, Platform: Platform.Ulp), + new Family("mat3 cultmath_snoise_hessian(vec3)", "ulp-bounded", r => new object[] { Point3(r) }), + new Family("mat3 cultmath_snoise_hessian(vec3)", "ulp-bounded, on the simplex tie set", r => new object[] { TiePoint3(r) }), + new Family("float cultmath_integral_clamp_linear(float, float, float, float, float, float)", "ulp-bounded", + r => + { + var lo = Uniform(r, -2.0f, 0.5f); + return new object[] { Uniform(r, -4.0f, 4.0f), r.Next(8) == 0 ? 0.0f : Uniform(r, -4.0f, 4.0f), lo, lo + Uniform(r, 0.0f, 3.0f), Uniform(r, -3.0f, 3.0f), Uniform(r, -3.0f, 3.0f) }; + }), }).ToArray(); // The fixture text, each case's arguments given by family and point index. @@ -693,8 +701,12 @@ public sealed class GlslMirrorTests _ => throw new InvalidOperationException(v.GetType().Name), })); + // A float3x3 flattens row by row; the matrices the library exports are symmetric, so the order is + // the same read as GLSL's columns. private static IEnumerable Scalars(object value) => value is float or int or uint ? new[] { value } + : value is float3x3 matrix + ? new object[] { matrix._m00, matrix._m01, matrix._m02, matrix._m10, matrix._m11, matrix._m12, matrix._m20, matrix._m21, matrix._m22 } : value.GetType().GetFields().Where(f => !f.IsStatic).SelectMany(f => Scalars(f.GetValue(value)!)); // ---- A small unified diff for the failure message ---- diff --git a/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs b/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs index 989ddf26..4ca598ad 100644 --- a/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs @@ -28,6 +28,8 @@ public sealed class HlslSourceCompatibilityTests Assert.NotNull(ShaderFunction(shader, "cultmath_smin_grad", typeof(float4), typeof(float4), typeof(float))); Assert.NotNull(ShaderFunction(shader, "cultmath_cellular", typeof(float3))); Assert.NotNull(ShaderFunction(shader, "cultmath_snoise_grad", typeof(float3))); + Assert.NotNull(ShaderFunction(shader, "cultmath_snoise_hessian", typeof(float3))); + Assert.NotNull(ShaderFunction(shader, "cultmath_integral_clamp_linear", typeof(float), typeof(float), typeof(float), typeof(float), typeof(float), typeof(float))); 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))); @@ -137,7 +139,7 @@ public sealed class HlslSourceCompatibilityTests } Assert.True(mismatches.Count == 0, $"{mismatches.Count} mismatches:{Environment.NewLine}{string.Join(Environment.NewLine, mismatches.Take(12))}"); - Assert.Equal(118, compared.Count); + Assert.Equal(120, compared.Count); } /// @@ -377,6 +379,9 @@ public sealed class HlslSourceCompatibilityTests /// private static Leaf[] Values(object value, string path = "") => Scalars.Contains(value.GetType()) ? new[] { new Leaf(path, value.GetType(), value) } + : value is float3x3 matrix + ? new[] { matrix._m00, matrix._m01, matrix._m02, matrix._m10, matrix._m11, matrix._m12, matrix._m20, matrix._m21, matrix._m22 } + .Select((entry, k) => new Leaf((path.Length == 0 ? "" : path + ".") + "_m" + k / 3 + k % 3, typeof(float), entry)).ToArray() : value.GetType().GetFields().Where(f => !f.IsStatic) .SelectMany(f => Values(f.GetValue(value)!, path.Length == 0 ? f.Name : path + "." + f.Name)) .ToArray(); diff --git a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs index 78a3b795..06b0ea3b 100644 --- a/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/NoiseBoundTests.cs @@ -489,25 +489,13 @@ public sealed class NoiseBoundTests // the quotient is not skewed by the rounding of p + h far from the origin. private static double Step(float x, float h, out double back) { back = x - (double)(x - h); return (x + h) - (double)x; } - // The spectral norm (largest absolute eigenvalue) of the symmetrised Hessian of snoise(float3) at p, - // from central differences of snoise_grad at h = 1e-3. + // The spectral norm (largest absolute eigenvalue) of snoise's Hessian at p: the library's own + // snoise_hessian, symmetric bit for bit, then the closed-form eigenvalues of a symmetric 3x3. private static double HessianNorm(float3 p) { - const float h = 1.0e-3f; - var m = new double[3, 3]; - for (var a = 0; a < 3; a++) - { - var e = new float3(a == 0 ? h : 0.0f, a == 1 ? h : 0.0f, a == 2 ? h : 0.0f); - var forward = Step(a == 0 ? p.x : a == 1 ? p.y : p.z, h, out var back); - var gp = snoise_grad(p + e); - var gm = snoise_grad(p - e); - m[0, a] = (gp.x - (double)gm.x) / (forward + back); - m[1, a] = (gp.y - (double)gm.y) / (forward + back); - m[2, a] = (gp.z - (double)gm.z) / (forward + back); - } - - var (a00, a11, a22) = (m[0, 0], m[1, 1], m[2, 2]); - var (a01, a02, a12) = ((m[0, 1] + m[1, 0]) / 2, (m[0, 2] + m[2, 0]) / 2, (m[1, 2] + m[2, 1]) / 2); + var m = snoise_hessian(p); + var (a00, a11, a22) = ((double)m._m00, (double)m._m11, (double)m._m22); + var (a01, a02, a12) = ((double)m._m01, (double)m._m02, (double)m._m12); var p1 = a01 * a01 + a02 * a02 + a12 * a12; if (p1 == 0.0) return Math.Max(Math.Abs(a00), Math.Max(Math.Abs(a11), Math.Abs(a22))); @@ -635,6 +623,180 @@ public sealed class NoiseBoundTests Assert.True(Math.Abs(ratio - 1.0) <= 0.005, $"SNOISE_HESSIAN is {ratio:R} times 1.10 x the re-measured {refined:R}; re-run MeasureHessian and re-pin"); } + // ---- The third-derivative constant (math.Span.cs) ---- + + // The symmetric third-derivative tensor of snoise at p, entry (i * 3 + j) * 3 + k: central differences + // of snoise_hessian along each axis at h = 1e-3 over the exact float32 steps, symmetrised over the six + // index orders (the Hessian is symmetric, so only the derivative index can differ). + private static double[] ThirdTensor(float3 p) + { + const float h = 1.0e-3f; + var d = new double[27]; + for (var k = 0; k < 3; k++) + { + var e = new float3(k == 0 ? h : 0.0f, k == 1 ? h : 0.0f, k == 2 ? h : 0.0f); + var forward = Step(k == 0 ? p.x : k == 1 ? p.y : p.z, h, out var back); + var plus = snoise_hessian(p + e); + var minus = snoise_hessian(p - e); + for (var i = 0; i < 3; i++) + for (var j = 0; j < 3; j++) + d[(i * 3 + j) * 3 + k] = (plus[i][j] - (double)minus[i][j]) / (forward + back); + } + + var t = new double[27]; + for (var i = 0; i < 3; i++) + for (var j = 0; j < 3; j++) + for (var k = 0; k < 3; k++) + t[(i * 3 + j) * 3 + k] = (d[(i * 3 + j) * 3 + k] + d[(i * 3 + k) * 3 + j] + d[(j * 3 + i) * 3 + k] + + d[(j * 3 + k) * 3 + i] + d[(k * 3 + i) * 3 + j] + d[(k * 3 + j) * 3 + i]) / 6.0; + return t; + } + + // T(u, u, u) for a symmetric 3-tensor. + private static double Cubic(double[] t, double[] u) + { + var sum = 0.0; + for (var i = 0; i < 3; i++) + for (var j = 0; j < 3; j++) + for (var k = 0; k < 3; k++) + sum += t[(i * 3 + j) * 3 + k] * u[i] * u[j] * u[k]; + return sum; + } + + // The 26 unit directions of the axes, the edge midpoints and the cube corners. + private static readonly double[][] TensorStarts = Enumerable.Range(0, 27).Select(c => new[] { c % 3 - 1.0, c / 3 % 3 - 1.0, c / 9 - 1.0 }) + .Where(v => v[0] != 0.0 || v[1] != 0.0 || v[2] != 0.0) + .Select(v => { var length = Math.Sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); return new[] { v[0] / length, v[1] / length, v[2] / length }; }).ToArray(); + + // The operator norm of a symmetric 3-tensor, max over unit u of |T(u, u, u)| (for a symmetric tensor + // that is the injective norm). From each of 26 starts, the shifted power ascent u <- normalise(T(u, u, .) + // +- alpha u), with alpha twice the largest absolute row sum, which makes the iteration monotone in + // the cubic form (Kolda and Mayo's SS-HOPM), the sign chosen as the start's. A local maximum, so a + // lower bound; MeasureThird checks it against a dense direction sweep at its witness. + private static double TensorNorm(double[] t) + { + var alpha = 0.0; + for (var i = 0; i < 3; i++) + { + var row = 0.0; + for (var jk = 0; jk < 9; jk++) + row += Math.Abs(t[i * 9 + jk]); + alpha = Math.Max(alpha, 2.0 * row); + } + + var best = 0.0; + var g = new double[3]; + foreach (var start in TensorStarts) + { + var u = (double[])start.Clone(); + var sign = Cubic(t, u) < 0.0 ? -1.0 : 1.0; + for (var iteration = 0; iteration < 60; iteration++) + { + for (var i = 0; i < 3; i++) + { + var sum = 0.0; + for (var j = 0; j < 3; j++) + for (var k = 0; k < 3; k++) + sum += t[(i * 3 + j) * 3 + k] * u[j] * u[k]; + g[i] = sign * sum + alpha * u[i]; + } + + var length = Math.Sqrt(g[0] * g[0] + g[1] * g[1] + g[2] * g[2]); + if (length == 0.0) + break; + var moved = 0.0; + for (var i = 0; i < 3; i++) + { + var next = g[i] / length; + moved += (next - u[i]) * (next - u[i]); + u[i] = next; + } + + if (moved < 1.0e-20) + break; + } + + best = Math.Max(best, Math.Abs(Cubic(t, u))); + } + + return best; + } + + // The operator norm of snoise's third-derivative tensor at p. + private static double ThirdNorm(float3 p) => TensorNorm(ThirdTensor(p)); + + [Fact] + public void TensorNormFindsKnownMaxima() + { + // T = e_x (x) e_x (x) e_x: the cubic form is u_x^3, norm 1. + var single = new double[27]; + single[0] = 1.0; + Assert.InRange(TensorNorm(single), 1.0 - 1.0e-9, 1.0 + 1.0e-9); + + // The six permutations of (0, 1, 2) at 1: the cubic form is 6 u_x u_y u_z, largest on the cube + // diagonal at 6 / (3 sqrt 3) = 2 / sqrt 3. + var permutations = new double[27]; + foreach (var (i, j, k) in new[] { (0, 1, 2), (0, 2, 1), (1, 0, 2), (1, 2, 0), (2, 0, 1), (2, 1, 0) }) + permutations[(i * 3 + j) * 3 + k] = 1.0; + Assert.InRange(TensorNorm(permutations), 2.0 / Math.Sqrt(3.0) - 1.0e-9, 2.0 / Math.Sqrt(3.0) + 1.0e-9); + + // A negative extremum counts by its absolute value. + var negative = new double[27]; + negative[13] = -3.0; + Assert.InRange(TensorNorm(negative), 3.0 - 1.0e-9, 3.0 + 1.0e-9); + } + + /// + /// The provenance of SNOISE_THIRD: the largest operator norm of snoise's third-derivative tensor + /// (ThirdNorm) over 1e6 seeded points of [-256, 256]^3, the 1e4 largest refined by ascent over the point, + /// times 1.10, on the tied kernel. Prints the start whose ascent reached the maximum, which + /// ThirdConstantPinsSampled climbs from again, and a sweep of 20,000 directions at the point the + /// ascent ends on (the tensor norm there against the largest cubic form the sweep finds, so the + /// tightness of the inner maximisation is on the record). 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 MeasureThird() + { + var (sampled, refined, witness) = MeasureMax(ThirdNorm, false, 0x4E57, 1_000_000, 10_000); + output.WriteLine($"IV-REPORT third: max |D3| sampled {sampled:R}, refined {refined:R}, x 1.10 = {(float)(refined * 1.10):R}, witness ({witness.x:R}, {witness.y:R}, {witness.z:R})"); + var (best, at) = AscendPoint(ThirdNorm, witness); + var tensor = ThirdTensor(at); + var sweep = 0.0; + const int directions = 20_000; + for (var k = 0; k < directions; k++) + { + var z = 1.0 - 2.0 * (k + 0.5) / directions; + var r = Math.Sqrt(1.0 - z * z); + var phi = k * Math.PI * (3.0 - Math.Sqrt(5.0)); + sweep = Math.Max(sweep, Math.Abs(Cubic(tensor, new[] { r * Math.Cos(phi), r * Math.Sin(phi), z }))); + } + + output.WriteLine($"IV-REPORT third witness ascent ends at ({at.x:R}, {at.y:R}, {at.z:R}): tensor norm {best:R}, dense sweep of {directions} directions {sweep:R} ({sweep / best:R} of it)"); + Assert.True(sweep <= best * (1.0 + 1.0e-6), $"a direction sweep finds {sweep:R} over the ascent's {best:R}: the inner maximisation misses a maximum"); + Assert.True(refined <= SNOISE_THIRD, $"measured {refined:R} exceeds SNOISE_THIRD {SNOISE_THIRD:R}"); + } + + // The start MeasureThird printed as its witness. + private static readonly float3 ThirdWitness = new(0.0f, 0.0f, 0.0f); + + /// + /// Pins SNOISE_THIRD: no third-derivative norm over 1e5 seeded points exceeds it, and the ascent from + /// MeasureThird's witness reaches a maximum that, times 1.10, is the constant to within half a + /// percent, so a constant moved by 1% fails here. + /// + [Fact] + public void ThirdConstantPinsSampled() + { + var (sampled, _, _) = MeasureMax(ThirdNorm, false, 0x4E58, 100_000, 0); + Assert.True(sampled <= SNOISE_THIRD, $"the third-derivative norm reaches {sampled:R}, above SNOISE_THIRD {SNOISE_THIRD:R}"); + var refined = Ascend(ThirdNorm, ThirdWitness); + Assert.True(refined <= SNOISE_THIRD, $"the third-derivative norm reaches {refined:R}, above SNOISE_THIRD {SNOISE_THIRD:R}"); + var ratio = SNOISE_THIRD / (1.10 * refined); + Assert.True(Math.Abs(ratio - 1.0) <= 0.005, $"SNOISE_THIRD is {ratio:R} times 1.10 x the re-measured {refined:R}; re-run MeasureThird and re-pin"); + } + /// /// SNOISE_HESSIAN where the curvature is realised. At the point the ascent from MeasureHessian's witness /// reaches, the absolute second difference of snoise along 4,000 seeded directions u at step 0.05 has a diff --git a/packages/cultmath/tests/CultMath.Tests/SpanPrimitiveTests.cs b/packages/cultmath/tests/CultMath.Tests/SpanPrimitiveTests.cs new file mode 100644 index 00000000..82da11b7 --- /dev/null +++ b/packages/cultmath/tests/CultMath.Tests/SpanPrimitiveTests.cs @@ -0,0 +1,329 @@ +using CultMath; +using Xunit; +using static CultMath.math; + +namespace CultMath.Tests; + +/// +/// The span primitives (math.Span.cs): snoise_hessian against the finite differences of snoise_grad and of +/// snoise, with each tolerance measured by Richardson's two step sizes rather than assumed; its symmetry and +/// its continuity across the simplex diagonal; and integral_clamp_linear against adaptive quadrature in +/// double. SNOISE_THIRD is measured and pinned in NoiseBoundTests (MeasureThird, ThirdConstantPinsSampled). +/// +/// Hand mutations these must kill: any term of the per-corner second derivative dropped or mis-signed +/// (48 m^2 s x_i x_j, 8 m^3 p_j x_i, 8 m^3 p_i x_j, 8 m^3 s delta_ij) fails +/// SnoiseHessianMatchesCentralDifferences; the 0.3.0 cyclic tie order fails +/// SnoiseHessianContinuousOnSimplexDiagonal; a crossing of integral_clamp_linear moved, a clamp piece +/// dropped or the reversed-span sign lost fails IntegralClampLinearMatchesQuadrature. +/// +public sealed class SpanPrimitiveTests +{ + private readonly ITestOutputHelper output; + + public SpanPrimitiveTests(ITestOutputHelper output) => this.output = output; + + private const double Unit = 1.0 / 16777216.0; + + // A point of [-50, 50]^3 on the 2^-10 lattice, so that p +- h is exact in float32 for h a power of two + // up to 2^-4 and the finite-difference quotients carry no step rounding. + private static float3 LatticePoint(System.Random random) + { + float Axis() => MathF.Round((random.NextSingle() * 100.0f - 50.0f) * 1024.0f) / 1024.0f; + return new float3(Axis(), Axis(), Axis()); + } + + // The simplex tie set of NoiseGradTests.SnoiseIsContinuousOnTheSimplexDiagonal: the line (t, t, t) and + // the lattice diagonals v + s (1, 1, 1), v = L - sum(L) / 6. + private static List TiePoints() + { + const int vertices = 160; + const int perVertex = 8; + var random = new System.Random(0x71E5); + var points = new List(); + for (var k = 0; k < 400; k++) + { + var t = (float)(random.NextDouble() * 128.0 - 64.0); + points.Add(new float3(t, t, t)); + } + + points.Add(new float3(0.0f, 0.0f, 0.0f)); + for (var k = 0; k < vertices; k++) + { + var lx = random.Next(-64, 65); + var ly = random.Next(-64, 65); + var lz = random.Next(-64, 65); + var shift = (lx + ly + lz) / 6.0f; + for (var j = 0; j < perVertex; j++) + { + var s = j == 0 ? 0.0f : random.NextSingle(); + points.Add(new float3(lx - shift + s, ly - shift + s, lz - shift + s)); + } + } + + return points; + } + + private static float3 Axis(int a, float h) => new(a == 0 ? h : 0.0f, a == 1 ? h : 0.0f, a == 2 ? h : 0.0f); + + // Column a of the Hessian by the central difference of snoise_grad at step h, in double. + private static double[] GradientDifference(float3 p, int a, float h) + { + var plus = snoise_grad(p + Axis(a, h)); + var minus = snoise_grad(p - Axis(a, h)); + return new[] { (plus.x - (double)minus.x) / (2.0 * h), (plus.y - (double)minus.y) / (2.0 * h), (plus.z - (double)minus.z) / (2.0 * h) }; + } + + // Worst ratio of |H - D(h/2)| to its tolerance, the number of entries compared and the largest |H| met. + private (double Worst, int Compared, double LargestEntry, double LargestTolerance) CompareToGradientDifferences(IEnumerable points) + { + const float h = 1.0f / 256.0f; + double worst = 0.0, largest = 0.0, largestTolerance = 0.0; + var compared = 0; + foreach (var p in points) + { + var hessian = snoise_hessian(p); + var gradient = snoise_grad(p); + var scale = Math.Max(Math.Max(Math.Abs(gradient.x), Math.Abs(gradient.y)), Math.Max(Math.Abs(gradient.z), 1.0)); + for (var a = 0; a < 3; a++) + { + var coarse = GradientDifference(p, a, h); + var fine = GradientDifference(p, a, h / 2.0f); + for (var i = 0; i < 3; i++) + { + // The truncation of the finer quotient is a third of the gap between the two, as the + // error is O(h^2); the roundoff is a few ulp of the gradient over the step. Both are + // measured per entry and tripled, so the tolerance follows the point's own roughness. + var tolerance = 3.0 * (Math.Abs(coarse[i] - fine[i]) + 64.0 * Unit * scale / (h / 2.0)); + var error = Math.Abs(hessian[i][a] - fine[i]); + worst = Math.Max(worst, error / tolerance); + largest = Math.Max(largest, Math.Abs(hessian[i][a])); + largestTolerance = Math.Max(largestTolerance, tolerance); + Assert.True(error <= tolerance, $"snoise_hessian({p})[{i}][{a}] = {hessian[i][a]:R} differs from the central difference of snoise_grad {fine[i]:R} by {error:R}, over the measured tolerance {tolerance:R}"); + compared++; + } + } + } + + return (worst, compared, largest, largestTolerance); + } + + [Fact] + public void SnoiseHessianMatchesCentralDifferences() + { + var random = new System.Random(0x5A11); + var points = Enumerable.Range(0, 10_000).Select(_ => LatticePoint(random)).ToList(); + var (worst, compared, largest, largestTolerance) = CompareToGradientDifferences(points); + output.WriteLine($"snoise_hessian vs central differences of snoise_grad over {points.Count} seeded points: {compared} entries, worst error {worst:F3} of its Richardson tolerance, largest |H| {largest:R}, largest tolerance {largestTolerance:R}"); + + // The comparison has teeth: curvature of the size the constant claims is present, and the tolerance + // is far under it, so a dropped term of the second derivative cannot hide in it. + Assert.True(largest > 10.0, $"the sample never meets curvature above 10 (largest {largest:R})"); + Assert.True(largestTolerance < 0.5, $"the tolerance reaches {largestTolerance:R}, too wide to see a dropped term"); + + var (tieWorst, tieCompared, tieLargest, tieTolerance) = CompareToGradientDifferences(TiePoints()); + output.WriteLine($"snoise_hessian vs central differences of snoise_grad on the simplex tie set: {tieCompared} entries, worst error {tieWorst:F3} of its Richardson tolerance, largest |H| {tieLargest:R}, largest tolerance {tieTolerance:R}"); + Assert.True(tieLargest > 1.0, $"the tie set never meets curvature above 1 (largest {tieLargest:R})"); + } + + // The second difference of snoise in double, h apart, for the entry (i, j): the three-point stencil + // for i = j, the four-point mixed stencil otherwise. + private static double SecondDifference(float3 p, int i, int j, float h) + { + if (i == j) + return (snoise(p + Axis(i, h)) - 2.0 * (double)snoise(p) + snoise(p - Axis(i, h))) / ((double)h * h); + var pp = snoise(p + Axis(i, h) + Axis(j, h)); + var pm = snoise(p + Axis(i, h) - Axis(j, h)); + var mp = snoise(p - Axis(i, h) + Axis(j, h)); + var mm = snoise(p - Axis(i, h) - Axis(j, h)); + return ((double)pp - pm - mp + mm) / (4.0 * h * h); + } + + [Fact] + public void SnoiseHessianMatchesSecondDifferencesOfSnoise() + { + var random = new System.Random(0x5A12); + double worst = 0.0, largest = 0.0; + var compared = 0; + for (var k = 0; k < 4000; k++) + { + var p = LatticePoint(random); + var hessian = snoise_hessian(p); + for (var i = 0; i < 3; i++) + for (var j = i; j < 3; j++) + { + // Richardson on the O(h^2) second difference: the extrapolation (4 S(h/2) - S(h)) / 3 is + // O(h^4); its own error is bounded by the gap between the two quotients, and the roundoff + // of snoise (a few ulp) over h^2 is added. + const float h = 1.0f / 16.0f; + var coarse = SecondDifference(p, i, j, h); + var fine = SecondDifference(p, i, j, h / 2.0f); + var extrapolated = (4.0 * fine - coarse) / 3.0; + var tolerance = 2.0 * Math.Abs(coarse - fine) + 64.0 * Unit / ((h / 2.0) * (h / 2.0)); + var error = Math.Abs(hessian[i][j] - extrapolated); + worst = Math.Max(worst, error / tolerance); + largest = Math.Max(largest, Math.Abs(hessian[i][j])); + compared++; + Assert.True(error <= tolerance, $"snoise_hessian({p})[{i}][{j}] = {hessian[i][j]:R} differs from the extrapolated second difference of snoise {extrapolated:R} by {error:R}, over the measured tolerance {tolerance:R}"); + } + } + + output.WriteLine($"snoise_hessian vs Richardson second differences of snoise over 4000 seeded points: {compared} entries, worst error {worst:F3} of its tolerance, largest |H| {largest:R}"); + Assert.True(largest > 10.0, $"the sample never meets curvature above 10 (largest {largest:R})"); + } + + [Fact] + public void SnoiseHessianIsSymmetric() + { + var random = new System.Random(0x5A13); + var points = Enumerable.Range(0, 20_000).Select(_ => new float3(random.NextSingle() * 512.0f - 256.0f, random.NextSingle() * 512.0f - 256.0f, random.NextSingle() * 512.0f - 256.0f)) + .Concat(TiePoints()).Append(new float3(0.0f, 0.0f, 0.0f)).ToList(); + foreach (var p in points) + { + var h = snoise_hessian(p); + Assert.True(BitConverter.SingleToInt32Bits(h._m01) == BitConverter.SingleToInt32Bits(h._m10) + && BitConverter.SingleToInt32Bits(h._m02) == BitConverter.SingleToInt32Bits(h._m20) + && BitConverter.SingleToInt32Bits(h._m12) == BitConverter.SingleToInt32Bits(h._m21), $"snoise_hessian({p}) is not its own transpose: {h}"); + } + } + + [Fact] + public void SnoiseHessianContinuousOnSimplexDiagonal() + { + // Within 1 or 2 ulp of a tie point along each axis the Hessian moves by at most the third derivative + // bound times the distance: SNOISE_THIRD x 2 sqrt 3 x 1e-4 covers a 1e-4 ball, and the points are + // far closer than that. The 0.3.0 cyclic tie order is a jump of order one here. + var bound = SNOISE_THIRD * 2.0 * Math.Sqrt(3.0) * 1.0e-4; + var worst = 0.0; + foreach (var p in TiePoints()) + { + var centre = snoise_hessian(p); + for (var axis = 0; axis < 3; axis++) + { + for (var direction = -1; direction <= 1; direction += 2) + { + var q = p; + for (var ulp = 1; ulp <= 2; ulp++) + { + var component = axis == 0 ? q.x : axis == 1 ? q.y : q.z; + component = direction < 0 ? MathF.BitDecrement(component) : MathF.BitIncrement(component); + q = axis == 0 ? new float3(component, q.y, q.z) : axis == 1 ? new float3(q.x, component, q.z) : new float3(q.x, q.y, component); + var moved = snoise_hessian(q); + for (var i = 0; i < 3; i++) + for (var j = 0; j < 3; j++) + { + var jump = Math.Abs(moved[i][j] - (double)centre[i][j]); + worst = Math.Max(worst, jump); + Assert.True(jump <= bound, $"snoise_hessian jumps {jump:R} in entry [{i}][{j}] within {ulp} ulp of {p}, over {bound:R}"); + } + } + } + } + } + + output.WriteLine($"largest Hessian jump within 2 ulp of the simplex tie set: {worst:R} against the bound {bound:R}"); + } + + // Adaptive Simpson of f over [a, b] in double, to an absolute tolerance, with a depth cap. + private static double Quadrature(Func f, double a, double b, double tolerance) + { + double Simpson(double fa, double fm, double fb, double width) => width / 6.0 * (fa + 4.0 * fm + fb); + + double Recurse(double x0, double x1, double f0, double fm, double f1, double whole, double eps, int depth) + { + var mid = 0.5 * (x0 + x1); + var lm = 0.5 * (x0 + mid); + var rm = 0.5 * (mid + x1); + var flm = f(lm); + var frm = f(rm); + var left = Simpson(f0, flm, fm, mid - x0); + var right = Simpson(fm, frm, f1, x1 - mid); + if (depth <= 0 || Math.Abs(left + right - whole) <= 15.0 * eps) + return left + right + (left + right - whole) / 15.0; + return Recurse(x0, mid, f0, flm, fm, left, eps / 2.0, depth - 1) + Recurse(mid, x1, fm, frm, f1, right, eps / 2.0, depth - 1); + } + + if (a == b) + return 0.0; + var fa = f(a); + var fb = f(b); + var fmid = f(0.5 * (a + b)); + return Recurse(a, b, fa, fmid, fb, Simpson(fa, fmid, fb, b - a), tolerance, 40); + } + + private static double ClampedLine(double a, double b, double lo, double hi, double t) => Math.Min(Math.Max(a + b * t, lo), hi); + + // The reference integral in double: the adaptive quadrature of the clamped line, sign for a reversed span. + private static double ReferenceIntegral(float a, float b, float lo, float hi, float t0, float t1) + { + double low = Math.Min(lo, hi); + var integral = t0 <= t1 + ? Quadrature(t => ClampedLine(a, b, low, hi, t), t0, t1, 1.0e-12) + : -Quadrature(t => ClampedLine(a, b, low, hi, t), t1, t0, 1.0e-12); + return integral; + } + + [Fact] + public void IntegralClampLinearMatchesQuadrature() + { + var random = new System.Random(0x1A7E); + float Uniform(float lo, float hi) => lo + (hi - lo) * random.NextSingle(); + double worstRatio = 0.0; + int zeroSlopes = 0, reversed = 0, both = 0, neither = 0, ramps = 0; + for (var k = 0; k < 100_000; k++) + { + var a = Uniform(-4.0f, 4.0f); + // One in eight slopes is zero, one in sixteen is tiny (a crossing far outside any span). + var b = (k % 8) switch { 0 => 0.0f, 1 => 1.0e-30f * (k % 16 == 1 ? 1.0f : -1.0f), _ => Uniform(-4.0f, 4.0f) }; + var lo = Uniform(-2.0f, 0.5f); + var hi = lo + Uniform(0.0f, 3.0f); + var t0 = Uniform(-3.0f, 3.0f); + var t1 = Uniform(-3.0f, 3.0f); + var got = integral_clamp_linear(a, b, lo, hi, t0, t1); + var want = ReferenceIntegral(a, b, lo, hi, t0, t1); + var span = Math.Abs((double)t1 - t0); + var error = Math.Abs(got - want); + worstRatio = Math.Max(worstRatio, error / Math.Max(1.0e-5 * Math.Abs(want), 1.0e-6 * span)); + Assert.True(error <= 1.0e-5 * Math.Abs(want) || error <= 1.0e-6 * span, $"integral_clamp_linear({a:R}, {b:R}, {lo:R}, {hi:R}, {t0:R}, {t1:R}) = {got:R}, quadrature {want:R}"); + if (b == 0.0f) + zeroSlopes++; + if (t1 < t0) + reversed++; + if (b != 0.0f) + { + double atLo = (lo - (double)a) / b, atHi = (hi - (double)a) / b; + double from = Math.Min(t0, t1), to = Math.Max(t0, t1); + var crossesLo = atLo > from && atLo < to; + var crossesHi = atHi > from && atHi < to; + if (crossesLo && crossesHi) + both++; + else if (!crossesLo && !crossesHi) + neither++; + if (hi > lo && (crossesLo || crossesHi)) + ramps++; + } + } + + output.WriteLine($"integral_clamp_linear vs adaptive quadrature over 100000 cases: {zeroSlopes} flat, {reversed} reversed, {both} crossing both bounds, {neither} crossing neither, {ramps} on a ramp; worst error {worstRatio:F3} of its tolerance"); + Assert.True(zeroSlopes > 10_000 && reversed > 40_000 && both > 2_000 && neither > 2_000, $"the sample meets a branch too rarely: {zeroSlopes} flat, {reversed} reversed, {both} both, {neither} neither"); + } + + [Fact] + public void IntegralClampLinearHasExactValuesOnTheBinaryFloats() + { + // The ramp t from 0 to 1 clamped to [0.25, 0.75] over [0, 1]: 0.25 x 0.25 below, 0.5 x 0.5 on the ramp, 0.75 x 0.25 above. + Assert.Equal(0.0625f + 0.25f + 0.1875f, integral_clamp_linear(0.0f, 1.0f, 0.25f, 0.75f, 0.0f, 1.0f)); + // A reversed span is the negative of the span. + Assert.Equal(-0.5f, integral_clamp_linear(0.0f, 1.0f, 0.25f, 0.75f, 1.0f, 0.0f)); + Assert.Equal(integral_clamp_linear(0.3f, -0.7f, -0.5f, 0.9f, -2.0f, 1.5f), -integral_clamp_linear(0.3f, -0.7f, -0.5f, 0.9f, 1.5f, -2.0f)); + // b = 0 is the constant, clamped. + Assert.Equal(2.0f * 3.0f, integral_clamp_linear(5.0f, 0.0f, 0.0f, 2.0f, 1.0f, 4.0f)); + Assert.Equal(-1.0f * 3.0f, integral_clamp_linear(-5.0f, 0.0f, -1.0f, 2.0f, 1.0f, 4.0f)); + // An empty span, and a span that never leaves the ramp (the trapezoid: 3 x mean of 1 and 4). + Assert.Equal(0.0f, integral_clamp_linear(0.0f, 1.0f, -9.0f, 9.0f, 2.0f, 2.0f)); + Assert.Equal(7.5f, integral_clamp_linear(0.0f, 1.0f, -9.0f, 9.0f, 1.0f, 4.0f)); + // A falling line starts at the upper bound: 1 for t in [0, 1], the ramp 1 to 0 over [1, 2], 0 after. + Assert.Equal(1.0f + 0.5f, integral_clamp_linear(2.0f, -1.0f, 0.0f, 1.0f, 0.0f, 3.0f)); + // lo above hi clamps to hi everywhere. + Assert.Equal(2.0f * 3.0f, integral_clamp_linear(0.0f, 1.0f, 5.0f, 2.0f, 0.0f, 3.0f)); + } +} diff --git a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Affine.hlsl b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Affine.hlsl index f165148b..8727081e 100644 --- a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Affine.hlsl +++ b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Affine.hlsl @@ -12,6 +12,10 @@ // math.Affine.cs (SNOISE_HESSIAN). static const float CULTMATH_SNOISE_HESSIAN = 56.050385; +// An upper bound on the operator norm of the third-derivative tensor of cultmath_snoise(float3), the +// next order after CULTMATH_SNOISE_HESSIAN; provenance in math.Span.cs (SNOISE_THIRD). +static const float CULTMATH_SNOISE_THIRD = 1.0; + float3 cultmath_af_point(float value) { return float3(value, 0.0, 0.0); } float3 cultmath_af_symbol(float x0, float a) { return float3(x0, a, 0.0); } diff --git a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl index c29fd603..18cd1fc9 100644 --- a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl +++ b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl @@ -228,6 +228,96 @@ float4 cultmath_snoise_grad(float3 value) return float4(grad, value2); } +// Analytic Hessian of cultmath_snoise(float3); see the C# math.snoise_hessian comment for the +// derivation. Every local up through p0..p3 and x0..x3 is cultmath_snoise_grad's own. A float3x3 built +// from nine scalars reads as rows here and as columns in GLSL's mat3: the matrix is symmetric. +float3x3 cultmath_snoise_hessian(float3 value) +{ + float2 c = float2(1.0 / 6.0, 1.0 / 3.0); + float3 i = floor(value + dot(value, c.yyy)); + float3 x0 = value - i + dot(i, c.xxx); + float3 g = float3(step(x0.y, x0.x), step(x0.z, x0.y), 1.0 - step(x0.z, x0.x)); + float3 l = 1.0 - g; + float3 i1 = min(g.xyz, l.zxy); + float3 i2 = max(g.xyz, l.zxy); + float3 x1 = x0 - i1 + c.xxx; + float3 x2 = x0 - i2 + c.yyy; + float3 x3 = x0 - 0.5; + + i = cultmath_snoise_mod289(i); + float4 p = cultmath_snoise_permute(cultmath_snoise_permute(cultmath_snoise_permute( + i.z + float4(0.0, i1.z, i2.z, 1.0)) + i.y + float4(0.0, i1.y, i2.y, 1.0)) + i.x + float4(0.0, i1.x, i2.x, 1.0)); + + const float n_ = 0.142857142857; + float3 ns = float3(2.0 * n_, 0.5 * n_ - 1.0, n_); + float4 j = p - 49.0 * floor(p * ns.z * ns.z); + float4 x_ = floor(j * ns.z); + float4 y_ = floor(j - 7.0 * x_); + float4 x = x_ * ns.x + ns.yyyy; + float4 y = y_ * ns.x + ns.yyyy; + float4 h = 1.0 - abs(x) - abs(y); + + float4 b0 = float4(x.xy, y.xy); + float4 b1 = float4(x.zw, y.zw); + float4 s0 = floor(b0) * 2.0 + 1.0; + float4 s1 = floor(b1) * 2.0 + 1.0; + float4 sh = -step(h, float4(0.0, 0.0, 0.0, 0.0)); + float4 a0 = b0.xzyw + s0.xzyw * sh.xxyy; + float4 a1 = b1.xzyw + s1.xzyw * sh.zzww; + + float3 p0 = float3(a0.xy, h.x); + float3 p1 = float3(a0.zw, h.y); + float3 p2 = float3(a1.xy, h.z); + float3 p3 = float3(a1.zw, h.w); + float4 norm = 1.79284291400159 - 0.85373472095314 * float4(dot(p0, p0), dot(p1, p1), dot(p2, p2), dot(p3, p3)); + p0 *= norm.x; + p1 *= norm.y; + p2 *= norm.z; + p3 *= norm.w; + + float4 m0 = max(0.5 - float4(dot(x0, x0), dot(x1, x1), dot(x2, x2), dot(x3, x3)), 0.0); + float4 m2 = m0 * m0; + float4 m3 = m2 * m0; + float4 px = float4(dot(p0, x0), dot(p1, x1), dot(p2, x2), dot(p3, x3)); + float4 wide = 48.0 * m2 * px; + float4 narrow = 8.0 * m3; + + float4 cxx = float4(x0.x, x1.x, x2.x, x3.x); + float4 cxy = float4(x0.y, x1.y, x2.y, x3.y); + float4 cxz = float4(x0.z, x1.z, x2.z, x3.z); + float4 cpx = float4(p0.x, p1.x, p2.x, p3.x); + float4 cpy = float4(p0.y, p1.y, p2.y, p3.y); + float4 cpz = float4(p0.z, p1.z, p2.z, p3.z); + + float hxx = 105.0 * (dot(wide * cxx, cxx) - dot(narrow, 2.0 * cpx * cxx + px)); + float hyy = 105.0 * (dot(wide * cxy, cxy) - dot(narrow, 2.0 * cpy * cxy + px)); + float hzz = 105.0 * (dot(wide * cxz, cxz) - dot(narrow, 2.0 * cpz * cxz + px)); + float hxy = 105.0 * (dot(wide * cxx, cxy) - dot(narrow, cpx * cxy + cpy * cxx)); + float hxz = 105.0 * (dot(wide * cxx, cxz) - dot(narrow, cpx * cxz + cpz * cxx)); + float hyz = 105.0 * (dot(wide * cxy, cxz) - dot(narrow, cpy * cxz + cpz * cxy)); + return float3x3(hxx, hxy, hxz, hxy, hyy, hyz, hxz, hyz, hzz); +} + +// The integral over [t0, t1] of clamp(a + b t, lo, hi); see the C# math.integral_clamp_linear comment. +float cultmath_integral_clamp_linear(float a, float b, float lo, float hi, float t0, float t1) +{ + if (b == 0.0) + return clamp(a, lo, hi) * (t1 - t0); + + lo = min(lo, hi); + float direction = t1 < t0 ? -1.0 : 1.0; + float start = min(t0, t1); + float stop = max(t0, t1); + float atLo = (lo - a) / b; + float atHi = (hi - a) / b; + float enter = clamp(min(atLo, atHi), start, stop); + float leave = clamp(max(atLo, atHi), start, stop); + float before = b > 0.0 ? lo : hi; + float after = b > 0.0 ? hi : lo; + float total = before * (enter - start) + (a + b * (0.5 * (enter + leave))) * (leave - enter) + after * (stop - leave); + return direction * total; +} + // Octave sum of cultmath_snoise_grad (invariant 8); see the C# math.fbm_grad comment for the // per-octave frequency/amplitude derivation and why octaves is clamped. float4 cultmath_fbm_grad(float3 p, int octaves, float lacunarity, float gain)