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

Merge hands/cultmath-snoise-hessian

Browse source
This commit is contained in:
Meta Vens 2026-10-14 13:35:45 +02:00
commit 9282836fa5
15 changed files with 2596 additions and 866 deletions
Download patch file Download diff file Expand all files Collapse all files

9
packages/cultmath/README.md
View file

@ -78,6 +78,15 @@ shader semantics.
then `CultMath.Phacelle.glsl` only if it calls `cultmath_phacelle` and
`CultMath.Cycles.glsl` only if it calls a `cultmath_cycles_` function, then its
own shader source.
- Span primitives (`math.Span.cs`, `cultmath_*` in HLSL and GLSL): `snoise_hessian`,
the analytic Hessian of `snoise(float3)` on the same lattice, hash and tie order as
`snoise_grad` (a symmetric `float3x3`; GLSL `mat3`), and `integral_clamp_linear(a,
b, lo, hi, t0, t1)`, the closed-form integral over `[t0, t1]` of `clamp(a + b t,
lo, hi)`. Three measured bounds describe `snoise`: `SNOISE_LIPSCHITZ` (gradient
length, `NoiseBoundTests.MeasureLipschitz`), `SNOISE_HESSIAN` (Hessian spectral
norm, `MeasureHessian`) and `SNOISE_THIRD` (third-derivative operator norm,
`MeasureThird`); each is 1.10 times a seeded maximum refined by ascent, and a
fast test pins it. They are empirical with a margin, not proofs.
- `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.

6
packages/cultmath/shaders/CultMath.Affine.hlsl
View file

@ -10,7 +10,11 @@
// An upper bound on the spectral norm of the Hessian of cultmath_snoise(float3); provenance in
// math.Affine.cs (SNOISE_HESSIAN).
static const float CULTMATH_SNOISE_HESSIAN = 56.050385;
static const float CULTMATH_SNOISE_HESSIAN = 55.642143;
// 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 = 536.9152;
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); }

97
packages/cultmath/shaders/CultMath.glsl
View file

@ -228,6 +228,97 @@ vec4 cultmath_snoise_grad(vec3 value)
return vec4(grad, value2);
}
// Analytic Hessian of cultmath_snoise(vec3); 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. The matrix is
// built from nine scalars, which read as rows in HLSL and as columns in GLSL; it is symmetric, so the
// two agree.
mat3 cultmath_snoise_hessian(vec3 value)
{
vec2 c = vec2(1.0 / 6.0, 1.0 / 3.0);
vec3 i = floor(value + dot(value, c.yyy));
vec3 x0 = value - i + dot(i, c.xxx);
vec3 g = vec3(step(x0.y, x0.x), step(x0.z, x0.y), 1.0 - step(x0.z, x0.x));
vec3 l = 1.0 - g;
vec3 i1 = min(g.xyz, l.zxy);
vec3 i2 = max(g.xyz, l.zxy);
vec3 x1 = x0 - i1 + c.xxx;
vec3 x2 = x0 - i2 + c.yyy;
vec3 x3 = x0 - 0.5;
i = cultmath_snoise_mod289(i);
vec4 p = cultmath_snoise_permute(cultmath_snoise_permute(cultmath_snoise_permute(
i.z + vec4(0.0, i1.z, i2.z, 1.0)) + i.y + vec4(0.0, i1.y, i2.y, 1.0)) + i.x + vec4(0.0, i1.x, i2.x, 1.0));
const float n_ = 0.142857142857;
vec3 ns = vec3(2.0 * n_, 0.5 * n_ - 1.0, n_);
vec4 j = p - 49.0 * floor(p * ns.z * ns.z);
vec4 x_ = floor(j * ns.z);
vec4 y_ = floor(j - 7.0 * x_);
vec4 x = x_ * ns.x + ns.yyyy;
vec4 y = y_ * ns.x + ns.yyyy;
vec4 h = 1.0 - abs(x) - abs(y);
vec4 b0 = vec4(x.xy, y.xy);
vec4 b1 = vec4(x.zw, y.zw);
vec4 s0 = floor(b0) * 2.0 + 1.0;
vec4 s1 = floor(b1) * 2.0 + 1.0;
vec4 sh = -step(h, vec4(0.0, 0.0, 0.0, 0.0));
vec4 a0 = b0.xzyw + s0.xzyw * sh.xxyy;
vec4 a1 = b1.xzyw + s1.xzyw * sh.zzww;
vec3 p0 = vec3(a0.xy, h.x);
vec3 p1 = vec3(a0.zw, h.y);
vec3 p2 = vec3(a1.xy, h.z);
vec3 p3 = vec3(a1.zw, h.w);
vec4 norm = 1.79284291400159 - 0.85373472095314 * vec4(dot(p0, p0), dot(p1, p1), dot(p2, p2), dot(p3, p3));
p0 *= norm.x;
p1 *= norm.y;
p2 *= norm.z;
p3 *= norm.w;
vec4 m0 = max(0.5 - vec4(dot(x0, x0), dot(x1, x1), dot(x2, x2), dot(x3, x3)), 0.0);
vec4 m2 = m0 * m0;
vec4 m3 = m2 * m0;
vec4 px = vec4(dot(p0, x0), dot(p1, x1), dot(p2, x2), dot(p3, x3));
vec4 wide = 48.0 * m2 * px;
vec4 narrow = 8.0 * m3;
vec4 cxx = vec4(x0.x, x1.x, x2.x, x3.x);
vec4 cxy = vec4(x0.y, x1.y, x2.y, x3.y);
vec4 cxz = vec4(x0.z, x1.z, x2.z, x3.z);
vec4 cpx = vec4(p0.x, p1.x, p2.x, p3.x);
vec4 cpy = vec4(p0.y, p1.y, p2.y, p3.y);
vec4 cpz = vec4(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 mat3(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.
vec4 cultmath_fbm_grad(vec3 p, int octaves, float lacunarity, float gain)
@ -653,7 +744,11 @@ vec4 cultmath_iv_frustum_ball(vec2 centreSlope, float z0, float z1, float footpr
// An upper bound on the spectral norm of the Hessian of cultmath_snoise(vec3); provenance in
// math.Affine.cs (SNOISE_HESSIAN).
const float CULTMATH_SNOISE_HESSIAN = 56.050385;
const float CULTMATH_SNOISE_HESSIAN = 55.642143;
// An upper bound on the operator norm of the third-derivative tensor of cultmath_snoise(vec3), the
// next order after CULTMATH_SNOISE_HESSIAN; provenance in math.Span.cs (SNOISE_THIRD).
const float CULTMATH_SNOISE_THIRD = 536.9152;
vec3 cultmath_af_point(float value) { return vec3(value, 0.0, 0.0); }
vec3 cultmath_af_symbol(float x0, float a) { return vec3(x0, a, 0.0); }

91
packages/cultmath/shaders/CultMath.hlsl
View file

@ -228,6 +228,97 @@ 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. The matrix is
// built from nine scalars, which read as rows in HLSL and as columns in GLSL; it is symmetric, so the
// two agree.
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)

11
packages/cultmath/src/CultMath/math.Affine.cs
View file

@ -21,9 +21,10 @@ public static partial class math
// plain function of float values, so it mirrors into the HLSL common subset.
// An upper bound on the spectral norm of the Hessian of snoise(float3). Provenance:
// NoiseBoundTests.MeasureHessian (slow, explicit) takes the largest Hessian norm (central
// differences of snoise_grad, largest absolute eigenvalue) over 1e6 seeded points, refines the 1e4
// largest by ascent, and multiplies the refined maximum by 1.10;
// NoiseBoundTests.MeasureHessian (slow, explicit) takes the largest Hessian norm (snoise_hessian,
// largest absolute eigenvalue) over 1e6 seeded points, refines the 1e4 largest by ascent, and
// multiplies the refined maximum by 1.10 (50.5838 measured, 55.642143 pinned; the central
// differences of snoise_grad it replaced read 50.955 here, their roundoff at h = 1e-3);
// NoiseBoundTests.HessianConstantPinsSampledCurvature pins this value against a fast re-measure and
// the ascent from MeasureHessian's witness. Empirical with a margin, not a proof. Random balls
// cannot see a constant 20% low, so the defence is those two curvature tests plus
@ -31,8 +32,8 @@ 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.
public const float SNOISE_HESSIAN = 56.050385f;
// kernel must re-run MeasureHessian and re-pin it. SNOISE_THIRD (math.Span.cs) is the next-order pin.
public const float SNOISE_HESSIAN = 55.642143f;
/// <summary>The point form (x, 0, 0): the value x at every eps. Exact.</summary>
public static float3 af_point(float value) => new(value, 0.0f, 0.0f);

123
packages/cultmath/src/CultMath/math.Span.cs Normal file
View file

@ -0,0 +1,123 @@
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. Measured 488.1047 (a dense sweep of 20,000 directions at the witness
// found 99.9998% of the tensor norm there). 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 = 536.9152f;
// 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;
}
}

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

@ -0,0 +1,2 @@
fileFormatVersion: 2
guid: 1466848dfbfc420ab4f3f67090944a67

7
packages/cultmath/tests/CultMath.Tests/GlslLowering.cs
View file

@ -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)\)");

14
packages/cultmath/tests/CultMath.Tests/GlslMirrorTests.cs
View file

@ -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<object> 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 ----

7
packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs
View file

@ -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);
}
/// <summary>
@ -377,6 +379,9 @@ public sealed class HlslSourceCompatibilityTests
/// </summary>
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();

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

@ -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)));
@ -617,7 +605,7 @@ public sealed class NoiseBoundTests
}
// The start MeasureHessian printed as its witness.
private static readonly float3 HessianWitness = new(-245.19885f, -53.01323f, -91.47363f);
private static readonly float3 HessianWitness = new(-13.7920685f, -9.672134f, 61.563934f);
/// <summary>
/// Pins SNOISE_HESSIAN: no Hessian norm over 1e5 seeded points exceeds it, and the ascent from
@ -635,6 +623,181 @@ 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 the six of 26 starts (axes, edge midpoints, cube corners) with
// the largest |T(u, u, u)|, 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.OrderByDescending(v => Math.Abs(Cubic(t, v))).Take(6))
{
var u = (double[])start.Clone();
var sign = Cubic(t, u) < 0.0 ? -1.0 : 1.0;
for (var iteration = 0; iteration < 40; 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);
}
/// <summary>
/// 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.
/// </summary>
[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(-168.3755f, -3.8205261f, -109.16948f);
/// <summary>
/// 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.
/// </summary>
[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");
}
/// <summary>
/// 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

341
packages/cultmath/tests/CultMath.Tests/SpanPrimitiveTests.cs Normal file
View file

@ -0,0 +1,341 @@
using CultMath;
using Xunit;
using static CultMath.math;
namespace CultMath.Tests;
/// <summary>
/// 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.
/// </summary>
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<float3> TiePoints()
{
const int vertices = 160;
const int perVertex = 8;
var random = new System.Random(0x71E5);
var points = new List<float3>();
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<float3> 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;
var outside = 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: (4 S(h/2) - S(h)) / 3 is O(h^4). Two such
// extrapolations, from (h, h/2) and from (h/2, h/4), differ by about fifteen times the
// error of the finer one, so the comparison is against the finer and its tolerance is
// four times that gap, plus the roundoff of snoise (a few ulp) over (h/4)^2. Near a
// kernel's edge (m = 0) snoise's fourth derivative jumps and the extrapolation loses its
// order: such an entry may exceed the tight tolerance, but never by more than the gap
// between the last two second differences (an error that shrinks at least linearly in h
// is at most that gap), and the tight tolerance must hold for all but a fraction.
const float h = 1.0f / 16.0f;
var s0 = SecondDifference(p, i, j, h);
var s1 = SecondDifference(p, i, j, h / 2.0f);
var s2 = SecondDifference(p, i, j, h / 4.0f);
var coarseExtrapolation = (4.0 * s1 - s0) / 3.0;
var extrapolated = (4.0 * s2 - s1) / 3.0;
var tolerance = 4.0 * Math.Abs(extrapolated - coarseExtrapolation) + 64.0 * Unit / ((h / 4.0) * (h / 4.0));
var loose = tolerance + 2.0 * Math.Abs(s1 - s2);
var error = Math.Abs(hessian[i][j] - extrapolated);
worst = Math.Max(worst, error / tolerance);
largest = Math.Max(largest, Math.Abs(hessian[i][j]));
compared++;
if (error > tolerance)
outside++;
Assert.True(error <= loose, $"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 {loose:R}");
}
}
output.WriteLine($"snoise_hessian vs Richardson second differences of snoise over 4000 seeded points: {compared} entries, {outside} over the tight tolerance, worst error {worst:F3} of it, largest |H| {largest:R}");
Assert.True(largest > 10.0, $"the sample never meets curvature above 10 (largest {largest:R})");
Assert.True(outside <= compared / 200, $"{outside} of {compared} entries are over the tight Richardson tolerance, more than 0.5%");
}
[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<double, double> 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));
}
}

2458
packages/cultmath/tests/CultMath.Tests/fixtures/glsl-parity.json
View file

File diff suppressed because it is too large Load diff

6
packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Affine.hlsl
View file

@ -10,7 +10,11 @@
// An upper bound on the spectral norm of the Hessian of cultmath_snoise(float3); provenance in
// math.Affine.cs (SNOISE_HESSIAN).
static const float CULTMATH_SNOISE_HESSIAN = 56.050385;
static const float CULTMATH_SNOISE_HESSIAN = 55.642143;
// 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 = 536.9152;
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); }

91
packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl
View file

@ -228,6 +228,97 @@ 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. The matrix is
// built from nine scalars, which read as rows in HLSL and as columns in GLSL; it is symmetric, so the
// two agree.
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)