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-stage-functions (stage-functions + stage-shaders)

Browse source 
Span stages in C#, HLSL and GLSL: SpanPath/SpanStage inputs, float-rounding allowances on curvature and fold, zero-length flat spans, pinned shader guards, exact noise pins, compile-hlsl.ps1 with -WX, struct-aware shader mirror. Soul verdicts cut-stage-functions.s2-s4 and cut-stage-shaders.s1 hold; main's interim changes are cultcache-rs and docs only, disjoint from the verified C#/shader tree.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Meta Vens 2026-10-16 04:12:29 +02:00
commit d7d1bcf38e
27 changed files with 3456 additions and 102 deletions
Download patch file Download diff file Expand all files Collapse all files

9
packages/cultmath/README.md
View file

@ -82,7 +82,11 @@ shader semantics.
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
lo, hi)`; and the octave stages `span_fbm(path, stage,
octaves, lacunarity, gain)` and `span_ridged` (shader twins `cultmath_span_fbm` and `cultmath_span_ridged`,
whose callers build the `SpanStage` themselves), over a `SpanPath` and a `SpanStage` (`SpanStage.Identity`
is the plain fbm in C#), which return a `SpanNoise` (`field` as `float4(gradient, value)`,
`curvature`, `fold`) for a stretched, offset, shifted octave sum over a span. 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
@ -164,6 +168,7 @@ CultMath keeps portable DXC outside git under `packages/cultmath/.tools/`:
```powershell
.\packages\cultmath\tools\get-dxc.ps1
.\packages\cultmath\tools\compile-hlsl.ps1
.\packages\cultmath\tools\get-spirv-cross.ps1
.\packages\cultmath\tools\compile-hlsl-spirv.ps1 `
-ShaderPath E:\Projects\Odin\crates\muninn-move-tracker\shaders\MoveSphereCandidate.comp.hlsl `
@ -174,6 +179,8 @@ CultMath keeps portable DXC outside git under `packages/cultmath/.tools/`:
-OutputPath E:\Projects\Odin\crates\muninn-move-tracker\artifacts\shader\MoveSphereCandidate.comp.metal
```
`compile-hlsl.ps1` compiles both `CultMath.hlsl` copies as a `lib_6_3` library with the x64 DXC; a warning fails it too (`-WX`).
DXC emits Vulkan SPIR-V with `-spirv`. `-TargetEnv vulkan1.2` is exposed by the
script, but the current Windows DXC release compiled this shader successfully
with DXC's default Vulkan target while rejecting the dotted target-env value.

102
packages/cultmath/docs/design.md
View file

@ -112,8 +112,10 @@ local) and writes through a swizzle component (`v.xy.x = s`, which C# rejects as
modifying a temporary).
The HLSL side is checked too: `tools/get-dxc.ps1` fetches Microsoft's
DirectXShaderCompiler, and
`dxc -T lib_6_3 -HV 2021 shaders/CultMath.hlsl` must compile cleanly.
DirectXShaderCompiler (x64), and `tools/compile-hlsl.ps1` compiles both copies of
`CultMath.hlsl` with it (`dxc -T lib_6_3 -HV 2021 -WX`, so a warning fails like an error and the
text must compile cleanly); the script throws when a shader file is missing, when dxc cannot
launch or when it exits non-zero, so a stale exit code never reads as a pass.
The mirror test compares bit patterns (any NaN equals any NaN; -0 differs from
0). It proves the text of `CultMath.hlsl` computes the same float32 results as
@ -123,13 +125,20 @@ compiled under `using static CultMath.math;`, so its intrinsics are C# `math`
itself: the test proves the composition in the text, not the intrinsic rules,
which only `HlslSemanticsTests` pins.
Every mirror function is compared, internals included: its counterpart is the `math`
function of the same name, public or private (`snoise_mod289` and `snoise_permute` are
private), and a function with no counterpart fails by name. A struct parameter maps to the
CultMath type of the same name, and only if both have the same instance fields by name, type
and order (`SpanPath`, `SpanStage`); the argument is built from the case's scalars field by
field and copied to its C# twin by field name.
The comparison also allows exactly one struct return shape: a mirror function
whose HLSL return type is a plain struct of `float`/`floatN`/`int`/`intN`
fields is compared to its C# counterpart field by field, recursing into each
field's own components, down to bit-for-bit scalars, matching each leaf's
declaration path (field names, not just position) and declared type as well
as its value. `CultCellular` (`math.cs`, `cellular`) and `CultPhasor`
(`math.Phacelle.cs`, `phacelle`) are its consumers. A struct return with a mismatched, reordered-by-name, or
as its value. `CultCellular` (`math.cs`, `cellular`), `CultPhasor`
(`math.Phacelle.cs`, `phacelle`) and `SpanNoise` (`math.Span.cs`, `span_fbm`) are its consumers. A struct return with a mismatched, reordered-by-name, or
retyped field does not get a second comparison path or a shape of its own;
it fails the same walk that already handles a bare vector return.
@ -777,11 +786,18 @@ differentiating each corner's gradient `-8 m^3 s x + m^4 p` once more. Each
off-diagonal entry is computed once and placed twice, so the matrix equals its
transpose bit for bit, and a clamped corner (`m = 0`) contributes nothing. Built
from nine scalars, the matrix reads as columns in GLSL (`mat3`) and rows in HLSL,
which is the same matrix because it is symmetric. `SNOISE_THIRD = 536.9152`
which is the same matrix because it is symmetric. `SNOISE_THIRD = 535.8009`
(`math.Span.cs`) is the next order: 1.10 times the maximum, over 1e6 seeded points,
of the operator norm of the third-derivative tensor of `snoise` (central
differences of `snoise_hessian` at h = 1e-3, then an ascent over unit vectors),
refined from the 1e4 largest; measured 488.1047 on the tied kernel.
of the operator norm of the third-derivative tensor of `snoise`, refined from the
1e4 largest; measured 487.0917 on the tied kernel. The tensor is closed form, in
double: `snoise_hessian`'s per-corner entry differentiated once more, on the
kernel's lattice, tie order and hash, and the test
`ThirdTensorMatchesDifferencesOfTheHessian` holds it against Richardson-extrapolated
differences of the double Hessian. The operator norm is a local ascent over unit
vectors (a dense sweep of 20,000 directions at the witness finds 99.9987% of it) and
the outer maximum a sampled search refined by ascent, so the constant is empirical
with a 10% margin, not a proof; an earlier pin from central differences of the
float Hessian read 488.1047, its roundoff.
`MeasureThird` takes the measurement and `ThirdConstantPinsSampled` pins it.
`integral_clamp_linear(a, b, lo, hi, t0, t1)` is the closed-form integral over
`[t0, t1]` of `clamp(a + b t, lo, hi)`. All three have HLSL and GLSL lowerings and
@ -799,6 +815,76 @@ constant of 40 misses there). The margin over the function itself: a scan of the
independently over every hash configuration gives 52.22, so 55.64 sits 7% above the
configuration bound. A change to the kernels must re-run the measurements.
### Span stages
`span_fbm` and `span_ridged` (`math.Span.cs`; `cultmath_span_fbm` and `cultmath_span_ridged` in HLSL over
the same struct types, `cultmath_span_stage` behind both, and GLSL derived) are the octave stages
of a warp chain, `span_fbm(path, stage, octaves, lacunarity, gain)`, each returning
one `SpanNoise { float4 field; float curvature; float fold; }`. The inputs are two
structs with named public fields and no positional constructor, so the three floats
of the span and the six fields of the transform cannot be swapped silently.
`SpanPath { midpoint, velocity, speed, acceleration, halfLength }` is the path
`q(t)`, `|t| <= halfLength`, through `midpoint` with `velocity` at `t = 0`,
`|q'| <= speed` and `|q''| <= acceleration` along it. `SpanStage { frequency,
stretch, axis, shift, offset, octaveStep }` makes octave `o` sample `snoise_grad`
at `x_o(q) = f_o (M q + shift) + offset + o octaveStep`, with `f_o =
lacunarity^o`, amplitude `gain^o`, and `M = frequency (I + (1/stretch - 1) a a^T)`,
`a = normalize(axis)`; a zero axis or `stretch` 1 means no stretch, and `stretch`
must be positive and finite (unguarded, like `lacunarity`). `default(SpanStage)`
has frequency 0; `SpanStage.Identity` is the stage whose field equals `fbm_grad`
and `ridged_grad` bit for bit. `octaves`, `lacunarity` and `gain` stay positional:
they are `fbm_grad`'s own trailing triple, and `octaves` is the one int among the
scalars. `field` is the octave sum and its gradient in `q` (`f_o M grad n_o`).
`span_ridged` sums `gain^o (1 - |n_o|)`, `ridged_grad`'s convention. The Span
Lab's stage is `1 - 2|n|` per octave: with `G` the sum of `gain^o` over the
octaves (signed, so the alternating sum for a negative gain), its value is `2 value
- G`, and its gradient, curvature and fold are twice `span_ridged`'s.
The shaders have no `SpanStage.Identity`: a shader caller builds every stage field itself (in GLSL,
`SpanStage(1.0, 1.0, vec3(0.0), vec3(0.0), vec3(0.0), vec3(0.0))` is the identity). A constant would
need a second transformation rule (C# has no const structs, GLSL wants a constructor) and a mirror
exemption, for one value; the mirror test proves instead that the shader text at a C#-built identity
returns `cultmath_fbm_grad` and `cultmath_ridged_grad`'s field bit for bit.
Bounds come from the three pins only. With `|M| = |frequency| max(1, 1/stretch)` and
`nu = f_o |M| speed` the noise-space speed, the Hessian of the noise on the span has
operator norm at most `Hl = min(SNOISE_HESSIAN, ||snoise_hessian(x_mid)||_F + e_h +
SNOISE_THIRD nu halfLength)`, so `curvature = sum gain^o f_o^2 |M|^2 Hl` bounds the
operator norm of the field's Hessian in `q` at every point of the span (the pin takes
over once the third-derivative term exceeds it, and at `halfLength` 0 the local tier
alone remains). For `span_ridged`, an octave leaves the branch it is on at the
midpoint (`1 - s n`, `s` the midpoint's sign of `n`) only where `n` changes sign, and
there `|n| <= V - |n_mid|`, with `V = (|grad n_mid . w| + e_g |w|) halfLength + (Hl nu |w|
+ Gl f_o |M| acceleration) halfLength^2 / 2 + 2 e_n`, `w = f_o M velocity` and `Gl =
min(SNOISE_LIPSCHITZ, |grad n_mid| + e_g + Hl nu halfLength)`; `fold = sum 2 gain^o
max(V - |n_mid|, 0)`, exactly 0 for `span_fbm` and for an octave that cannot reach zero.
The float32 `snoise`, `snoise_grad` and `snoise_hessian` are not the exact functions, so
the quantities read at the midpoint carry allowances, by the derivation of `af_snoise`'s
comment. `S_o = |f_o| (|M q|_1 + |shift|_1) + |offset|_1 + |o octaveStep|_1 + sqrt(3) nu
halfLength` is the sum of the magnitudes the octave coordinate's roundings are relative
to, and `e_n = 2^-14 + L 2^-20 S_o` (value), `e_g = 2^-12 + H 2^-20 S_o` (gradient
length) and `e_h = 2^-10 + T 2^-20 S_o` (Frobenius norm of the Hessian), with `L`, `H`,
`T` the three pins; they add no pin. Measured against a double kernel with the float
hash (6000 points per band, `|x|` from 2^-2 to 2^23), two samples observed the error's
maxima as fractions of `e_n`, `e_g` and `e_h`: 0.057, 0.046 and 0.075 in the cut's own
sample (seed not recorded; the 0.075 is against `T 2^-20 |x|_1` alone, at the smallest
band, where the absolute Hessian error is 2.3e-5 and the absolute term carries it), and
0.064, 0.0485 and 0.057 in an audit's (seed 0x501). They are samples, not bounds: a new
seed moves them. They say the allowances enclose the error with a margin of about 15.
Errors grow linearly in `|x|`: the value's is 2e-4 at 2^8, 1.5e-2 at 2^14 and 0.94 at 2^20.
Rounding in the Frobenius-to-operator-norm step needs no allowance: the Frobenius norm
exceeds the operator norm.
Domain. No coordinate range is checked, clamped or refused. The allowances grow with
`S_o`, so `curvature` saturates at `SNOISE_HESSIAN` (true everywhere) by `S_o` of about
1.1e5, and `fold` turns positive for every octave whose `V + 2 e_n` reaches `|n_mid|`:
the bounds stay sound at every magnitude and stop being informative where float32 noise
is (see "Precision domain" above). `SpanStageTests` encloses the exact (double) Hessian
norm and both the float and the exact branch departure at 64 points of seeded spans in
bands of `S_o` from under 1e2 to over 1e6, and in a fixed span where the local tier
with its third-derivative term halved fails.
The slice. `af_frustum_ball(centreSlope, z0, z1, footprintPerDepth,
warpVariation)` returns the centre `(m_c z_m, z_m)` and a radius `z1
footprintPerDepth + warpVariation`, with `af_frustum_axis` the axis `(m_c, 1) h`:

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

@ -14,7 +14,7 @@ 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;
static const float CULTMATH_SNOISE_THIRD = 535.8009;
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); }

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

@ -302,10 +302,10 @@ mat3 cultmath_snoise_hessian(vec3 value)
// 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);
if (b == 0.0)
return t1 != t0 ? clamp(a, lo, hi) * (t1 - t0) : 0.0;
float direction = t1 < t0 ? -1.0 : 1.0;
float start = min(t0, t1);
float stop = max(t0, t1);
@ -315,7 +315,9 @@ float cultmath_integral_clamp_linear(float a, float b, float lo, float hi, float
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);
float held = enter > start ? before * (enter - start) : 0.0;
float kept = stop > leave ? after * (stop - leave) : 0.0;
float total = held + (a + b * (0.5 * (enter + leave))) * (leave - enter) + kept;
return direction * total;
}
@ -748,7 +750,7 @@ 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;
const float CULTMATH_SNOISE_THIRD = 535.8009;
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); }
@ -870,6 +872,119 @@ vec4 cultmath_af_frustum_ball(vec2 centreSlope, float z0, float z1, float footpr
// Span stages (invariant 8): the SpanPath, SpanStage and SpanNoise of math.Span.cs, field for field in
// C#'s order, and cultmath_span_stage, its octave stage, statement for statement so the text is a bit-for-bit
// mirror. They follow the includes because the bounds read CULTMATH_SNOISE_LIPSCHITZ, _HESSIAN and _THIRD.
// A shader caller builds a SpanStage itself (SpanStage(1.0, 1.0, vec3(0.0), vec3(0.0), vec3(0.0), vec3(0.0))
// is the identity of math.SpanStage.Identity, which has no shader form); see math.span_stage's comments for
// the span, the transform, the bounds and the allowances.
struct SpanPath
{
vec3 midpoint;
vec3 velocity;
float speed;
float acceleration;
float halfLength;
};
struct SpanStage
{
float frequency;
float stretch;
vec3 axis;
vec3 shift;
vec3 offset;
vec3 octaveStep;
};
struct SpanNoise
{
vec4 field;
float curvature;
float fold;
};
// The octave stage behind cultmath_span_fbm and cultmath_span_ridged, as math.span_stage is behind span_fbm
// and span_ridged. The Frobenius norm of the Hessian is written out in C#'s order (GLSL has no _mRC member
// access); the C# local named step is octaveShift here (it would shadow GLSL's step()).
SpanNoise cultmath_span_stage(bool ridged, SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain)
{
vec3 q = path.midpoint;
float halfLength = path.halfLength;
octaves = clamp(octaves, 0, 16);
bool stretched = dot(stage.axis, stage.axis) > 0.0 && stage.stretch != 1.0;
vec3 direction = stretched ? normalize(stage.axis) : vec3(0.0, 0.0, 0.0);
float squeeze = stretched ? 1.0 / stage.stretch - 1.0 : 0.0;
float norm = abs(stage.frequency) * (stretched ? max(1.0, 1.0 / stage.stretch) : 1.0);
vec3 mq = stretched ? stage.frequency * (q + squeeze * dot(direction, q) * direction) : stage.frequency * q;
vec3 mv = stretched ? stage.frequency * (path.velocity + squeeze * dot(direction, path.velocity) * direction) : stage.frequency * path.velocity;
float amplitude = 1.0;
float scale = 1.0;
float value = 0.0;
vec3 gradient = vec3(0.0, 0.0, 0.0);
float curvature = 0.0;
float fold = 0.0;
for (int i = 0; i < octaves; i++)
{
vec3 octaveShift = float(i) * stage.octaveStep;
vec3 x = scale * (mq + stage.shift) + (stage.offset + octaveShift);
vec4 n = cultmath_snoise_grad(x);
vec3 g = vec3(n.x, n.y, n.z);
vec3 mg = stretched ? stage.frequency * (g + squeeze * dot(direction, g) * direction) : stage.frequency * g;
if (ridged)
{
float s = sign(n.w);
value += amplitude * (1.0 - abs(n.w));
gradient += -amplitude * scale * s * mg;
}
else
{
value += amplitude * n.w;
gradient += amplitude * scale * mg;
}
float weight = abs(amplitude);
float factor = abs(scale) * norm;
float nu = factor * path.speed;
float rounding = 9.5367431640625e-7 * (abs(scale) * (cultmath_csum(abs(mq)) + cultmath_csum(abs(stage.shift))) + cultmath_csum(abs(stage.offset)) + cultmath_csum(abs(octaveShift)) + 1.7320508 * nu * halfLength);
float valueSlack = 6.103515625e-5 + CULTMATH_SNOISE_LIPSCHITZ * rounding;
float gradientSlack = 2.44140625e-4 + CULTMATH_SNOISE_HESSIAN * rounding;
float hessianSlack = 9.765625e-4 + CULTMATH_SNOISE_THIRD * rounding;
mat3 hm = cultmath_snoise_hessian(x);
float frobenius = sqrt(hm[0][0] * hm[0][0] + hm[0][1] * hm[0][1] + hm[0][2] * hm[0][2] + hm[1][0] * hm[1][0] + hm[1][1] * hm[1][1] + hm[1][2] * hm[1][2]
+ hm[2][0] * hm[2][0] + hm[2][1] * hm[2][1] + hm[2][2] * hm[2][2]);
float hessianBound = min(CULTMATH_SNOISE_HESSIAN, frobenius + hessianSlack + CULTMATH_SNOISE_THIRD * nu * halfLength);
curvature += weight * factor * factor * hessianBound;
if (ridged)
{
vec3 w = scale * mv;
float reach = hessianBound * nu * halfLength;
float gl = min(CULTMATH_SNOISE_LIPSCHITZ, length(g) + gradientSlack + reach);
float travel = (abs(dot(g, w)) + gradientSlack * length(w)) * halfLength + (hessianBound * nu * length(w) + gl * factor * path.acceleration) * (0.5 * halfLength * halfLength) + 2.0 * valueSlack;
fold += 2.0 * weight * max(travel - abs(n.w), 0.0);
}
scale *= lacunarity;
amplitude *= gain;
}
SpanNoise result;
result.field = vec4(gradient, value);
result.curvature = curvature;
result.fold = fold;
return result;
}
SpanNoise cultmath_span_fbm(SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain)
{
return cultmath_span_stage(false, path, stage, octaves, lacunarity, gain);
}
SpanNoise cultmath_span_ridged(SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain)
{
return cultmath_span_stage(true, path, stage, octaves, lacunarity, gain);
}
float cultmath_value_noise(vec2 position)
{
vec2 cell = floor(position);

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

@ -302,10 +302,10 @@ float3x3 cultmath_snoise_hessian(float3 value)
// 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);
if (b == 0.0)
return t1 != t0 ? clamp(a, lo, hi) * (t1 - t0) : 0.0;
float direction = t1 < t0 ? -1.0 : 1.0;
float start = min(t0, t1);
float stop = max(t0, t1);
@ -315,7 +315,9 @@ float cultmath_integral_clamp_linear(float a, float b, float lo, float hi, float
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);
float held = enter > start ? before * (enter - start) : 0.0;
float kept = stop > leave ? after * (stop - leave) : 0.0;
float total = held + (a + b * (0.5 * (enter + leave))) * (leave - enter) + kept;
return direction * total;
}
@ -603,6 +605,119 @@ CultCellular cultmath_cellular(float3 p)
#include "CultMath.Cycles.hlsl"
#include "CultMath.Affine.hlsl"
// Span stages (invariant 8): the SpanPath, SpanStage and SpanNoise of math.Span.cs, field for field in
// C#'s order, and cultmath_span_stage, its octave stage, statement for statement so the text is a bit-for-bit
// mirror. They follow the includes because the bounds read CULTMATH_SNOISE_LIPSCHITZ, _HESSIAN and _THIRD.
// A shader caller builds a SpanStage itself (SpanStage(1.0, 1.0, vec3(0.0), vec3(0.0), vec3(0.0), vec3(0.0))
// is the identity of math.SpanStage.Identity, which has no shader form); see math.span_stage's comments for
// the span, the transform, the bounds and the allowances.
struct SpanPath
{
float3 midpoint;
float3 velocity;
float speed;
float acceleration;
float halfLength;
};
struct SpanStage
{
float frequency;
float stretch;
float3 axis;
float3 shift;
float3 offset;
float3 octaveStep;
};
struct SpanNoise
{
float4 field;
float curvature;
float fold;
};
// The octave stage behind cultmath_span_fbm and cultmath_span_ridged, as math.span_stage is behind span_fbm
// and span_ridged. The Frobenius norm of the Hessian is written out in C#'s order (GLSL has no _mRC member
// access); the C# local named step is octaveShift here (it would shadow GLSL's step()).
SpanNoise cultmath_span_stage(bool ridged, SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain)
{
float3 q = path.midpoint;
float halfLength = path.halfLength;
octaves = clamp(octaves, 0, 16);
bool stretched = dot(stage.axis, stage.axis) > 0.0 && stage.stretch != 1.0;
float3 direction = stretched ? normalize(stage.axis) : float3(0.0, 0.0, 0.0);
float squeeze = stretched ? 1.0 / stage.stretch - 1.0 : 0.0;
float norm = abs(stage.frequency) * (stretched ? max(1.0, 1.0 / stage.stretch) : 1.0);
float3 mq = stretched ? stage.frequency * (q + squeeze * dot(direction, q) * direction) : stage.frequency * q;
float3 mv = stretched ? stage.frequency * (path.velocity + squeeze * dot(direction, path.velocity) * direction) : stage.frequency * path.velocity;
float amplitude = 1.0;
float scale = 1.0;
float value = 0.0;
float3 gradient = float3(0.0, 0.0, 0.0);
float curvature = 0.0;
float fold = 0.0;
for (int i = 0; i < octaves; i++)
{
float3 octaveShift = (float)i * stage.octaveStep;
float3 x = scale * (mq + stage.shift) + (stage.offset + octaveShift);
float4 n = cultmath_snoise_grad(x);
float3 g = float3(n.x, n.y, n.z);
float3 mg = stretched ? stage.frequency * (g + squeeze * dot(direction, g) * direction) : stage.frequency * g;
if (ridged)
{
float s = sign(n.w);
value += amplitude * (1.0 - abs(n.w));
gradient += -amplitude * scale * s * mg;
}
else
{
value += amplitude * n.w;
gradient += amplitude * scale * mg;
}
float weight = abs(amplitude);
float factor = abs(scale) * norm;
float nu = factor * path.speed;
float rounding = 9.5367431640625e-7 * (abs(scale) * (cultmath_csum(abs(mq)) + cultmath_csum(abs(stage.shift))) + cultmath_csum(abs(stage.offset)) + cultmath_csum(abs(octaveShift)) + 1.7320508 * nu * halfLength);
float valueSlack = 6.103515625e-5 + CULTMATH_SNOISE_LIPSCHITZ * rounding;
float gradientSlack = 2.44140625e-4 + CULTMATH_SNOISE_HESSIAN * rounding;
float hessianSlack = 9.765625e-4 + CULTMATH_SNOISE_THIRD * rounding;
float3x3 hm = cultmath_snoise_hessian(x);
float frobenius = sqrt(hm[0][0] * hm[0][0] + hm[0][1] * hm[0][1] + hm[0][2] * hm[0][2] + hm[1][0] * hm[1][0] + hm[1][1] * hm[1][1] + hm[1][2] * hm[1][2]
+ hm[2][0] * hm[2][0] + hm[2][1] * hm[2][1] + hm[2][2] * hm[2][2]);
float hessianBound = min(CULTMATH_SNOISE_HESSIAN, frobenius + hessianSlack + CULTMATH_SNOISE_THIRD * nu * halfLength);
curvature += weight * factor * factor * hessianBound;
if (ridged)
{
float3 w = scale * mv;
float reach = hessianBound * nu * halfLength;
float gl = min(CULTMATH_SNOISE_LIPSCHITZ, length(g) + gradientSlack + reach);
float travel = (abs(dot(g, w)) + gradientSlack * length(w)) * halfLength + (hessianBound * nu * length(w) + gl * factor * path.acceleration) * (0.5 * halfLength * halfLength) + 2.0 * valueSlack;
fold += 2.0 * weight * max(travel - abs(n.w), 0.0);
}
scale *= lacunarity;
amplitude *= gain;
}
SpanNoise result;
result.field = float4(gradient, value);
result.curvature = curvature;
result.fold = fold;
return result;
}
SpanNoise cultmath_span_fbm(SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain)
{
return cultmath_span_stage(false, path, stage, octaves, lacunarity, gain);
}
SpanNoise cultmath_span_ridged(SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain)
{
return cultmath_span_stage(true, path, stage, octaves, lacunarity, gain);
}
float cultmath_value_noise(float2 position)
{
float2 cell = floor(position);

29
packages/cultmath/src/CultMath/SpanNoise.cs Normal file
View file

@ -0,0 +1,29 @@
namespace CultMath;
/// <summary>
/// The result of <see cref="math.span_fbm"/> and <see cref="math.span_ridged"/> (invariant 8's one struct
/// return shape): an octave stage's value and gradient at the midpoint of a span, with the bounds that
/// hold over the whole span.
/// </summary>
public struct SpanNoise
{
/// <summary>(&#8711;, value): the stage's gradient in the span's space and its value, in the <c>float4(gradient.xyz, value.w)</c> layout of a bare value-and-gradient return.</summary>
public float4 field;
/// <summary>An upper bound, at every point the span can reach, on the operator norm of the field's Hessian in the span's space.</summary>
public float curvature;
/// <summary>
/// An upper bound on how far a ridged sum can leave its midpoint branch over the span (each octave's
/// 1 - |n| against the same octave continued with the midpoint's sign of n); 0 for
/// <see cref="math.span_fbm"/> and for any octave that cannot cross zero over the span.
/// </summary>
public float fold;
public SpanNoise(float4 field, float curvature, float fold)
{
this.field = field;
this.curvature = curvature;
this.fold = fold;
}
}

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

@ -0,0 +1,2 @@
fileFormatVersion: 2
guid: fe0da4c343e24c3781b076943f08ddbc

26
packages/cultmath/src/CultMath/SpanPath.cs Normal file
View file

@ -0,0 +1,26 @@
namespace CultMath;
/// <summary>
/// The span <see cref="math.span_fbm"/> and <see cref="math.span_ridged"/> bound over: the path q(t),
/// |t| &lt;= <see cref="halfLength"/>, through <see cref="midpoint"/> at t = 0 with <see cref="velocity"/> there,
/// |q'| &lt;= <see cref="speed"/> and |q''| &lt;= <see cref="acceleration"/> along it. speed is at least
/// |velocity|; speed, acceleration and halfLength are non-negative. Name every field: the three floats
/// have no order a caller could mistake for a convention.
/// </summary>
public struct SpanPath
{
/// <summary>q(0), the point the field and its gradient are evaluated at.</summary>
public float3 midpoint;
/// <summary>q'(0).</summary>
public float3 velocity;
/// <summary>An upper bound on |q'(t)| over the span; at least |velocity|.</summary>
public float speed;
/// <summary>An upper bound on |q''(t)| over the span.</summary>
public float acceleration;
/// <summary>The span covers |t| &lt;= halfLength; 0 is the midpoint alone.</summary>
public float halfLength;
}

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

@ -0,0 +1,2 @@
fileFormatVersion: 2
guid: 3a8c1d52e7b94f06a1d5c2b9e04f7a31

33
packages/cultmath/src/CultMath/SpanStage.cs Normal file
View file

@ -0,0 +1,33 @@
namespace CultMath;
/// <summary>
/// The transform of one octave stage of a warp chain, as <see cref="math.span_fbm"/> and
/// <see cref="math.span_ridged"/> read it. Octave o samples snoise at x_o = f_o (M q + shift) + offset +
/// o octaveStep, with f_o = lacunarity^o and M = <see cref="frequency"/> (I + (1/<see cref="stretch"/> - 1) a a^T),
/// a = normalize(<see cref="axis"/>); a zero axis, or a stretch of 1, means no stretch. stretch must be
/// positive and finite, and nothing guards it. <c>default(SpanStage)</c> has frequency 0, which samples
/// one point: start from <see cref="Identity"/>.
/// </summary>
public struct SpanStage
{
/// <summary>The stage whose span_fbm/span_ridged field equals fbm_grad/ridged_grad: frequency 1, stretch 1, every vector zero.</summary>
public static readonly SpanStage Identity = new SpanStage { frequency = 1.0f, stretch = 1.0f };
/// <summary>The scale of M; its sign is kept.</summary>
public float frequency;
/// <summary>The factor the axis direction is stretched by (M scales it by 1/stretch); 1 is none.</summary>
public float stretch;
/// <summary>The direction stretch applies along; any length, zero for none.</summary>
public float3 axis;
/// <summary>Added to M q before the per-octave scaling.</summary>
public float3 shift;
/// <summary>Added after the per-octave scaling.</summary>
public float3 offset;
/// <summary>Added to the octave coordinate once per octave index: o octaveStep.</summary>
public float3 octaveStep;
}

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

@ -0,0 +1,2 @@
fileFormatVersion: 2
guid: 9d46b0f27c1e4a85b3e8f1a60c2d5b97

168
packages/cultmath/src/CultMath/math.Span.cs
View file

@ -5,15 +5,18 @@ 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;
// the largest such norm over 1e6 seeded points of [-256, 256]^3 of the closed-form tensor (snoise_hessian's
// per-corner derivation differentiated once more, in double, with the kernel's lattice, tie order and
// hash; NoiseBoundTests.ThirdTensorMatchesDifferencesOfTheHessian pins it against differences of the
// double Hessian), 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:
// refined maximum 487.0917, a dense sweep of 20,000 directions at the witness's ascent end found
// 99.9987% of the tensor norm there. The inner maximisation is a local ascent and the outer a sampled
// search refined by ascent: 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 = 535.8009f;
// 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
@ -100,14 +103,17 @@ public static partial class math
// 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.
// clamps to hi everywhere, as clamp does, and is normalised (lo = min(lo, hi)) before anything reads it,
// so no clamp ever sees lo above hi (GLSL leaves that undefined). A crossing that overflows (b tiny) lies past the span and
// cuts nothing. A held stretch contributes only when its length is positive, so an infinite bound
// (relu has hi = +inf) times a zero length never forms: the integral of a ramp clamped on one side
// only is finite.
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);
if (b == 0.0f)
return t1 != t0 ? clamp(a, lo, hi) * (t1 - t0) : 0.0f;
var direction = t1 < t0 ? -1.0f : 1.0f;
var start = min(t0, t1);
var stop = max(t0, t1);
@ -117,7 +123,139 @@ public static partial class math
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);
var held = enter > start ? before * (enter - start) : 0.0f;
var kept = stop > leave ? after * (stop - leave) : 0.0f;
var total = held + (a + b * (0.5f * (enter + leave))) * (leave - enter) + kept;
return direction * total;
}
private static float span_frobenius(float3x3 m) =>
sqrt(m._m00 * m._m00 + m._m01 * m._m01 + m._m02 * m._m02 + m._m10 * m._m10 + m._m11 * m._m11 + m._m12 * m._m12
+ m._m20 * m._m20 + m._m21 * m._m21 + m._m22 * m._m22);
// Plain fbm of snoise over a span: the octave stage of a warp chain, its value and gradient at the
// span's midpoint and the bounds that hold over the whole span. See span_stage for the span, the
// transform and the bound rules. At SpanStage.Identity the field equals fbm_grad(path.midpoint,
// octaves, lacunarity, gain) bit for bit.
public static SpanNoise span_fbm(SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain) =>
span_stage(false, path, stage, octaves, lacunarity, gain);
// Ridged fbm over a span: octave o contributes gain^o (1 - |n_o|), ridged_grad's convention. fold
// bounds how far the sum can leave the branch it is on at the midpoint. At SpanStage.Identity the field
// equals ridged_grad bit for bit.
//
// The Span Lab's stage is 1 - 2|n| per octave. With G the sum of gain^o over the octaves (signed: a
// negative gain gives its alternating sum), the lab's value is 2 value - G, and its gradient,
// curvature and fold are twice this one's. A consumer in the lab form applies that mapping and does
// not derive its own bounds.
public static SpanNoise span_ridged(SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain) =>
span_stage(true, path, stage, octaves, lacunarity, gain);
// The span is SpanPath: q(t), |t| <= halfLength, through midpoint with velocity at t = 0, |q'| <= speed
// and |q''| <= acceleration along it (speed at least |velocity|; all three and halfLength non-negative).
// The transform is SpanStage: octave o samples snoise_grad at x_o(q) = f_o (M q + shift) + offset +
// o octaveStep with f_o = lacunarity^o, amplitude gain^o and M = frequency (I + (1/stretch - 1) a a^T),
// a = normalize(axis); a zero axis, or stretch 1, means no stretch. stretch must be positive and
// finite, and nothing guards it, as nothing guards lacunarity. octaves clamps to [0, 16] as fbm_grad's
// does. The field is the octave sum in q: gradient f_o M grad n_o (M is symmetric), the same octave
// loop as fbm_grad and ridged_grad.
//
// Bounds, per octave, from SNOISE_LIPSCHITZ (L), SNOISE_HESSIAN (H) and SNOISE_THIRD (T) only. With
// |M| = |frequency| max(1, 1/stretch) (1 without stretch) and nu = f_o |M| speed, the noise-space speed:
// * |x(t) - x_mid| <= nu halfLength, so the Hessian of n at x(t) has operator norm at most
// Hl = min(H, ||hessian(x_mid)||_F + e_h + T nu halfLength) (the Hessian is T-Lipschitz in operator
// norm, and the Frobenius norm bounds the operator norm), and |grad n(x(t))| at most
// Gl = min(L, |grad n_mid| + e_g + Hl nu halfLength).
// * curvature = sum gain^o f_o^2 |M|^2 Hl: the Hessian in q of n(x_o(q)) is f_o^2 M H M.
// * fold, ridged only. Over the span the octave's value departs from its midpoint branch
// gain^o (1 - s n) (s the sign of n at the midpoint) only where n has changed sign, and there
// |n(t)| <= V - |n_mid|, V bounding |n(t) - n_mid|. Writing n(t) - n_mid as the integral of
// grad n . x', x' = w + rho with w = f_o M velocity and |rho| <= f_o |M| acceleration |t|, and
// grad n = grad n_mid + delta with |delta| <= Hl nu |t|:
// V = (|grad n_mid . w| + e_g |w|) halfLength + (Hl nu |w| + Gl f_o |M| acceleration) halfLength^2 / 2
// + 2 e_n,
// and the octave's fold is 2 gain^o max(V - |n_mid|, 0), zero when it cannot reach zero.
//
// The float32 snoise, snoise_grad and snoise_hessian are not the exact functions, and the field a
// consumer reads is the float one, so what is read at the midpoint carries an allowance, by the
// derivation af_snoise's comment gives (math.Affine.cs): the skew and the corner offsets round the
// coordinate by about 2u of the magnitudes it is formed from, which moves a result by its derivative
// times that, and the caller's scaling of the point and a GPU's evaluation take 2u to 16u = 2^-20.
// S_o is the sum of those magnitudes, |f_o| (|M q|_1 + |shift|_1) + |offset|_1 + |o octaveStep|_1, plus
// the span's reach sqrt(3) nu halfLength (the coordinate of a path point). The allowances are the
// value's e_n = 2^-14 + L 2^-20 S_o, the gradient length's e_g = 2^-12 + H 2^-20 S_o (H is the
// gradient's own derivative's pin) and the Hessian's Frobenius norm's e_h = 2^-10 + T 2^-20 S_o (T is
// the Hessian's derivative's pin): af_snoise's absolute terms for the first two, and 2^-10 for the
// Hessian's, which the kernel's rounding at small coordinates (2.3e-5 absolute, measured against a
// double kernel) sits well under. Two samples of that probe (6000 points per band, |x| from 2^-2 to
// 2^23) observed the error's maxima, as fractions of e_n, e_g and e_h: 0.057, 0.046 and 0.075 in the
// cut's own sample (seed not recorded; the 0.075 is against T 2^-20 |x|_1 alone, at the smallest
// band), and 0.064, 0.0485 and 0.057 in an audit's (seed 0x501). They are samples, not bounds: a new
// seed moves them. What they say is that the allowances enclose the error with a margin of about 15.
// No new pin enters.
//
// Domain: no coordinate range is checked, clamped or refused. The allowances grow with S_o, so
// curvature saturates at H (true everywhere) by S_o of about 1.1e5, and fold turns positive for every
// octave whose V + 2 e_n reaches |n_mid|: the bounds stay sound at every magnitude and stop being
// informative where float32 noise is (design.md, Precision domain).
private static SpanNoise span_stage(bool ridged, SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain)
{
var q = path.midpoint;
var halfLength = path.halfLength;
octaves = clamp(octaves, 0, 16);
var stretched = dot(stage.axis, stage.axis) > 0.0f && stage.stretch != 1.0f;
var direction = stretched ? normalize(stage.axis) : new float3(0.0f, 0.0f, 0.0f);
var squeeze = stretched ? 1.0f / stage.stretch - 1.0f : 0.0f;
var norm = abs(stage.frequency) * (stretched ? max(1.0f, 1.0f / stage.stretch) : 1.0f);
var mq = stretched ? stage.frequency * (q + squeeze * dot(direction, q) * direction) : stage.frequency * q;
var mv = stretched ? stage.frequency * (path.velocity + squeeze * dot(direction, path.velocity) * direction) : stage.frequency * path.velocity;
var amplitude = 1.0f;
var scale = 1.0f;
var value = 0.0f;
var gradient = new float3(0.0f, 0.0f, 0.0f);
var curvature = 0.0f;
var fold = 0.0f;
for (var i = 0; i < octaves; i++)
{
var step = (float)i * stage.octaveStep;
var x = scale * (mq + stage.shift) + (stage.offset + step);
var n = snoise_grad(x);
var g = new float3(n.x, n.y, n.z);
var mg = stretched ? stage.frequency * (g + squeeze * dot(direction, g) * direction) : stage.frequency * g;
if (ridged)
{
var s = (float)sign(n.w);
value += amplitude * (1.0f - abs(n.w));
gradient += -amplitude * scale * s * mg;
}
else
{
value += amplitude * n.w;
gradient += amplitude * scale * mg;
}
var weight = abs(amplitude);
var factor = abs(scale) * norm;
var nu = factor * path.speed;
var rounding = 9.5367431640625e-7f * (abs(scale) * (csum(abs(mq)) + csum(abs(stage.shift))) + csum(abs(stage.offset)) + csum(abs(step)) + 1.7320508f * nu * halfLength);
var valueSlack = 6.103515625e-5f + SNOISE_LIPSCHITZ * rounding;
var gradientSlack = 2.44140625e-4f + SNOISE_HESSIAN * rounding;
var hessianSlack = 9.765625e-4f + SNOISE_THIRD * rounding;
var local = min(SNOISE_HESSIAN, span_frobenius(snoise_hessian(x)) + hessianSlack + SNOISE_THIRD * nu * halfLength);
curvature += weight * factor * factor * local;
if (ridged)
{
var w = scale * mv;
var reach = local * nu * halfLength;
var gl = min(SNOISE_LIPSCHITZ, length(g) + gradientSlack + reach);
var travel = (abs(dot(g, w)) + gradientSlack * length(w)) * halfLength + (local * nu * length(w) + gl * factor * path.acceleration) * (0.5f * halfLength * halfLength) + 2.0f * valueSlack;
fold += 2.0f * weight * max(travel - abs(n.w), 0.0f);
}
scale *= lacunarity;
amplitude *= gain;
}
return new SpanNoise(new float4(gradient, value), curvature, fold);
}
}

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

@ -45,6 +45,54 @@ public sealed class GlslMirrorTests
}
}
// Counting stand-ins for the mirror's clamp, one per overload, declared in the mirror's class so they
// hide math.clamp for every function in it. Each counts a call whose minimum exceeds its maximum in any
// component and forwards to math.clamp. (The int vector overloads only forward: no mirror function
// clamps integer vectors with swapped bounds, and counting them would be unreachable.)
private const string CountingClampShadows = @"
public static int ClampMinAboveMax;
static float clamp(float x, float a, float b) { if (a > b) ClampMinAboveMax++; return math.clamp(x, a, b); }
static float2 clamp(float2 x, float2 a, float2 b) { if (a.x > b.x || a.y > b.y) ClampMinAboveMax++; return math.clamp(x, a, b); }
static float3 clamp(float3 x, float3 a, float3 b) { if (a.x > b.x || a.y > b.y || a.z > b.z) ClampMinAboveMax++; return math.clamp(x, a, b); }
static float4 clamp(float4 x, float4 a, float4 b) { if (a.x > b.x || a.y > b.y || a.z > b.z || a.w > b.w) ClampMinAboveMax++; return math.clamp(x, a, b); }
static int clamp(int x, int a, int b) { if (a > b) ClampMinAboveMax++; return math.clamp(x, a, b); }
static int2 clamp(int2 x, int2 a, int2 b) { if (a.x > b.x || a.y > b.y) ClampMinAboveMax++; return math.clamp(x, a, b); }
static int3 clamp(int3 x, int3 a, int3 b) { if (a.x > b.x || a.y > b.y || a.z > b.z) ClampMinAboveMax++; return math.clamp(x, a, b); }
static int4 clamp(int4 x, int4 a, int4 b) { if (a.x > b.x || a.y > b.y || a.z > b.z || a.w > b.w) ClampMinAboveMax++; return math.clamp(x, a, b); }
";
/// <summary>
/// The lowering keeps HLSL's statement order (CommittedGlslEqualsLowering pins the GLSL as that
/// lowering), and GLSL leaves clamp(x, minVal, maxVal) undefined for minVal above maxVal, so no clamp
/// reached by cultmath_integral_clamp_linear may see lo above hi, on any path and for any spelling of the
/// normalisation. The mirror is compiled with counting clamps and the function is run over the edge grid
/// (lo above hi, infinite bounds, zero-length spans, b = 0 and b != 0): the count of such clamp calls
/// stays 0 after every call. The cases are HlslSourceCompatibilityTests.IntegralClampLinearCases: the grid fixes
/// a in {-2.5, 0, 1.5}, so a normalisation keyed on a past 1.75 never clamps wrongly there; the committed
/// fixture's random cases reach it.
/// </summary>
[Fact]
public void IntegralClampLinearNeverClampsWithMinAboveMax()
{
var (assembly, errors) = HlslSourceCompatibilityTests.CompileShaderMirror(CountingClampShadows);
Assert.True(assembly is not null, string.Join(Environment.NewLine, errors));
var shaderType = assembly!.GetType("CultMathHlsl.HlslShader")!;
var shader = Activator.CreateInstance(shaderType);
var count = shaderType.GetField("ClampMinAboveMax")!;
var function = shaderType.GetMethod("cultmath_integral_clamp_linear", BindingFlags.Instance | BindingFlags.NonPublic)!;
var bad = new List<string>();
foreach (var c in HlslSourceCompatibilityTests.IntegralClampLinearCases)
{
count.SetValue(null, 0);
function.Invoke(shader, c.Cast<object>().ToArray());
if ((int)count.GetValue(null)! > 0)
bad.Add($"({string.Join(", ", c)})");
}
Assert.True(bad.Count == 0, $"{bad.Count} cases clamp with the minimum above the maximum:{Environment.NewLine}{string.Join(Environment.NewLine, bad.Take(12))}");
}
// Identifier boundaries: cultmath_lerp is a CultMath function, lerp( is HLSL.
private static readonly Regex HlslOnly = new(
@"\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");
@ -638,8 +686,68 @@ public sealed class GlslMirrorTests
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) };
}),
// The edge cases the family above never draws: lo above hi (b = 0 and b != 0; the lowering must
// normalise before any clamp sees it) and infinite bounds (hi = +inf, lo = -inf, both).
new Family("float cultmath_integral_clamp_linear(float, float, float, float, float, float)", "ulp-bounded, lo above hi and infinite bounds",
r =>
{
var a = Uniform(r, -4.0f, 4.0f);
var slope = Uniform(r, -4.0f, 4.0f);
var (t0, t1) = (Uniform(r, -3.0f, 3.0f), Uniform(r, -3.0f, 3.0f));
return r.Next(5) switch
{
0 => new object[] { a, 0.0f, Uniform(r, 0.5f, 2.0f), Uniform(r, -2.0f, 0.4f), t0, t1 },
1 => new object[] { a, slope, Uniform(r, 0.5f, 2.0f), Uniform(r, -2.0f, 0.4f), t0, t1 },
2 => new object[] { a, slope, Uniform(r, -2.0f, 1.0f), float.PositiveInfinity, t0, t1 },
3 => new object[] { a, slope, float.NegativeInfinity, Uniform(r, -1.0f, 2.0f), t0, t1 },
_ => new object[] { a, slope, float.NegativeInfinity, float.PositiveInfinity, t0, t1 },
};
}),
// Zero-length spans (t1 = t0 three times in four) with a flat or sloped line, and infinite clamp bounds
// on either side or both: an infinite bound times a zero length must give +0, and a flat line clamped
// to an infinite bound over a nonzero span is infinite, never NaN.
new Family("float cultmath_integral_clamp_linear(float, float, float, float, float, float)", "ulp-bounded, zero-length spans and flat infinite clamps",
r =>
{
float Line() => r.Next(3) switch { 0 => 0.0f, 1 => -0.0f, _ => Uniform(r, -4.0f, 4.0f) };
float Bound() => r.Next(3) switch { 0 => float.NegativeInfinity, 1 => float.PositiveInfinity, _ => Uniform(r, -2.0f, 2.0f) };
var (a, b, lo, hi) = (Line(), Line(), Bound(), Bound());
var t0 = Uniform(r, -3.0f, 3.0f);
var t1 = r.Next(4) != 0 ? t0 : Uniform(r, -3.0f, 3.0f);
return new object[] { a, b, lo, hi, t0, t1 };
}),
// The span stages over their struct arguments: realistic stages and spans (SpanStageTests' own draw), the
// midpoint a Point3, one span in four a point.
new Family("SpanNoise cultmath_span_fbm(SpanPath, SpanStage, int, float, float)", "ulp-bounded", SpanArguments),
new Family("SpanNoise cultmath_span_ridged(SpanPath, SpanStage, int, float, float)", "ulp-bounded", SpanArguments),
}).ToArray();
private static object[] SpanArguments(System.Random r)
{
var stage = SpanStageTests.Stage.Draw(r);
var halfLength = r.Next(4) == 0 ? 0.0f : new[] { 0.02f, 0.1f, 0.5f }[r.Next(3)];
var span = SpanStageTests.Span.Draw(r, Point3(r), halfLength, Uniform(r, 0.0f, 2.0f));
return new object[] { span.Path, stage.S, stage.Octaves, stage.Lacunarity, stage.Gain };
}
/// <summary>
/// The committed fixture's cases of the three cultmath_integral_clamp_linear families, decoded into the
/// C# arguments (a, b, lo, hi, t0, t1): the random cases the clamp-normalisation tests also run.
/// </summary>
internal static IReadOnlyList<float[]> IntegralClampLinearFixtureCases
{
get
{
const string clamp = "float cultmath_integral_clamp_linear(float, float, float, float, float, float)";
var method = Resolve(clamp);
return JsonNode.Parse(File.ReadAllText(FixturePath))!["functions"]!.AsArray()
.Where(family => (string)family!["glsl"]! == clamp)
.SelectMany(family => family!["cases"]!.AsArray())
.Select(entry => Arguments(method, ((string)entry!).Split(" -> ")[0]).Select(a => (float)a).ToArray())
.ToArray();
}
}
// The fixture text, each case's arguments given by family and point index.
private static string GoldenFixture(Func<int, int, object[]> arguments)
{
@ -682,6 +790,7 @@ public sealed class GlslMirrorTests
{
["float"] = typeof(float), ["int"] = typeof(int), ["vec2"] = typeof(float2), ["vec3"] = typeof(float3),
["vec4"] = typeof(float4), ["ivec3"] = typeof(int3), ["ivec4"] = typeof(int4),
["SpanPath"] = typeof(SpanPath), ["SpanStage"] = typeof(SpanStage),
};
// The C# math method a GLSL signature names: cultmath_ dropped, parameter types mapped.

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

@ -39,21 +39,65 @@ public sealed class HlslSourceCompatibilityTests
Assert.NotNull(ShaderFunction(shader, "cultmath_af_snoise", typeof(float3), typeof(float3), typeof(float)));
Assert.NotNull(ShaderFunction(shader, "cultmath_af_fbm", typeof(float3), typeof(float3), typeof(float), typeof(int), typeof(float), typeof(float)));
Assert.NotNull(ShaderFunction(shader, "cultmath_af_frustum_ball", typeof(float2), typeof(float), typeof(float), typeof(float), typeof(float)));
var (spanPath, spanStage) = SpanTypes(shader);
Assert.NotNull(ShaderFunction(shader, "cultmath_span_fbm", spanPath, spanStage, typeof(int), typeof(float), typeof(float)));
Assert.NotNull(ShaderFunction(shader, "cultmath_span_ridged", spanPath, spanStage, typeof(int), typeof(float), typeof(float)));
Assert.NotNull(ShaderFunction(shader, "cultmath_span_stage", typeof(bool), spanPath, spanStage, typeof(int), typeof(float), typeof(float)));
}
// The mirror's struct types are private nested types of HlslShader (a struct without a modifier).
private static (Type Path, Type Stage) SpanTypes(Type shader) =>
(shader.GetNestedType("SpanPath", BindingFlags.Public | BindingFlags.NonPublic)!, shader.GetNestedType("SpanStage", BindingFlags.Public | BindingFlags.NonPublic)!);
// HLSL functions carry no access modifier, so they compile as private instance methods.
private static MethodInfo? ShaderFunction(Type shader, string name, params Type[] parameters) =>
shader.GetMethod(name, BindingFlags.Instance | BindingFlags.NonPublic, null, parameters, null);
// Mirror internals with no public C# math counterpart. Any other cultmath_* function must match one.
private static readonly string[] MirrorOnly = { "cultmath_snoise_mod289", "cultmath_snoise_permute" };
private static readonly float[] Specials = { 0.0f, -0.0f, 1.0f, -1.0f, 0.5f, -2.75f, 3.0f, 1.0e-7f, 1.0e4f, -1.0e4f, 1.0e7f, float.NaN };
// Shader constants with no C# twin of the same name (empty when every constant has one).
private static readonly string[] ConstantsWithoutTwin = Array.Empty<string>();
/// <summary>
/// Every mirror function with a C# math counterpart returns the same bits (any NaN equals any NaN;
/// Every static const of the compiled HLSL mirror (a static literal CULTMATH_* field, whatever its access) is the
/// same float32, bit for bit, as the math constant of the same name without the prefix (private ones
/// included), so a pin edited in one runtime and not the other fails here whichever constant it is.
/// A constant with neither a twin nor a listing in ConstantsWithoutTwin fails too.
/// </summary>
[Fact]
public void EveryShaderConstantEqualsItsCSharpTwin()
{
var (assembly, errors) = CompileShaderMirror();
Assert.True(assembly is not null, string.Join(Environment.NewLine, errors));
var shader = assembly!.GetType("CultMathHlsl.HlslShader")!;
var constants = shader.GetFields(BindingFlags.Public | BindingFlags.NonPublic | BindingFlags.Static)
.Where(f => f.IsLiteral && f.Name.StartsWith("CULTMATH_", StringComparison.Ordinal))
.ToArray();
Assert.True(constants.Length >= 10, $"only {constants.Length} shader constants were found by reflection");
foreach (var constant in constants)
{
var name = constant.Name["CULTMATH_".Length..];
var twin = typeof(math).GetField(name, BindingFlags.Public | BindingFlags.NonPublic | BindingFlags.Static);
if (twin is null)
{
Assert.True(ConstantsWithoutTwin.Contains(constant.Name), $"{constant.Name} has no math.{name} and is not listed in ConstantsWithoutTwin");
continue;
}
var shaderValue = (float)constant.GetRawConstantValue()!;
var mathValue = (float)twin.GetValue(null)!;
Assert.True(BitConverter.SingleToInt32Bits(shaderValue) == BitConverter.SingleToInt32Bits(mathValue),
$"{constant.Name} is {shaderValue:R} in the shader and math.{name} is {mathValue:R}");
}
}
// Special values every mirror argument takes: zeros, ties, negatives, large values, NaN and both infinities.
private static readonly float[] Specials = { 0.0f, -0.0f, 1.0f, -1.0f, 0.5f, -2.75f, 3.0f, 1.0e-7f, 1.0e4f, -1.0e4f, 1.0e7f, float.NaN, float.PositiveInfinity, float.NegativeInfinity };
/// <summary>
/// Every mirror function matches the math function of the same name, public or private (a function with no
/// twin fails, and nothing is exempt), with a struct parameter mapped by name and fields to the CultMath
/// type of the same name, and returns the same bits (any NaN equals any NaN;
/// -0 differs from 0) over random inputs, every special value in every argument (ties, zeros,
/// negatives, NaN, large values, integer lattice points), special pairs, and vectors whose first two
/// negatives, NaN, infinities, large values, integer lattice points), special pairs, and vectors whose first two
/// components tie (x0.x == x0.y in snoise). This proves the text of CultMath.hlsl computes the same
/// float32 results as C# math on the CPU. It does not prove GPU agreement: driver sin precision alone
/// makes cultmath_hash and value noise differ bit for bit on hardware. The mirror's intrinsics are C#
@ -116,20 +160,15 @@ public sealed class HlslSourceCompatibilityTests
foreach (var mirror in shaderType.GetMethods(BindingFlags.Instance | BindingFlags.NonPublic).Where(m => m.Name.StartsWith("cultmath_")))
{
var types = mirror.GetParameters().Select(p => p.ParameterType).ToArray();
var counterpart = typeof(math).GetMethod(mirror.Name["cultmath_".Length..], BindingFlags.Public | BindingFlags.Static, null, types, null);
if (counterpart is null)
{
Assert.Contains(mirror.Name, MirrorOnly);
continue;
}
var counterpart = Twin(shaderType, mirror.Name["cultmath_".Length..], types);
Assert.True(counterpart is not null, $"{mirror.Name} has no C# math function of the same name and parameter types, public or private");
compared.Add(mirror.Name);
foreach (var input in cases)
{
var slot = 0;
var args = types.Select(t => Scalars.Contains(t) ? Scalar(t, input(slot++))
: Activator.CreateInstance(t, Components(t).Select(f => Scalar(f.FieldType, input(slot++))).ToArray())!).ToArray();
var expected = Values(counterpart.Invoke(null, args)!);
var args = types.Select(t => Build(t, () => input(slot++))).ToArray();
var expected = Values(counterpart!.Invoke(null, args.Select(a => CopyFields(a, TwinType(shaderType, a.GetType())!)).ToArray())!);
var actual = Values(mirror.Invoke(shader, args)!);
if (!BitwiseEqual(expected, actual))
{
@ -139,7 +178,160 @@ public sealed class HlslSourceCompatibilityTests
}
Assert.True(mismatches.Count == 0, $"{mismatches.Count} mismatches:{Environment.NewLine}{string.Join(Environment.NewLine, mismatches.Take(12))}");
Assert.Equal(120, compared.Count);
Assert.Equal(125, compared.Count);
}
/// <summary>
/// The 1500 arguments (a, b, lo, hi, t0, t1) at which integral_clamp_linear's zero-length rule, its lo
/// normalisation and its infinite bounds all bite: a in {-2.5, 0, 1.5}, b in {-1.5, -0, 0, 2}, lo in
/// {-inf, -1, 0.5, 3, +inf}, hi in {-inf, -0.5, 1, 4, +inf}, and (t0, t1) in {(0, 0), (1.25, 1.25), (-1, 2),
/// (2, -1), (0, 1e-30)}. One spelling, read by the mirror grid test and by the GLSL clamp shadow test.
/// </summary>
internal static IReadOnlyList<float[]> IntegralClampLinearEdgeGrid { get; } = BuildIntegralClampLinearEdgeGrid();
private static List<float[]> BuildIntegralClampLinearEdgeGrid()
{
var inf = float.PositiveInfinity;
var grid = new List<float[]>();
foreach (var a in new[] { -2.5f, 0.0f, 1.5f })
foreach (var b in new[] { -1.5f, -0.0f, 0.0f, 2.0f })
foreach (var lo in new[] { -inf, -1.0f, 0.5f, 3.0f, inf })
foreach (var hi in new[] { -inf, -0.5f, 1.0f, 4.0f, inf })
foreach (var (t0, t1) in new[] { (0.0f, 0.0f), (1.25f, 1.25f), (-1.0f, 2.0f), (2.0f, -1.0f), (0.0f, 1.0e-30f) })
grid.Add(new[] { a, b, lo, hi, t0, t1 });
return grid;
}
/// <summary>
/// The edge grid followed by the committed GLSL fixture's cases of the three integral_clamp_linear families
/// (GlslMirrorTests.IntegralClampLinearFixtureCases, decoded): one list, read by the mirror test and by the
/// clamp-shadow test. The grid fixes a in {-2.5, 0, 1.5}, so a guard keyed on a past 1.75 hides from it
/// (Soul s4); the fixture's random a in [-4, 4] closes that.
/// </summary>
internal static IReadOnlyList<float[]> IntegralClampLinearCases => IntegralClampLinearEdgeGrid.Concat(GlslMirrorTests.IntegralClampLinearFixtureCases).ToArray();
/// <summary>
/// The compiled mirror's cultmath_integral_clamp_linear returns C#'s bits on every edge case (the grid and
/// the fixture's cases), and C#'s result is never NaN there. The generic case set catches each of the function's guards (the held and
/// kept stretch guards, the flat branch's zero-length guard) only by its rotation: with the kept guard
/// removed it finds 0 mismatches, with the held guard removed 1 once infinities join the specials, and the
/// edge grid finds 246 for either every time.
/// </summary>
[Fact]
public void IntegralClampLinearMirrorMatchesCSharpOnTheEdgeCases()
{
var (assembly, errors) = CompileShaderMirror();
Assert.True(assembly is not null, string.Join(Environment.NewLine, errors));
var shaderType = assembly!.GetType("CultMathHlsl.HlslShader")!;
var shader = Activator.CreateInstance(shaderType);
var mirror = ShaderFunction(shaderType, "cultmath_integral_clamp_linear", Enumerable.Repeat(typeof(float), 6).ToArray())!;
var mismatches = new List<string>();
var cases = IntegralClampLinearCases;
foreach (var c in cases)
{
var expected = integral_clamp_linear(c[0], c[1], c[2], c[3], c[4], c[5]);
var actual = (float)mirror.Invoke(shader, c.Cast<object>().ToArray())!;
Assert.False(float.IsNaN(expected), $"C# returns NaN at ({string.Join(", ", c)})");
if (BitConverter.SingleToInt32Bits(expected) != BitConverter.SingleToInt32Bits(actual))
mismatches.Add($"({string.Join(", ", c)}): C# {expected:R}, mirror {actual:R}");
}
Assert.True(mismatches.Count == 0, $"{mismatches.Count} of {cases.Count} cases differ:{Environment.NewLine}{string.Join(Environment.NewLine, mismatches.Take(12))}");
}
/// <summary>
/// At a SpanStage the shader has no identity constant for, so C#'s identity (copied field by field) is the
/// argument: the mirror's cultmath_span_fbm and cultmath_span_ridged return the field of the mirror's own
/// cultmath_fbm_grad and cultmath_ridged_grad, bit for bit, for 0 to 20 octaves over random and simplex-tie
/// points. The shader-side statement of math.SpanStage.Identity (design decision a: no shader identity).
/// </summary>
[Fact]
public void SpanStagesAtIdentityAreFbmAndRidgedInTheMirror()
{
var (assembly, errors) = CompileShaderMirror();
Assert.True(assembly is not null, string.Join(Environment.NewLine, errors));
var shaderType = assembly!.GetType("CultMathHlsl.HlslShader")!;
var shader = Activator.CreateInstance(shaderType);
var (pathType, stageType) = SpanTypes(shaderType);
var spanArguments = new[] { pathType, stageType, typeof(int), typeof(float), typeof(float) };
var pairs = new[]
{
(Span: ShaderFunction(shaderType, "cultmath_span_fbm", spanArguments)!, Plain: ShaderFunction(shaderType, "cultmath_fbm_grad", typeof(float3), typeof(int), typeof(float), typeof(float))!),
(Span: ShaderFunction(shaderType, "cultmath_span_ridged", spanArguments)!, Plain: ShaderFunction(shaderType, "cultmath_ridged_grad", typeof(float3), typeof(int), typeof(float), typeof(float))!),
};
var identity = CopyFields(SpanStage.Identity, stageType);
var random = new System.Random(0x57A6E2);
var compared = 0;
foreach (var p in SpanStageTests.Points(random))
{
var octaves = random.Next(0, 21);
var lacunarity = 1.4f + 1.2f * random.NextSingle();
var gain = 0.3f + 0.5f * random.NextSingle();
var path = CopyFields(new SpanPath { midpoint = p, velocity = new float3(0.6f, 0.0f, 0.8f), speed = 1.0f, acceleration = 0.5f, halfLength = 0.25f }, pathType);
foreach (var (span, plain) in pairs)
{
var field = FieldOf(span.Invoke(shader, new object[] { path, identity, octaves, lacunarity, gain })!, "field");
var expected = plain.Invoke(shader, new object[] { p, octaves, lacunarity, gain })!;
Assert.True(BitwiseEqual(Values(expected), Values(field)), $"{span.Name} at {p}, {octaves} octaves: {string.Join(" ", Values(field).Select(v => v.Value))} is not {plain.Name} {string.Join(" ", Values(expected).Select(v => v.Value))}");
compared++;
}
}
Assert.True(compared >= 460, $"only {compared} identity comparisons ran");
}
/// <summary>
/// The generic comparison draws every argument from [-100, 100] (and specials), where the octave loop's
/// frequency, stretch and span terms saturate: 0 or 16 octaves, a Hessian bound pinned at its global value.
/// This sweep draws SpanStageTests' realistic stages and spans (4000, one in four with a zero axis, one in
/// four with a point span, halfLength 0), both stages, and the mirror must return C#'s bits for the whole
/// SpanNoise. Each of stretched and unstretched stages, and of positive and zero fold, must be met at least
/// 500 times, or the sweep proves nothing about that branch.
/// </summary>
[Fact]
public void SpanStageMirrorMatchesCSharpBitForBitOnRealisticSpans()
{
var (assembly, errors) = CompileShaderMirror();
Assert.True(assembly is not null, string.Join(Environment.NewLine, errors));
var shaderType = assembly!.GetType("CultMathHlsl.HlslShader")!;
var shader = Activator.CreateInstance(shaderType);
var (pathType, stageType) = SpanTypes(shaderType);
var spanArguments = new[] { pathType, stageType, typeof(int), typeof(float), typeof(float) };
var mirrors = new[] { ShaderFunction(shaderType, "cultmath_span_fbm", spanArguments)!, ShaderFunction(shaderType, "cultmath_span_ridged", spanArguments)! };
var random = new System.Random(0x5EED4);
float Next(float extent) => (random.NextSingle() * 2.0f - 1.0f) * extent;
int stretched = 0, unstretched = 0, folded = 0, flat = 0;
for (var i = 0; i < 4000; i++)
{
var drawn = SpanStageTests.Stage.Draw(random);
var stage = drawn.S;
if (i % 4 == 1)
stage.axis = new float3(0.0f, 0.0f, 0.0f);
var halfLength = i % 4 == 2 ? 0.0f : new[] { 0.02f, 0.1f, 0.5f, 2.0f }[random.Next(4)];
var span = SpanStageTests.Span.Draw(random, new float3(Next(30.0f), Next(30.0f), Next(30.0f)), halfLength, random.NextSingle() * 2.0f);
var arguments = new object[] { CopyFields(span.Path, pathType), CopyFields(stage, stageType), drawn.Octaves, drawn.Lacunarity, drawn.Gain };
if (dot(stage.axis, stage.axis) > 0.0f && stage.stretch != 1.0f)
stretched++;
else
unstretched++;
for (var m = 0; m < 2; m++)
{
var expected = m == 0 ? span_fbm(span.Path, stage, drawn.Octaves, drawn.Lacunarity, drawn.Gain) : span_ridged(span.Path, stage, drawn.Octaves, drawn.Lacunarity, drawn.Gain);
var actual = mirrors[m].Invoke(shader, arguments)!;
Assert.True(BitwiseEqual(Values(expected), Values(actual)),
$"{mirrors[m].Name} differs from C# at midpoint {span.Path.midpoint}, halfLength {halfLength:R}, {drawn.Octaves} octaves: C# {string.Join(" ", Values(expected).Select(v => v.Value))}, mirror {string.Join(" ", Values(actual).Select(v => v.Value))}");
if (m == 1 && expected.fold > 0.0f)
folded++;
else if (m == 1)
flat++;
}
}
Assert.True(stretched >= 500 && unstretched >= 500 && folded >= 500 && flat >= 500,
$"the sweep meets a branch too rarely: stretched {stretched}, unstretched {unstretched}, fold > 0 {folded}, fold = 0 {flat} of 4000");
}
/// <summary>
@ -352,14 +544,66 @@ public sealed class HlslSourceCompatibilityTests
}
// Integer arguments take the bit pattern of the float input, so specials become 0, 0x80000000, NaN bits, ...
private static readonly Type[] Scalars = { typeof(float), typeof(int), typeof(uint) };
// A bool argument is true when its float input is positive.
private static readonly Type[] Scalars = { typeof(float), typeof(int), typeof(uint), typeof(bool) };
private static object Scalar(Type t, float f) => t == typeof(float) ? f
: t == typeof(bool) ? f > 0.0f
: t == typeof(int) ? BitConverter.SingleToInt32Bits(f) : (object)(uint)BitConverter.SingleToInt32Bits(f);
// Scoped to scalar-typed fields: this only ever builds function ARGUMENTS (float2/3/4, int2/3/4, ...),
// never a struct-of-struct return, so it must not recurse.
private static FieldInfo[] Components(Type t) => t.GetFields().Where(f => !f.IsStatic && Scalars.Contains(f.FieldType)).ToArray();
private static FieldInfo[] InstanceFields(Type t) => t.GetFields(BindingFlags.Public | BindingFlags.Instance);
// One function ARGUMENT from the case's scalars, in field order: a scalar from its float, a vector by the
// constructor taking its scalar fields, a struct field by field, recursively.
private static object Build(Type t, Func<float> next)
{
if (Scalars.Contains(t))
return Scalar(t, next());
var fields = InstanceFields(t);
var constructor = t.GetConstructor(fields.Select(f => f.FieldType).ToArray());
if (constructor is not null)
return constructor.Invoke(fields.Select(f => Build(f.FieldType, next)).ToArray());
var value = Activator.CreateInstance(t)!;
foreach (var field in fields)
field.SetValue(value, Build(field.FieldType, next));
return value;
}
// The CultMath type a mirror parameter type stands for: itself, or for a type nested in HlslShader the
// CultMath type of the same name, only if both have the same instance fields by name, type (mapped
// likewise) and order. Null when there is no such twin.
private static Type? TwinType(Type shaderType, Type mirror)
{
if (mirror.DeclaringType != shaderType)
return mirror;
var twin = typeof(math).Assembly.GetType("CultMath." + mirror.Name);
if (twin is null)
return null;
var mine = InstanceFields(mirror);
var theirs = InstanceFields(twin);
return mine.Length == theirs.Length && mine.Zip(theirs).All(p => p.First.Name == p.Second.Name && TwinType(shaderType, p.First.FieldType) == p.Second.FieldType) ? twin : null;
}
// The math function of the mirror function's name, with its parameter types mapped, public or private.
private static MethodInfo? Twin(Type shaderType, string name, Type[] mirrorTypes)
{
var types = mirrorTypes.Select(t => TwinType(shaderType, t)).ToArray();
return types.Any(t => t is null) ? null
: typeof(math).GetMethod(name, BindingFlags.Public | BindingFlags.NonPublic | BindingFlags.Static, null, types!, null);
}
// A struct value as another struct type of the same fields by name: itself when the type is the same.
private static object CopyFields(object source, Type target)
{
if (source.GetType() == target)
return source;
var copy = Activator.CreateInstance(target)!;
foreach (var field in InstanceFields(target))
field.SetValue(copy, CopyFields(source.GetType().GetField(field.Name)!.GetValue(source)!, field.FieldType));
return copy;
}
private static object FieldOf(object value, string name) => value.GetType().GetField(name)!.GetValue(value)!;
/// <summary>
/// One flattened scalar leaf of a (possibly nested) return value: its declaration path (field
@ -412,8 +656,11 @@ public sealed class HlslSourceCompatibilityTests
: string.Empty);
}
/// <summary>The documented HLSL-to-C# transformations, and nothing else.</summary>
internal static string TransformHlslToCSharp(string hlsl)
/// <summary>
/// The documented HLSL-to-C# transformations, and nothing else. <paramref name="members"/> is C# source a
/// test supplies, placed first in the class body: it is not a transformation of the HLSL.
/// </summary>
internal static string TransformHlslToCSharp(string hlsl, string members = "")
{
// 1. Include guards (#ifndef/#define/#endif) are dropped; C# has no textual include.
var source = Regex.Replace(hlsl, @"(?m)^\s*#(ifndef|define|endif)\b.*$", string.Empty);
@ -430,18 +677,19 @@ public sealed class HlslSourceCompatibilityTests
// 5. Struct fields gain `public`; HLSL has no field access-modifier concept, and a C# struct's
// fields default to private, which would hide them from the mirror test's reflection-based
// field comparison (invariant 8's one struct return shape, e.g. CultCellular).
// field comparison (invariant 8's one struct return shape, e.g. CultCellular) and the struct
// arguments' copies (SpanPath, SpanStage), which reflection reads and writes by field name.
source = Regex.Replace(source, @"(?ms)(struct\s+\w+\s*\{)(.*?)(\};)", m =>
m.Groups[1].Value + Regex.Replace(m.Groups[2].Value, @"(?m)^(\s*)(\w+\s+\w+;)", "$1public $2") + m.Groups[3].Value);
// 6. The file body is wrapped in a class; C# has no free functions.
return "using CultMath;\nusing static CultMath.math;\nnamespace CultMathHlsl;\npublic class HlslShader\n{\n" + source + "\n}\n";
return "using CultMath;\nusing static CultMath.math;\nnamespace CultMathHlsl;\npublic class HlslShader\n{\n" + members + source + "\n}\n";
}
internal static (Assembly? Assembly, IReadOnlyList<string> Errors) CompileShaderMirror()
internal static (Assembly? Assembly, IReadOnlyList<string> Errors) CompileShaderMirror(string members = "")
{
var hlsl = ReadShaderSource(FindCultMathRoot());
var tree = CSharpSyntaxTree.ParseText(TransformHlslToCSharp(hlsl));
var tree = CSharpSyntaxTree.ParseText(TransformHlslToCSharp(hlsl, members));
var references = ((string)AppContext.GetData("TRUSTED_PLATFORM_ASSEMBLIES")!)
.Split(Path.PathSeparator)
.Append(typeof(math).Assembly.Location)

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

@ -469,8 +469,8 @@ public sealed class NoiseBoundTests
/// <summary>
/// Pins SNOISE_LIPSCHITZ: no |snoise_grad| over 1e5 seeded points exceeds it, and the ascent from
/// MeasureLipschitz's witness reaches a maximum that, times 1.10, is the constant to within half a
/// percent, so a constant lowered (or raised) by 1% fails here.
/// MeasureLipschitz's witness reaches a maximum that, times 1.10, equals the constant bit for bit
/// (the expression MeasureLipschitz prints), so any moved constant fails here.
/// </summary>
[Fact]
public void LipschitzConstantPinsSampledGradients()
@ -479,8 +479,8 @@ public sealed class NoiseBoundTests
Assert.True(sampled <= SNOISE_LIPSCHITZ, $"|snoise_grad| reaches {sampled:R}, above SNOISE_LIPSCHITZ {SNOISE_LIPSCHITZ:R}");
var refined = AscendGradientNorm(LipschitzWitness);
Assert.True(refined <= SNOISE_LIPSCHITZ, $"|snoise_grad| reaches {refined:R}, above SNOISE_LIPSCHITZ {SNOISE_LIPSCHITZ:R}");
var ratio = SNOISE_LIPSCHITZ / (1.10f * refined);
Assert.True(MathF.Abs(ratio - 1.0f) <= 0.005f, $"SNOISE_LIPSCHITZ is {ratio:R} times 1.10 x the re-measured {refined:R}; re-run MeasureLipschitz and re-pin");
var measured = 1.10f * refined;
Assert.True(SNOISE_LIPSCHITZ == measured, $"SNOISE_LIPSCHITZ is {SNOISE_LIPSCHITZ:R}, the re-measured 1.10 x {refined:R} is {measured:R}; re-run MeasureLipschitz and re-pin");
}
// ---- The curvature and 2D Lipschitz constants (math.Affine.cs, math.Interval.cs) ----
@ -609,8 +609,8 @@ public sealed class NoiseBoundTests
/// <summary>
/// Pins SNOISE_HESSIAN: no Hessian norm over 1e5 seeded points exceeds it, and the ascent from
/// MeasureHessian'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.
/// MeasureHessian's witness reaches a maximum that, times 1.10, equals the constant bit for bit
/// (the expression MeasureHessian prints), so any moved constant fails here.
/// </summary>
[Fact]
public void HessianConstantPinsSampledCurvature()
@ -619,39 +619,202 @@ public sealed class NoiseBoundTests
Assert.True(sampled <= SNOISE_HESSIAN, $"the Hessian norm reaches {sampled:R}, above SNOISE_HESSIAN {SNOISE_HESSIAN:R}");
var refined = Ascend(HessianNorm, HessianWitness);
Assert.True(refined <= SNOISE_HESSIAN, $"the Hessian norm reaches {refined:R}, above SNOISE_HESSIAN {SNOISE_HESSIAN:R}");
var ratio = SNOISE_HESSIAN / (1.10 * refined);
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");
var measured = (float)(refined * 1.10);
Assert.True(SNOISE_HESSIAN == measured, $"SNOISE_HESSIAN is {SNOISE_HESSIAN:R}, the re-measured 1.10 x {refined:R} is {measured: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)
private static float3 CultMod289(float3 v) => v - floor(v * (1.0f / 289.0f)) * 289.0f;
private static float4 CultMod289(float4 v) => v - floor(v * (1.0f / 289.0f)) * 289.0f;
private static float4 CultPermute(float4 v) => CultMod289(((v * 34.0f) + 10.0f) * v);
// The four simplex corners of snoise at (vx, vy, vz): each corner's offset x_c from the point, and its
// gradient direction p_c, in double, with snoise_hessian's lattice, step() tie order and (float) hash,
// so one corner table serves the Hessian and its derivative. Key names the lattice cell and traversal
// order: two points with the same key are differentiated through the same corners.
internal static (double[][] X, double[][] P, double[] Key) Corners(double vx, double vy, double vz)
{
const float h = 1.0e-3f;
var d = new double[27];
for (var k = 0; k < 3; k++)
const double cx = 1.0 / 6.0;
const double cy = 1.0 / 3.0;
var skew = (vx + vy + vz) * cy;
double ix = Math.Floor(vx + skew), iy = Math.Floor(vy + skew), iz = Math.Floor(vz + skew);
var unskew = (ix + iy + iz) * cx;
double x0 = vx - ix + unskew, y0 = vy - iy + unskew, z0 = vz - iz + unskew;
double gx = x0 >= y0 ? 1.0 : 0.0, gy = y0 >= z0 ? 1.0 : 0.0, gz = 1.0 - (x0 >= z0 ? 1.0 : 0.0);
double lx = 1.0 - gx, ly = 1.0 - gy, lz = 1.0 - gz;
double i1x = Math.Min(gx, lz), i1y = Math.Min(gy, lx), i1z = Math.Min(gz, ly);
double i2x = Math.Max(gx, lz), i2y = Math.Max(gy, lx), i2z = Math.Max(gz, ly);
var x = new[]
{
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);
new[] { x0, y0, z0 },
new[] { x0 - i1x + cx, y0 - i1y + cx, z0 - i1z + cx },
new[] { x0 - i2x + cy, y0 - i2y + cy, z0 - i2z + cy },
new[] { x0 - 0.5, y0 - 0.5, z0 - 0.5 },
};
var i = CultMod289(new float3((float)ix, (float)iy, (float)iz));
var hash = CultPermute(CultPermute(CultPermute(
i.z + new float4(0.0f, (float)i1z, (float)i2z, 1.0f)) + i.y + new float4(0.0f, (float)i1y, (float)i2y, 1.0f)) + i.x + new float4(0.0f, (float)i1x, (float)i2x, 1.0f));
const float n = 0.142857142857f;
var ns = new float3(2.0f * n, 0.5f * n - 1.0f, n);
var j = hash - 49.0f * floor(hash * ns.z * ns.z);
var xs = floor(j * ns.z);
var ys = floor(j - 7.0f * xs);
var gxs = xs * ns.x + ns.y;
var gys = ys * ns.x + ns.y;
var h = 1.0f - abs(gxs) - abs(gys);
var b0 = new float4(gxs.x, gxs.y, gys.x, gys.y);
var b1 = new float4(gxs.z, gxs.w, gys.z, gys.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;
double[] Of(float3 v) => new double[] { v.x, v.y, v.z };
return (x, new[] { Of(p0), Of(p1), Of(p2), Of(p3) }, new[] { ix, iy, iz, i1x, i1y, i1z, i2x, i2y, i2z });
}
// The Hessian of snoise at v in double, entry i * 3 + j: snoise_hessian's per-corner derivation, summed
// over Corners. A corner with m = 0.5 - |x|^2 <= 0 contributes nothing.
internal static double[] HessianDouble(double vx, double vy, double vz)
{
var (x, p, _) = Corners(vx, vy, vz);
var hessian = new double[9];
for (var c = 0; c < 4; c++)
{
var m = 0.5 - (x[c][0] * x[c][0] + x[c][1] * x[c][1] + x[c][2] * x[c][2]);
if (m <= 0.0)
continue;
var s = p[c][0] * x[c][0] + p[c][1] * x[c][1] + p[c][2] * x[c][2];
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);
for (var k = 0; k < 3; k++)
hessian[i * 3 + k] += 105.0 * (48.0 * m * m * s * x[c][i] * x[c][k]
- 8.0 * m * m * m * (p[c][k] * x[c][i] + p[c][i] * x[c][k] + (i == k ? s : 0.0)));
}
return hessian;
}
// The value of snoise at v in double: 105 sum m^4 s over Corners, the kernel's exact value (the tied
// kernel, as HessianDouble is its exact Hessian).
internal static double ValueDouble(double vx, double vy, double vz)
{
var (x, p, _) = Corners(vx, vy, vz);
var total = 0.0;
for (var c = 0; c < 4; c++)
{
var m = 0.5 - (x[c][0] * x[c][0] + x[c][1] * x[c][1] + x[c][2] * x[c][2]);
if (m <= 0.0)
continue;
total += 105.0 * m * m * m * m * (p[c][0] * x[c][0] + p[c][1] * x[c][1] + p[c][2] * x[c][2]);
}
return total;
}
private static double Delta(int a, int b) => a == b ? 1.0 : 0.0;
// The symmetric third-derivative tensor of snoise at p, entry (i * 3 + j) * 3 + k, in closed form:
// snoise_hessian's per-corner entry 105 (48 m^2 s x_i x_j - 8 m^3 (p_j x_i + p_i x_j + s delta_ij))
// differentiated by x_k (dm = -2 x_k, ds = p_k, every corner's x moves one for one with p), which sums
// to 105 (-192 m s x_i x_j x_k + 48 m^2 (p_k x_i x_j + p_j x_i x_k + p_i x_j x_k)
// + 48 m^2 s (delta_ik x_j + delta_jk x_i + delta_ij x_k) - 8 m^3 (p_j delta_ik + p_i delta_jk + p_k delta_ij)).
// Symmetric in i, j and k by construction, so there is no symmetrising pass.
private static double[] ThirdTensor(float3 p)
{
var (x, g, _) = Corners(p.x, p.y, p.z);
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;
for (var c = 0; c < 4; c++)
{
var m = 0.5 - (x[c][0] * x[c][0] + x[c][1] * x[c][1] + x[c][2] * x[c][2]);
if (m <= 0.0)
continue;
var s = g[c][0] * x[c][0] + g[c][1] * x[c][1] + g[c][2] * x[c][2];
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] += 105.0 * (-192.0 * m * s * x[c][i] * x[c][j] * x[c][k]
+ 48.0 * m * m * (g[c][k] * x[c][i] * x[c][j] + g[c][j] * x[c][i] * x[c][k] + g[c][i] * x[c][j] * x[c][k])
+ 48.0 * m * m * s * (Delta(i, k) * x[c][j] + Delta(j, k) * x[c][i] + Delta(i, j) * x[c][k])
- 8.0 * m * m * m * (g[c][j] * Delta(i, k) + g[c][i] * Delta(j, k) + g[c][k] * Delta(i, j)));
}
return t;
}
/// <summary>
/// ThirdTensor is the third derivative of snoise: HessianDouble is the kernel's Hessian (within float
/// rounding of snoise_hessian over 200 seeded points, which proves Corners is the kernel's lattice),
/// and the closed form equals central differences of HessianDouble, Richardson-extrapolated from steps
/// h = 1e-3 and 2e-3 (error O(h^4)), at every one of 27 entries of every point whose neighbours share
/// its lattice cell and clear every corner's support edge. Dropping or mis-signing any term of the closed form moves some entry by far more
/// than the tolerance.
/// </summary>
[Fact]
public void ThirdTensorMatchesDifferencesOfTheHessian()
{
var random = new System.Random(0x7D1);
var compared = 0;
var largest = 0.0;
var worst = 0.0;
for (var n = 0; n < 400; n++)
{
var p = new float3(Uniform(random, -256.0f, 256.0f), Uniform(random, -256.0f, 256.0f), Uniform(random, -256.0f, 256.0f));
var hd = HessianDouble(p.x, p.y, p.z);
var kernel = snoise_hessian(p);
if (n < 200)
for (var i = 0; i < 3; i++)
for (var j = 0; j < 3; j++)
Assert.True(Math.Abs(hd[i * 3 + j] - kernel[i][j]) <= 0.1, $"HessianDouble({p}) entry {i},{j} is {hd[i * 3 + j]:R}, snoise_hessian {kernel[i][j]:R}");
double[] At(int axis, double delta) =>
HessianDouble(p.x + (axis == 0 ? delta : 0.0), p.y + (axis == 1 ? delta : 0.0), p.z + (axis == 2 ? delta : 0.0));
var (corner, _, key) = Corners(p.x, p.y, p.z);
// A corner whose support edge (m = 0) lies within the differencing window has a kink in the
// Hessian's second derivative there, which the O(h^4) extrapolation does not cover.
if (corner.Any(c => Math.Abs(0.5 - (c[0] * c[0] + c[1] * c[1] + c[2] * c[2])) < 0.01))
continue;
var sameCell = true;
for (var axis = 0; axis < 3; axis++)
foreach (var delta in new[] { -2.0e-3, 2.0e-3 })
sameCell &= Corners(p.x + (axis == 0 ? delta : 0.0), p.y + (axis == 1 ? delta : 0.0), p.z + (axis == 2 ? delta : 0.0)).Key.SequenceEqual(key);
if (!sameCell)
continue;
var closed = ThirdTensor(p);
compared++;
for (var k = 0; k < 3; k++)
{
var (plus1, minus1, plus2, minus2) = (At(k, 1.0e-3), At(k, -1.0e-3), At(k, 2.0e-3), At(k, -2.0e-3));
for (var i = 0; i < 3; i++)
for (var j = 0; j < 3; j++)
{
var d1 = (plus1[i * 3 + j] - minus1[i * 3 + j]) / 2.0e-3;
var d2 = (plus2[i * 3 + j] - minus2[i * 3 + j]) / 4.0e-3;
var richardson = (4.0 * d1 - d2) / 3.0;
var entry = closed[(i * 3 + j) * 3 + k];
largest = Math.Max(largest, Math.Abs(entry));
worst = Math.Max(worst, Math.Abs(entry - richardson) / (1.0 + Math.Abs(entry)));
}
}
}
output.WriteLine($"IV-REPORT third oracle: {compared} points, largest entry {largest:R}, worst relative difference from the Richardson differences {worst:R}");
Assert.True(compared >= 250, $"only {compared} of 400 points kept their lattice cell and cleared the support edges");
Assert.True(largest >= 200.0, $"the largest third-derivative entry compared is {largest:R}: the check has no teeth");
Assert.True(worst <= 1.0e-6, $"the closed-form third tensor differs from the differenced Hessian by {worst:R} (relative)");
}
// T(u, u, u) for a symmetric 3-tensor.
private static double Cubic(double[] t, double[] u)
{
@ -748,8 +911,8 @@ public sealed class NoiseBoundTests
}
/// <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,
/// The provenance of SNOISE_THIRD: the largest operator norm of snoise's closed-form third-derivative
/// tensor (ThirdNorm, ThirdTensor) 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
@ -780,12 +943,12 @@ public sealed class NoiseBoundTests
}
// The start MeasureThird printed as its witness.
private static readonly float3 ThirdWitness = new(-168.3755f, -3.8205261f, -109.16948f);
private static readonly float3 ThirdWitness = new(-184.81363f, -34.250977f, 68.26865f);
/// <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.
/// MeasureThird's witness reaches a maximum that, times 1.10, equals the constant bit for bit
/// (the expression MeasureThird prints), so any moved constant fails here.
/// </summary>
[Fact]
public void ThirdConstantPinsSampled()
@ -794,8 +957,8 @@ public sealed class NoiseBoundTests
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");
var measured = (float)(refined * 1.10);
Assert.True(SNOISE_THIRD == measured, $"SNOISE_THIRD is {SNOISE_THIRD:R}, the re-measured 1.10 x {refined:R} is {measured:R}; re-run MeasureThird and re-pin");
}
/// <summary>
@ -856,7 +1019,7 @@ public sealed class NoiseBoundTests
/// <summary>
/// Pins SNOISE2_LIPSCHITZ as LipschitzConstantPinsSampledGradients pins SNOISE_LIPSCHITZ: 1e5 seeded
/// points under it, and the ascent from MeasureLipschitz2's witness within half a percent of it / 1.10.
/// points under it, and the ascent from MeasureLipschitz2's witness, times 1.10, equals it bit for bit (the expression MeasureLipschitz2 prints).
/// </summary>
[Fact]
public void Lipschitz2ConstantPins()
@ -865,8 +1028,8 @@ public sealed class NoiseBoundTests
Assert.True(sampled <= SNOISE2_LIPSCHITZ, $"|grad snoise2| reaches {sampled:R}, above SNOISE2_LIPSCHITZ {SNOISE2_LIPSCHITZ:R}");
var refined = Ascend2(Gradient2Norm, Lipschitz2Witness);
Assert.True(refined <= SNOISE2_LIPSCHITZ, $"|grad snoise2| reaches {refined:R}, above SNOISE2_LIPSCHITZ {SNOISE2_LIPSCHITZ:R}");
var ratio = SNOISE2_LIPSCHITZ / (1.10 * refined);
Assert.True(Math.Abs(ratio - 1.0) <= 0.005, $"SNOISE2_LIPSCHITZ is {ratio:R} times 1.10 x the re-measured {refined:R}; re-run MeasureLipschitz2 and re-pin");
var measured = (float)(refined * 1.10);
Assert.True(SNOISE2_LIPSCHITZ == measured, $"SNOISE2_LIPSCHITZ is {SNOISE2_LIPSCHITZ:R}, the re-measured 1.10 x {refined:R} is {measured:R}; re-run MeasureLipschitz2 and re-pin");
}
/// <summary>

57
packages/cultmath/tests/CultMath.Tests/SpanPrimitiveTests.cs
View file

@ -319,6 +319,63 @@ public sealed class SpanPrimitiveTests
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");
}
/// <summary>
/// An infinite bound never reached by the line contributes nothing, not infinity times zero: relu
/// (hi = +inf) and its mirror (lo = -inf) give the finite integral, whether the line crosses the finite
/// bound inside the span or not. Both branches give +0 for a zero-length span.
/// </summary>
[Fact]
public void IntegralClampLinearIsFiniteForInfiniteBounds()
{
var inf = float.PositiveInfinity;
// t over [0, 1] clamped below at 0: all ramp, 1/2.
Assert.Equal(0.5f, integral_clamp_linear(0.0f, 1.0f, 0.0f, inf, 0.0f, 1.0f));
// t - 1/2 over [0, 1] clamped below at 0: 0 until 1/2, then the ramp's 1/8.
Assert.Equal(0.125f, integral_clamp_linear(-0.5f, 1.0f, 0.0f, inf, 0.0f, 1.0f));
// The falling line 1/2 - t clamped below at 0 is the mirror: 1/8 over [0, 1].
Assert.Equal(0.125f, integral_clamp_linear(0.5f, -1.0f, 0.0f, inf, 0.0f, 1.0f));
// lo = -inf: clamp(t, -inf, 1/2) over [0, 1] is t up to 1/2 then 1/2 (1/8 + 1/4), and 1 - t clamped
// above at 1/2 is 1/2 until t = 1/2 and then falls (1/4 + 1/8).
Assert.Equal(0.375f, integral_clamp_linear(0.0f, 1.0f, -inf, 0.5f, 0.0f, 1.0f));
Assert.Equal(0.375f, integral_clamp_linear(1.0f, -1.0f, -inf, 0.5f, 0.0f, 1.0f));
// Both infinite: the line itself, the trapezoid, in either direction.
Assert.Equal(7.5f, integral_clamp_linear(0.0f, 1.0f, -inf, inf, 1.0f, 4.0f));
Assert.Equal(-7.5f, integral_clamp_linear(0.0f, 1.0f, -inf, inf, 4.0f, 1.0f));
// A line starting on the finite bound and rising into the infinite one: all ramp, the integral of t over [2, 4].
Assert.Equal(6.0f, integral_clamp_linear(0.0f, 1.0f, 2.0f, inf, 2.0f, 4.0f));
// The flat branch obeys the same zero-length rule: a stretch of zero length contributes +0 however
// infinite the bound it would be clamped to, and a nonzero length of an infinite constant is infinite.
// The bit patterns are compared because Assert.Equal treats -0 as 0.
Assert.Equal(0, BitConverter.SingleToInt32Bits(integral_clamp_linear(-2.3244874f, -0.0f, 0.0f, -inf, 1.429815f, 1.429815f)));
Assert.Equal(0, BitConverter.SingleToInt32Bits(integral_clamp_linear(1.5f, 0.0f, inf, inf, 2.0f, 2.0f)));
Assert.Equal(0, BitConverter.SingleToInt32Bits(integral_clamp_linear(1.5f, -0.0f, -inf, -inf, 0.0f, 0.0f)));
Assert.Equal(inf, integral_clamp_linear(1.5f, 0.0f, inf, inf, 0.0f, 2.0f));
Assert.Equal(inf, integral_clamp_linear(1.5f, 0.0f, -inf, -inf, 2.0f, 0.0f));
// The sloped branch gives +0 for a zero-length span too, not -0: the direction sign never meets a
// zero total, whichever way the line leans and whatever the bounds.
Assert.Equal(0, BitConverter.SingleToInt32Bits(integral_clamp_linear(1.0f, 2.0f, -inf, inf, 1.0f, 1.0f)));
Assert.Equal(0, BitConverter.SingleToInt32Bits(integral_clamp_linear(-3.0f, -2.0f, 0.0f, 5.0f, 1.0f, 1.0f)));
Assert.Equal(0, BitConverter.SingleToInt32Bits(integral_clamp_linear(0.5f, 1.0f, 0.0f, inf, 0.0f, 0.0f)));
}
/// <summary>
/// lo above hi is normalised before the flat branch reads it: with b = 0 the answer is hi times the
/// span (what clamp(a, lo, hi) gives with lo = min(lo, hi)), and with b != 0 the same. The GLSL
/// lowering would otherwise call clamp with minVal above maxVal, which is undefined there.
/// </summary>
[Fact]
public void IntegralClampLinearNormalisesBeforeTheFlatBranch()
{
// b = 0: a above both, below both and between them all give hi.
foreach (var a in new[] { 9.0f, -9.0f, 3.5f })
Assert.Equal(2.0f * 3.0f, integral_clamp_linear(a, 0.0f, 5.0f, 2.0f, 0.0f, 3.0f));
// b != 0: hi everywhere, whichever way the line runs.
Assert.Equal(2.0f * 3.0f, integral_clamp_linear(0.0f, 1.0f, 5.0f, 2.0f, 0.0f, 3.0f));
Assert.Equal(2.0f * 3.0f, integral_clamp_linear(9.0f, -1.0f, 5.0f, 2.0f, 0.0f, 3.0f));
}
[Fact]
public void IntegralClampLinearHasExactValuesOnTheBinaryFloats()
{

893
packages/cultmath/tests/CultMath.Tests/SpanStageTests.cs Normal file
View file

@ -0,0 +1,893 @@
using Xunit;
using static CultMath.math;
namespace CultMath.Tests;
/// <summary>
/// span_fbm and span_ridged (math.Span.cs): the octave stages of a warp chain over a span. At the identity
/// transform they are fbm_grad and ridged_grad; their gradient is the derivative of their value through the
/// stretch and the per-octave scaling; curvature bounds the Hessian of the field in q and fold the departure
/// of a ridged sum from its midpoint branch, at every point of every span the tests draw, in the regime
/// where the local third-derivative term is what holds the bound and in the one where the global pin is.
/// Lines starting IV-REPORT carry the numbers the cut report quotes.
/// </summary>
public sealed class SpanStageTests
{
private readonly ITestOutputHelper output;
public SpanStageTests(ITestOutputHelper output) => this.output = output;
private static float Uniform(System.Random random, float lo, float hi) => lo + (hi - lo) * random.NextSingle();
private static float3 Vector(System.Random random, float lo, float hi) => new(Uniform(random, lo, hi), Uniform(random, lo, hi), Uniform(random, lo, hi));
private static float3 UnitVector(System.Random random)
{
float3 v;
do
{
v = Vector(random, -1.0f, 1.0f);
}
while (dot(v, v) is > 1.0f or < 0.01f);
return normalize(v);
}
private static float LogUniform(System.Random random, float lo, float hi) => MathF.Exp(Uniform(random, MathF.Log(lo), MathF.Log(hi)));
private static float Signed(System.Random random, float magnitude) => random.Next(2) == 0 ? magnitude : -magnitude;
private static readonly float3 None = new(0.0f, 0.0f, 0.0f);
private static double L1(double[] v) => Math.Abs(v[0]) + Math.Abs(v[1]) + Math.Abs(v[2]);
// One draw of the stage parameters: every one of them off its identity value, over the ranges the stage
// takes in use (stage-parameter-range): |frequency| 1e-3..10, stretch 1e-2..1e2, |lacunarity| 0.5..3,
// |gain| 0.1..1.2, the last three and frequency with random signs where a sign is meaningful,
// octaveStep nonzero in every component.
internal readonly record struct Stage(SpanStage S, int Octaves, float Lacunarity, float Gain)
{
public static Stage Draw(System.Random random, int maxOctaves = 8, int minOctaves = 1) => new(
new SpanStage
{
frequency = Signed(random, LogUniform(random, 1.0e-3f, 10.0f)),
stretch = LogUniform(random, 1.0e-2f, 1.0e2f),
axis = UnitVector(random) * Uniform(random, 0.5f, 3.0f),
shift = Vector(random, -8.0f, 8.0f),
offset = Vector(random, -8.0f, 8.0f),
octaveStep = new float3(Signed(random, Uniform(random, 0.1f, 3.0f)), Signed(random, Uniform(random, 0.1f, 3.0f)), Signed(random, Uniform(random, 0.1f, 3.0f))),
},
random.Next(minOctaves, maxOctaves + 1),
Signed(random, Uniform(random, 0.5f, 3.0f)),
Signed(random, Uniform(random, 0.1f, 1.2f)));
// A draw whose octave coordinates stay under `limit` for every q with |q|_1 <= qL1: the fastest
// octave's |f_o| |M| is capped by rejection. For the derivative checks, where float rounding of the
// coordinate must not swamp a difference quotient.
public static Stage DrawSmooth(System.Random random, float limit, float qL1, int maxOctaves = 8, int minOctaves = 1)
{
while (true)
{
var stage = Draw(random, maxOctaves, minOctaves);
var worst = 0.0;
for (var o = 0; o < stage.Octaves; o++)
worst = Math.Max(worst, Math.Abs(stage.Scale(o)) * (stage.Norm() * qL1 + L1(Components(stage.S.shift))) + L1(Components(stage.S.offset)) + o * L1(Components(stage.S.octaveStep)));
if (worst < limit)
return stage;
}
}
public SpanNoise Fbm(SpanPath path) => span_fbm(path, S, Octaves, Lacunarity, Gain);
public SpanNoise Ridged(SpanPath path) => span_ridged(path, S, Octaves, Lacunarity, Gain);
public SpanNoise Fbm(SpanPath path, int octaves) => span_fbm(path, S, octaves, Lacunarity, Gain);
public SpanNoise Ridged(SpanPath path, int octaves) => span_ridged(path, S, octaves, Lacunarity, Gain);
// M = frequency (I + (1/stretch - 1) a a^T) in double, a the normalised axis (none for a zero axis).
public double[,] Matrix()
{
var m = new double[3, 3];
var a = Components(S.axis);
var length = Math.Sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]);
for (var i = 0; i < 3; i++)
for (var j = 0; j < 3; j++)
m[i, j] = S.frequency * ((i == j ? 1.0 : 0.0) + (length > 0.0 && S.stretch != 1.0f ? (1.0 / S.stretch - 1.0) * a[i] * a[j] / (length * length) : 0.0));
return m;
}
// Octave o's noise-space point for the span point q, in double from the same inputs.
public double[] PointD(float3 q, int octave)
{
var m = Matrix();
var scale = Math.Pow(Lacunarity, octave);
var x = new double[3];
var qv = Components(q);
var shift = Components(S.shift);
var offset = Components(S.offset);
var step = Components(S.octaveStep);
for (var i = 0; i < 3; i++)
{
var mq = 0.0;
for (var j = 0; j < 3; j++)
mq += m[i, j] * qv[j];
x[i] = scale * (mq + shift[i]) + offset[i] + octave * step[i];
}
return x;
}
public float3 Point(float3 q, int octave)
{
var x = PointD(q, octave);
return new float3((float)x[0], (float)x[1], (float)x[2]);
}
// S_o: the magnitudes octave o's coordinate roundings are relative to, plus the span's reach.
public double Magnitude(SpanPath path, int octave)
{
var m = Matrix();
var qv = Components(path.midpoint);
var mq = new double[3];
for (var i = 0; i < 3; i++)
for (var j = 0; j < 3; j++)
mq[i] += m[i, j] * qv[j];
var scale = Math.Abs(Scale(octave));
return scale * (L1(mq) + L1(Components(S.shift))) + L1(Components(S.offset)) + octave * L1(Components(S.octaveStep))
+ Math.Sqrt(3.0) * scale * Norm() * path.speed * path.halfLength;
}
public double MaxMagnitude(SpanPath path)
{
var largest = 0.0;
for (var o = 0; o < Octaves; o++)
largest = Math.Max(largest, Magnitude(path, o));
return largest;
}
// Octave o as a stage of its own: one octave at the scaled frequency, shift and offset.
public Stage OneOctave(int octave) => new(
new SpanStage { frequency = S.frequency * (float)Scale(octave), stretch = S.stretch, axis = S.axis, shift = S.shift * (float)Scale(octave), offset = S.offset + (float)octave * S.octaveStep, octaveStep = None },
1, 1.0f, 1.0f);
public double Amplitude(int octave) => Math.Pow(Gain, octave);
public double Scale(int octave) => Math.Pow(Lacunarity, octave);
// The spectral norm of M: frequency times the larger of 1 and 1/stretch.
public double Norm() => Math.Abs(S.frequency) * Math.Max(1.0, 1.0 / S.stretch);
}
private static double[] Components(float3 v) => new double[] { v.x, v.y, v.z };
private static SpanPath PathAt(float3 q) => new() { midpoint = q, velocity = None, speed = 0.0f, acceleration = 0.0f, halfLength = 0.0f };
// A span through q: velocity v, constant acceleration c of norm `acceleration`, so q(t) = q + v t + c t^2 / 2;
// the speed bound is |v| plus what c can add over halfLength.
internal readonly record struct Span(SpanPath Path, float3 Acceleration)
{
public static Span Draw(System.Random random, float3 q, float halfLength, float acceleration)
{
var velocity = UnitVector(random) * Uniform(random, 0.2f, 3.0f);
return new Span(new SpanPath { midpoint = q, velocity = velocity, speed = length(velocity) + acceleration * halfLength, acceleration = acceleration, halfLength = halfLength }, UnitVector(random) * acceleration);
}
public float3 At(float t) => Path.midpoint + Path.velocity * t + Acceleration * (0.5f * t * t);
}
// The operator norm of a symmetric 3x3 matrix (entries in row-major order): the larger |eigenvalue|.
private static double SymmetricNorm(double[] a)
{
var p1 = a[1] * a[1] + a[2] * a[2] + a[5] * a[5];
var q = (a[0] + a[4] + a[8]) / 3.0;
if (p1 == 0.0)
return Math.Max(Math.Abs(a[0]), Math.Max(Math.Abs(a[4]), Math.Abs(a[8])));
var p2 = (a[0] - q) * (a[0] - q) + (a[4] - q) * (a[4] - q) + (a[8] - q) * (a[8] - q) + 2.0 * p1;
var s = Math.Sqrt(p2 / 6.0);
var b = new[] { (a[0] - q) / s, a[1] / s, a[2] / s, a[3] / s, (a[4] - q) / s, a[5] / s, a[6] / s, a[7] / s, (a[8] - q) / s };
var r = Math.Clamp((b[0] * (b[4] * b[8] - b[5] * b[7]) - b[1] * (b[3] * b[8] - b[5] * b[6]) + b[2] * (b[3] * b[7] - b[4] * b[6])) / 2.0, -1.0, 1.0);
var phi = Math.Acos(r) / 3.0;
return Math.Max(Math.Abs(q + 2.0 * s * Math.Cos(phi)), Math.Abs(q + 2.0 * s * Math.Cos(phi + 2.0 * Math.PI / 3.0)));
}
// The exact Hessian in q of the stage field at span point q (signed per octave for a ridged sum), by the
// chain rule from the double kernel's Hessian: sum over octaves of weight_o f_o^2 M H_o M, with every
// octave point computed in double from the same SpanStage.
private static double[] FieldHessian(Stage stage, float3 q, bool ridged)
{
var m = stage.Matrix();
var total = new double[9];
for (var o = 0; o < stage.Octaves; o++)
{
var x = stage.PointD(q, o);
var h = NoiseBoundTests.HessianDouble(x[0], x[1], x[2]);
var weight = stage.Amplitude(o) * stage.Scale(o) * stage.Scale(o) * (ridged ? -Math.Sign(NoiseBoundTests.ValueDouble(x[0], x[1], x[2])) : 1.0);
for (var i = 0; i < 3; i++)
for (var j = 0; j < 3; j++)
{
var sum = 0.0;
for (var k = 0; k < 3; k++)
for (var l = 0; l < 3; l++)
sum += m[i, k] * h[k * 3 + l] * m[l, j];
total[i * 3 + j] += weight * sum;
}
}
return total;
}
private static bool SameBits(float4 a, float4 b) =>
BitConverter.SingleToInt32Bits(a.x) == BitConverter.SingleToInt32Bits(b.x) && BitConverter.SingleToInt32Bits(a.y) == BitConverter.SingleToInt32Bits(b.y)
&& BitConverter.SingleToInt32Bits(a.z) == BitConverter.SingleToInt32Bits(b.z) && BitConverter.SingleToInt32Bits(a.w) == BitConverter.SingleToInt32Bits(b.w);
internal static IEnumerable<float3> Points(System.Random random)
{
for (var k = 0; k < 200; k++)
yield return Vector(random, -60.0f, 60.0f);
// The simplex tie set: x = y = z, and integer lattice points.
for (var k = 0; k < 30; k++)
{
var t = Uniform(random, -50.0f, 50.0f);
yield return new float3(t, t, t);
yield return new float3(random.Next(-40, 41), random.Next(-40, 41), random.Next(-40, 41));
}
}
/// <summary>
/// The exact reference the enclosure tests lean on is the kernel: the double ValueDouble is snoise to
/// within float rounding at small coordinates, so a sign taken from it is the kernel's sign.
/// </summary>
[Fact]
public void ExactReferenceIsTheKernel()
{
var random = new System.Random(0xE4AC7);
for (var k = 0; k < 300; k++)
{
var p = Vector(random, -20.0f, 20.0f);
Assert.InRange(NoiseBoundTests.ValueDouble(p.x, p.y, p.z) - snoise(p), -2.0e-4, 2.0e-4);
}
}
/// <summary>
/// At SpanStage.Identity (also with an axis, stretch 1) and zero shift, offset and octaveStep, the field of
/// span_fbm and span_ridged is fbm_grad's and ridged_grad's, bit for bit, for 0 to 20 octaves (the clamp
/// included), every lacunarity and gain drawn, over random and simplex-tie points: any change to the octave
/// loop's order, scaling or accumulation, and any bit the rounding allowances put in the field, fails.
/// </summary>
[Fact]
public void SpanStagesAtIdentityAreFbmAndRidged()
{
var random = new System.Random(0x57A6E1);
foreach (var p in Points(random))
{
var octaves = random.Next(0, 21);
var lacunarity = Uniform(random, 1.4f, 2.6f);
var gain = Uniform(random, 0.3f, 0.8f);
var stage = SpanStage.Identity;
if (random.Next(2) == 1)
stage.axis = UnitVector(random) * 2.0f;
var path = new SpanPath { midpoint = p, velocity = UnitVector(random), speed = 1.0f, acceleration = 0.5f, halfLength = 0.25f };
var fbm = span_fbm(path, stage, octaves, lacunarity, gain);
var ridged = span_ridged(path, stage, octaves, lacunarity, gain);
Assert.True(SameBits(fbm.field, fbm_grad(p, octaves, lacunarity, gain)), $"span_fbm at {p}, {octaves} octaves: {fbm.field} is not fbm_grad {fbm_grad(p, octaves, lacunarity, gain)}");
Assert.True(SameBits(ridged.field, ridged_grad(p, octaves, lacunarity, gain)), $"span_ridged at {p}, {octaves} octaves: {ridged.field} is not ridged_grad {ridged_grad(p, octaves, lacunarity, gain)}");
Assert.Equal(0.0f, fbm.fold);
}
}
/// <summary>
/// The gradient is the derivative of the value through the whole transform: with every parameter off its
/// identity (frequency of either sign, stretch both ways, axis, shift, offset, octaveStep, lacunarity and
/// gain of either sign, 4 to 8 octaves) field.xyz equals the fourth-order central difference of field.w
/// along each axis, within a tolerance built from the third-derivative pin (truncation) and float32
/// rounding of the coordinates, and both agree with an independent octave sum. A dropped M or lacunarity
/// factor, M applied to the wrong side, a clamp on frequency, stretch or the octaveStep index, or a dropped
/// sign (abs of the scale or the amplitude) moves the field by the gradient itself.
/// </summary>
[Theory]
[InlineData(false)]
[InlineData(true)]
public void SpanStageGradientMatchesDifferences(bool ridged)
{
var random = new System.Random(ridged ? 0x6A1D : 0x6A1C);
var worst = 0.0;
var checkedCases = 0;
var spans = 0;
for (var k = 0; k < 6000 && spans < 80; k++)
{
var stage = Stage.DrawSmooth(random, 300.0f, 9.0f, 8, 4);
var q = Vector(random, -3.0f, 3.0f);
// The fold's crease at n = 0 is not differentiable: keep every octave's stencil clear of it.
if (ridged && Enumerable.Range(0, stage.Octaves).Any(o => abs(snoise(stage.Point(q, o))) < 0.15f))
continue;
spans++;
var span = Span.Draw(random, q, 0.1f, 0.1f);
SpanNoise At(float3 p)
{
var path = span.Path;
path.midpoint = p;
return ridged ? stage.Ridged(path) : stage.Fbm(path);
}
var centre = At(q);
// The value and gradient from an independent assembly of the octaves: pins shift, offset, octaveStep,
// M, lacunarity and gain themselves, which a value/gradient self-consistency check cannot see.
var matrix = stage.Matrix();
var referenceValue = 0.0;
var referenceGradient = new double[3];
for (var o = 0; o < stage.Octaves; o++)
{
var n = snoise_grad(stage.Point(q, o));
var sign = ridged ? -Math.Sign(n.w) : 1.0;
referenceValue += stage.Amplitude(o) * (ridged ? 1.0 - Math.Abs(n.w) : n.w);
var g = new double[] { n.x, n.y, n.z };
for (var i = 0; i < 3; i++)
for (var j = 0; j < 3; j++)
referenceGradient[i] += sign * stage.Amplitude(o) * stage.Scale(o) * matrix[i, j] * g[j];
}
var scale = 0.0;
for (var o = 0; o < stage.Octaves; o++)
scale += Math.Abs(stage.Amplitude(o) * stage.Scale(o)) * stage.Norm();
var amplitudes = Enumerable.Range(0, stage.Octaves).Sum(o => Math.Abs(stage.Amplitude(o)));
Assert.True(Math.Abs(referenceValue - centre.field.w) <= 2.0e-3 * (1.0 + amplitudes), $"stage {stage}: the value is {centre.field.w:R}, an independent octave sum gives {referenceValue:R}");
for (var i = 0; i < 3; i++)
Assert.True(Math.Abs(referenceGradient[i] - centre.field[i]) <= 1.0e-2 * (1.0 + scale), $"stage {stage}: gradient {"xyz"[i]} is {centre.field[i]:R}, an independent octave sum gives {referenceGradient[i]:R}");
// A step of 0.02 in the noise space of the fastest octave.
var fastest = stage.Norm() * Enumerable.Range(0, stage.Octaves).Max(o => Math.Abs(stage.Scale(o)));
var h = (float)(0.02 / fastest);
for (var axis = 0; axis < 3; axis++)
{
var e = new float3(axis == 0 ? h : 0.0f, axis == 1 ? h : 0.0f, axis == 2 ? h : 0.0f);
double F(float3 p) => At(p).field.w;
var step = (double)(q + e)[axis] - (q - e)[axis];
var d1 = (F(q + e) - F(q - e)) / step;
var d2 = (F(q + e * 2.0f) - F(q - e * 2.0f)) / ((double)(q + e * 2.0f)[axis] - (q - e * 2.0f)[axis]);
var estimate = (4.0 * d1 - d2) / 3.0;
var tolerance = 4.0e-2 * scale;
var error = Math.Abs(estimate - centre.field[axis]);
worst = Math.Max(worst, error / scale);
Assert.True(error <= tolerance, $"stage {stage}: d/d{"xyz"[axis]} of the value is {estimate:R} by differences and {centre.field[axis]:R} by span_{(ridged ? "ridged" : "fbm")} (tolerance {tolerance:R})");
checkedCases++;
}
}
Assert.True(checkedCases >= 120, $"only {checkedCases} derivatives were compared");
output.WriteLine($"IV-REPORT span gradient ({(ridged ? "ridged" : "fbm")}): {checkedCases} derivatives, worst error {worst:R} of the summed octave gain");
}
/// <summary>
/// stretch 1, or a zero axis, is the isotropic stage bit for bit whatever the other of the two is; a non-unit
/// axis is its normalisation (the tested axes normalise exactly); a stretch that does apply changes the
/// result.
/// </summary>
[Fact]
public void SpanStageStretchAndAxis()
{
var random = new System.Random(0x57E7C4);
for (var k = 0; k < 60; k++)
{
var stage = Stage.DrawSmooth(random, 300.0f, 15.0f, 3, 3);
var q = Vector(random, -5.0f, 5.0f);
var span = Span.Draw(random, q, 0.2f, 0.3f);
var plain = stage with { S = stage.S with { stretch = 1.0f } };
var noAxis = stage with { S = stage.S with { axis = None, stretch = 3.0f } };
var reference = plain.Ridged(span.Path);
var isotropic = (plain with { S = plain.S with { axis = None } }).Ridged(span.Path);
var viaNoAxis = noAxis.Ridged(span.Path);
var narrow = (stage with { S = stage.S with { axis = None, stretch = 0.25f } }).Ridged(span.Path);
Assert.True(SameBits(reference.field, isotropic.field) && reference.curvature == isotropic.curvature && reference.fold == isotropic.fold, "stretch 1 with an axis differs from no axis");
Assert.True(SameBits(isotropic.field, viaNoAxis.field) && isotropic.curvature == viaNoAxis.curvature && isotropic.fold == viaNoAxis.fold, "a zero axis with stretch 3 differs from the isotropic stage");
Assert.True(SameBits(isotropic.field, narrow.field) && isotropic.curvature == narrow.curvature && isotropic.fold == narrow.fold, "a zero axis with stretch 0.25 differs from the isotropic stage");
var unit = stage with { S = stage.S with { axis = new float3(0.0f, 1.0f, 0.0f) } };
var scaled = stage with { S = stage.S with { axis = new float3(0.0f, 2.0f, 0.0f) } };
var a = unit.Fbm(span.Path);
var b = scaled.Fbm(span.Path);
Assert.True(SameBits(a.field, b.field) && a.curvature == b.curvature, "a non-unit axis differs from its normalisation");
var stretched = unit.Fbm(span.Path);
Assert.False(SameBits(stretched.field, (unit with { S = unit.S with { stretch = 1.0f } }).Fbm(span.Path).field), "a stretch of " + unit.S.stretch + " did not change the field");
}
}
// The coordinate bands the enclosure tests draw: a large offset (0 keeps the drawn one) lifts S_o through
// the decades; Bucket names the decade a span's largest S_o is in.
private static readonly float[] OffsetScales = { 0.0f, 40.0f, 400.0f, 4.0e3f, 4.0e4f, 4.0e5f, 4.0e6f };
private static int Bucket(double magnitude) => magnitude < 1.0e2 ? 0 : magnitude < 1.0e3 ? 1 : magnitude < 1.0e4 ? 2 : magnitude < 1.0e5 ? 3 : magnitude < 1.0e6 ? 4 : 5;
// byStep lifts the coordinate through octaveStep (octave o adds o octaveStep) instead of offset.
private static Stage Banded(System.Random random, float offsetScale, int octaves, bool byStep)
{
var stage = Stage.Draw(random, octaves);
if (offsetScale > 0.0f)
stage = stage with { S = byStep ? stage.S with { octaveStep = Vector(random, -offsetScale, offsetScale) } : stage.S with { offset = Vector(random, -offsetScale, offsetScale) } };
return stage;
}
private static void CheckCurvature(Stage stage, bool ridged, Span span, SpanNoise result, int samples, ref double tightest, ref int checkedPoints)
{
for (var s = 0; s < samples; s++)
{
var t = span.Path.halfLength * (2.0f * s / (samples - 1) - 1.0f);
var norm = SymmetricNorm(FieldHessian(stage, span.At(t), ridged));
Assert.True(norm <= result.curvature * (1.0 + 1.0e-3), $"stage {stage}, span {span}: the Hessian norm at t = {t:R} is {norm:R}, above curvature {result.curvature:R}");
tightest = Math.Min(tightest, result.curvature > 0.0 ? 1.0 - norm / result.curvature : 1.0);
checkedPoints++;
}
}
/// <summary>
/// span-bounds-enclose, curvature: over seeded spans the operator norm of the exact field q-Hessian (chain
/// rule from the double kernel's Hessian, every octave point in double, M and the octave scales) at 64
/// points per span never exceeds curvature, in the local regime (short spans, where the third-derivative
/// term is the bound) and the global one (long spans, where SNOISE_HESSIAN is), and in every coordinate
/// band from S_o under 1e2 to over 1e6, where the float32 kernel is no longer the exact one; at
/// halfLength 0 curvature is the widened local tier, the Frobenius norm of the midpoint Hessians plus
/// e_h, capped at the pin, scaled by amplitude, f_o^2 and |M|^2; and a span through MeasureHessian's
/// witness stays enclosed in 200 directions, which the Frobenius norm at the midpoint alone does not
/// guarantee.
/// </summary>
[Theory]
[InlineData(false)]
[InlineData(true)]
public void SpanCurvatureEnclosesTheHessianAlongTheSpan(bool ridged)
{
var random = new System.Random(ridged ? 0xC0DE2 : 0xC0DE1);
var tightest = 1.0;
var checkedPoints = 0;
var buckets = new int[6];
foreach (var offsetScale in OffsetScales)
foreach (var halfLength in new[] { 0.02f, 0.5f, 8.0f })
for (var k = 0; k < 8; k++)
{
var stage = Banded(random, offsetScale, 3, k % 2 == 1);
var span = Span.Draw(random, Vector(random, -30.0f, 30.0f), halfLength, Uniform(random, 0.0f, 2.0f));
var result = ridged ? stage.Ridged(span.Path) : stage.Fbm(span.Path);
buckets[Bucket(stage.MaxMagnitude(span.Path))]++;
CheckCurvature(stage, ridged, span, result, 64, ref tightest, ref checkedPoints);
}
Assert.True(buckets.All(count => count >= 3), $"the coordinate bands drawn are {string.Join(", ", buckets)} spans: a band is not reached");
// halfLength 0: the widened local tier.
for (var k = 0; k < 40; k++)
{
var stage = Stage.DrawSmooth(random, 100.0f, 90.0f, 3);
var q = Vector(random, -30.0f, 30.0f);
var span = Span.Draw(random, q, 0.0f, 1.0f);
var result = ridged ? stage.Ridged(span.Path) : stage.Fbm(span.Path);
var expected = 0.0;
for (var o = 0; o < stage.Octaves; o++)
{
var h = snoise_hessian(stage.Point(q, o));
var frobenius = Math.Sqrt(Enumerable.Range(0, 9).Sum(i => (double)h[i / 3][i % 3] * h[i / 3][i % 3]));
var slack = Math.Pow(2.0, -10) + SNOISE_THIRD * Math.Pow(2.0, -20) * stage.Magnitude(span.Path, o);
expected += Math.Abs(stage.Amplitude(o)) * stage.Scale(o) * stage.Scale(o) * stage.Norm() * stage.Norm() * Math.Min(SNOISE_HESSIAN, frobenius + slack);
}
Assert.InRange(result.curvature, expected * (1.0 - 1.0e-3), expected * (1.0 + 1.0e-3));
}
// Through the Hessian witness: one octave at the identity transform, the span's midpoint 0.1 to 0.5 away,
// moving toward the witness at speed 1, so the witness lies on the span.
var witness = new float3(-13.7920685f, -9.672134f, 61.563934f);
var identity = new Stage(SpanStage.Identity, 1, 2.0f, 0.5f);
var through = 0.0;
for (var k = 0; k < 200; k++)
{
var u = UnitVector(random);
var distance = Uniform(random, 0.1f, 0.5f);
var span = new Span(new SpanPath { midpoint = witness + u * distance, velocity = -u, speed = 1.0f, acceleration = 0.0f, halfLength = 0.6f }, None);
var result = ridged ? identity.Ridged(span.Path) : identity.Fbm(span.Path);
CheckCurvature(identity, ridged, span, result, 64, ref tightest, ref checkedPoints);
through = Math.Max(through, SymmetricNorm(FieldHessian(identity, witness, ridged)));
}
Assert.True(through > 40.0, $"the witness Hessian norm is {through:R}: the through-the-witness spans do not reach the curvature");
output.WriteLine($"IV-REPORT span curvature ({(ridged ? "ridged" : "fbm")}): {checkedPoints} points enclosed, bands {string.Join("/", buckets)}, least relative slack {tightest:R}, witness norm {through:R}");
}
/// <summary>
/// local-tier-coefficients-unpinned: a fixed span (found by an audit, kept) at SpanStage.Identity, one
/// octave, where the Hessian's operator norm at the end of the span (11.768) exceeds the local tier built
/// with half of SNOISE_THIRD's term (about 10.1) and the exact norm along the 64 sampled points stays
/// under curvature (12.526 before the allowance).
/// </summary>
[Theory]
[InlineData(false)]
[InlineData(true)]
public void SpanCurvatureEnclosesThePinnedRankOneSpan(bool ridged)
{
var stage = new Stage(SpanStage.Identity, 1, 2.0f, 0.5f);
var span = new Span(new SpanPath { midpoint = new float3(144.93297f, -286.08800f, -53.59360f), velocity = new float3(-0.02139f, 0.17122f, -0.98500f), speed = 1.0f, acceleration = 0.0f, halfLength = 0.01f }, None);
var result = ridged ? stage.Ridged(span.Path) : stage.Fbm(span.Path);
var tightest = 1.0;
var points = 0;
CheckCurvature(stage, ridged, span, result, 64, ref tightest, ref points);
var end = SymmetricNorm(FieldHessian(stage, span.At(span.Path.halfLength), ridged));
Assert.True(end > 10.5 && result.curvature >= end, $"the end-of-span norm is {end:R} and curvature {result.curvature:R}: the span no longer pins the local tier");
}
// Octave o's amplitude-weighted value (n_o gain^o, or gain^o (1 - |n_o|) for ridged) at `point`, read from
// production: the difference of span_fbm's (span_ridged's) field at o + 1 and at o octaves.
private static double[] OctaveValues(Stage stage, float3 point, bool ridged)
{
var path = PathAt(point);
var values = new double[stage.Octaves];
var previous = 0.0;
for (var k = 1; k <= stage.Octaves; k++)
{
var total = ridged ? stage.Ridged(path, k).field.w : stage.Fbm(path, k).field.w;
values[k - 1] = total - previous;
previous = total;
}
return values;
}
// The departure of the octave sum from its midpoint branch, sum_o gain^o (s_o n_o - |n_o|) with s_o the
// sign of n_o at the midpoint, peaked over `samples` points of the span: Float for the field span_ridged
// and span_fbm themselves evaluate (read octave by octave through OctaveValues), Exact for the double
// kernel. Rounding bounds what the float differences of the accumulating sums can lose.
private static (double Float, double Exact, double Rounding) BranchDeparture(Stage stage, Span span, int samples)
{
var amplitudes = new float[stage.Octaves];
var weight = 1.0f;
var sum = 0.0;
for (var o = 0; o < stage.Octaves; o++)
{
amplitudes[o] = weight;
sum += Math.Abs(weight);
weight *= stage.Gain;
}
var middle = OctaveValues(stage, span.Path.midpoint, false);
var signs = Enumerable.Range(0, stage.Octaves).Select(o => Math.Sign(middle[o]) * Math.Sign(amplitudes[o])).ToArray();
var peakFloat = 0.0;
var peakExact = 0.0;
for (var s = 0; s < samples; s++)
{
var t = span.Path.halfLength * (2.0f * s / (samples - 1) - 1.0f);
var point = span.At(t);
var n = OctaveValues(stage, point, false);
var r = OctaveValues(stage, point, true);
var departureFloat = 0.0;
var departureExact = 0.0;
for (var o = 0; o < stage.Octaves; o++)
{
departureFloat += signs[o] * n[o] - (amplitudes[o] - r[o]);
var x = stage.PointD(point, o);
var exact = NoiseBoundTests.ValueDouble(x[0], x[1], x[2]);
departureExact += stage.Amplitude(o) * (signs[o] * exact - Math.Abs(exact));
}
peakFloat = Math.Max(peakFloat, Math.Abs(departureFloat));
peakExact = Math.Max(peakExact, Math.Abs(departureExact));
}
return (peakFloat, peakExact, stage.Octaves * Math.Pow(2.0, -22) * sum);
}
// Newton steps from q to the exact kernel's n = 0 (the float gradient steers, the double value decides).
private static float3 ProjectToZero(float3 q)
{
for (var iteration = 0; iteration < 4; iteration++)
{
var n = snoise_grad(q);
var gg = (double)n.x * n.x + (double)n.y * n.y + (double)n.z * n.z;
if (gg < 1.0e-12)
break;
var v = NoiseBoundTests.ValueDouble(q.x, q.y, q.z);
q = new float3((float)(q.x - v * n.x / gg), (float)(q.y - v * n.y / gg), (float)(q.z - v * n.z / gg));
}
return q;
}
/// <summary>
/// span-bounds-enclose, fold: for span_ridged, the value's departure from the midpoint branch (each
/// octave's gain^o (s n - |n|), s the sign of n at the midpoint, summed) never exceeds fold at 64 points of
/// every span, measured twice: for the float field span_ridged itself evaluates at the path point, and for
/// the exact octave sum (double kernel); in every coordinate band from S_o under 1e2 to over 1e6. fold is
/// positive on spans that do cross zero (and the departure there is not trivially zero), exactly 0 for
/// span_fbm and when no octave can reach zero.
/// </summary>
[Fact]
public void SpanFoldEnclosesTheBranchDeviation()
{
var random = new System.Random(0xF01D);
var crossing = 0;
var largest = 0.0;
var tightestFloat = 1.0;
var tightestExact = 1.0;
var buckets = new int[6];
foreach (var offsetScale in OffsetScales)
for (var k = 0; k < 60; k++)
{
var stage = Banded(random, offsetScale, 3, k / 2 % 2 == 1);
var halfLength = new[] { 0.02f, 0.1f, 0.5f, 2.0f }[random.Next(4)];
// Draw midpoints with small noise values so spans cross zero, as well as arbitrary ones.
var q = Vector(random, -30.0f, 30.0f);
if (k % 2 == 0)
{
var best = q;
var smallest = float.MaxValue;
for (var trial = 0; trial < 40; trial++)
{
var candidate = Vector(random, -30.0f, 30.0f);
var value = abs(snoise(stage.Point(candidate, 0)));
if (value < smallest)
(smallest, best) = (value, candidate);
}
q = best;
}
var span = Span.Draw(random, q, halfLength, Uniform(random, 0.0f, 2.0f));
buckets[Bucket(stage.MaxMagnitude(span.Path))]++;
var fold = stage.Ridged(span.Path).fold;
Assert.Equal(0.0f, stage.Fbm(span.Path).fold);
var (peakFloat, peakExact, rounding) = BranchDeparture(stage, span, 64);
Assert.True(peakFloat <= fold * (1.0 + 1.0e-3) + 1.0e-6 + rounding, $"stage {stage}, span {span}: the float branch departure is {peakFloat:R}, above fold {fold:R}");
Assert.True(peakExact <= fold * (1.0 + 1.0e-3) + 1.0e-6, $"stage {stage}, span {span}: the exact branch departure is {peakExact:R}, above fold {fold:R}");
if (fold > 0.0)
crossing++;
largest = Math.Max(largest, Math.Max(peakFloat, peakExact));
if (fold > 0.0 && peakFloat > 0.0)
tightestFloat = Math.Min(tightestFloat, peakFloat / fold);
if (fold > 0.0 && peakExact > 0.0)
tightestExact = Math.Min(tightestExact, peakExact / fold);
}
Assert.True(buckets.All(count => count >= 3), $"the coordinate bands drawn are {string.Join(", ", buckets)} spans: a band is not reached");
Assert.True(crossing >= 100, $"only {crossing} of {OffsetScales.Length * 60} spans had a positive fold");
Assert.True(largest > 0.1, $"the largest branch departure drawn is {largest:R}: the check has no teeth");
output.WriteLine($"IV-REPORT span fold: {crossing} spans with fold > 0, bands {string.Join("/", buckets)}, largest departure {largest:R}, tightest departure / fold {tightestFloat:R} (float) {tightestExact:R} (exact)");
// An octave that cannot reach zero: a midpoint where |n| is large and a short span give fold exactly 0.
var calm = new Stage(SpanStage.Identity, 1, 2.0f, 0.5f);
var found = 0;
for (var k = 0; k < 400 && found < 20; k++)
{
var q = Vector(random, -30.0f, 30.0f);
if (abs(snoise(q)) < 0.6f)
continue;
found++;
Assert.Equal(0.0f, calm.Ridged(new SpanPath { midpoint = q, velocity = new float3(1.0f, 0.0f, 0.0f), speed = 1.0f, acceleration = 0.0f, halfLength = 0.01f }).fold);
}
Assert.True(found >= 5, $"only {found} calm midpoints found");
}
/// <summary>
/// The Gl reach term over a long span: a straight path from rest along -s grad n (s the sign of n at the
/// midpoint), identity stage, one octave, half-length 6, found by an audit and kept. The path reaches
/// n = 0 and the branch departure (about 0.387) is covered only if the gradient bound grows by
/// Hl nu halfLength over the span; a reach that shrinks with halfLength (nu / halfLength) puts fold at
/// 0.384 and fails here, as does any reach right only for half-lengths up to sqrt 2.
/// </summary>
[Fact]
public void SpanFoldEnclosesTheLongPathFromRest()
{
var stage = new Stage(SpanStage.Identity, 1, 2.0f, 0.5f);
const float halfLength = 6.0f;
const float acceleration = 0.0055555557f;
var span = new Span(new SpanPath { midpoint = new float3(7.634474f, 24.228996f, -3.7595012f), velocity = None, speed = acceleration * halfLength, acceleration = acceleration, halfLength = halfLength }, new float3(-0.9375727f, -0.07116994f, 0.3404294f) * acceleration);
var fold = stage.Ridged(span.Path).fold;
var (departureFloat, departureExact, rounding) = BranchDeparture(stage, span, 129);
Assert.True(Math.Max(departureFloat, departureExact) > 0.3, $"the branch departure is {Math.Max(departureFloat, departureExact):R}: the span no longer crosses zero");
Assert.True(departureFloat <= fold * (1.0 + 1.0e-3) + 1.0e-6 + rounding, $"the float branch departure is {departureFloat:R}, above fold {fold:R}");
Assert.True(departureExact <= fold * (1.0 + 1.0e-3) + 1.0e-6, $"the exact branch departure is {departureExact:R}, above fold {fold:R}");
}
/// <summary>
/// A midpoint within float noise of n = 0, on a path that does not move: halfLength 0 (a point), and
/// halfLength 1 under a speed bound far above the path's motion (speed is an upper bound, so that is a
/// legal SpanPath, and composed chains pass loose ones). The exact kernel's sign at such a midpoint is
/// opposite to the float one in about a fifth of the draws, so fold must still cover the exact
/// departure there: a fold forced to 0 at halfLength 0, or an S_o whose reach term drives e_n below
/// zero under the loose bound, leaves it 0 and fails. span_fbm's fold is 0 throughout.
/// Near the origin S_o is tiny (scales 1e-4 and 1e-2 of the coordinate), and only the value allowance's
/// absolute term 2^-14 carries the kernel's own rounding: the exact departure there is the coordinate's
/// rounding times the gradient (1e-11 to 1e-9), so the exact check's additive tolerance is 1e-13, under it.
/// A larger tolerance absorbs a dropped absolute term whole.
/// </summary>
[Theory]
[InlineData(false)]
[InlineData(true)]
public void SpanFoldEnclosesTheExactBranchAtANearZeroMidpoint(bool loose)
{
var random = new System.Random(loose ? 0x100E : 0x4A11);
var stage = new Stage(SpanStage.Identity, 1, 2.0f, 0.5f);
var crossing = 0;
foreach (var scale in new[] { 1.0e-4f, 1.0e-2f, 30.0f, 300.0f, 3000.0f })
for (var k = 0; k < 300; k++)
{
var q = ProjectToZero(Vector(random, -scale, scale));
var speed = loose ? 4.0f * (abs(q.x) + abs(q.y) + abs(q.z) + 20.0f) : 0.0f;
var span = new Span(new SpanPath { midpoint = q, velocity = None, speed = speed, acceleration = 0.0f, halfLength = loose ? 1.0f : 0.0f }, None);
var fold = stage.Ridged(span.Path).fold;
var (departureFloat, departureExact, rounding) = BranchDeparture(stage, span, 2);
Assert.Equal(0.0f, stage.Fbm(span.Path).fold);
Assert.True(departureExact <= fold * (1.0 + 1.0e-3) + 1.0e-13, $"midpoint {q}, speed {speed:R}: the exact branch departure is {departureExact:R}, above fold {fold:R}");
Assert.True(departureFloat <= fold * (1.0 + 1.0e-3) + 1.0e-9 + rounding, $"midpoint {q}, speed {speed:R}: the float branch departure is {departureFloat:R}, above fold {fold:R}");
if (departureExact > 0.0)
crossing++;
}
Assert.True(crossing >= 80, $"only {crossing} of 1500 midpoints had an exact sign opposite to the float one: the check has no teeth");
}
/// <summary>
/// The step term of S_o, alone: an octaveStep of random direction and magnitude 1e3 to 1e5 (log-uniform) at
/// the identity stage, two octaves, the midpoint moved until octave 1's coordinate (2 q + step) sits on
/// the kernel's n = 0 (a Newton projection of that coordinate, then q = (x - step) / 2). There the
/// coordinate's float rounding (ulp 0.004 to 0.008) is what separates the float kernel from the exact one,
/// and only |o octaveStep|_1 in S_o carries it: the exact branch departure stays within fold, span_fbm's
/// fold is 0, and at least 100 draws have an exact sign opposite to the float one. Without the step term
/// fold is the 2^-14 and the small-coordinate rounding alone, and the departure passes it.
/// </summary>
[Fact]
public void SpanFoldEnclosesTheExactBranchAtLargeOctaveSteps()
{
var random = new System.Random(0x57E9);
var crossing = 0;
for (var k = 0; k < 600; k++)
{
var step = UnitVector(random) * LogUniform(random, 1.0e3f, 1.0e5f);
var transform = SpanStage.Identity;
transform.octaveStep = step;
var stage = new Stage(transform, 2, 2.0f, 0.5f);
var q = (ProjectToZero(Vector(random, -30.0f, 30.0f) + step) - step) * 0.5f;
var span = new Span(new SpanPath { midpoint = q, velocity = None, speed = 0.0f, acceleration = 0.0f, halfLength = 0.0f }, None);
var fold = stage.Ridged(span.Path).fold;
var (_, departureExact, _) = BranchDeparture(stage, span, 2);
Assert.Equal(0.0f, stage.Fbm(span.Path).fold);
Assert.True(departureExact <= fold * (1.0 + 1.0e-3) + 1.0e-9, $"octaveStep {step}, midpoint {q}: the exact branch departure is {departureExact:R}, above fold {fold:R}");
if (departureExact > 0.0)
crossing++;
}
Assert.True(crossing >= 100, $"only {crossing} of 600 midpoints had an exact sign opposite to the float one: the check has no teeth");
}
/// <summary>
/// The fold's acceleration term, alone, at the identity stage and under a transform (frequency 2.5, stretch
/// 0.5, shift and offset): the midpoint velocity is zero and the acceleration points along the gradient
/// (in q), against the sign of n, so n reaches zero and crosses it on the way to the span's ends by
/// acceleration only. Without the acceleration term, with the gradient bound under the gradient the path
/// reaches, or with the transform's norm f_o |M| misplaced in it, the departure from the midpoint branch
/// exceeds fold.
/// </summary>
[Theory]
[InlineData(false)]
[InlineData(true)]
public void SpanFoldCoversAccelerationAlone(bool transformed)
{
var random = new System.Random(transformed ? 0xACCE2 : 0xACCE1);
var stage = transformed
? new Stage(new SpanStage { frequency = 2.5f, stretch = 0.5f, axis = new float3(0.0f, 1.0f, 0.0f), shift = new float3(0.3f, -1.2f, 2.1f), offset = new float3(5.0f, -3.0f, 7.0f), octaveStep = None }, 1, 2.0f, 0.5f)
: new Stage(SpanStage.Identity, 1, 2.0f, 0.5f);
var m = stage.Matrix();
var cases = 0;
var crossings = 0;
var tightest = double.MaxValue;
for (var trial = 0; trial < 400_000 && cases < 150; trial++)
{
var q = Vector(random, -60.0f, 60.0f);
var n = snoise_grad(stage.Point(q, 0));
var gradient = new float3((float)(m[0, 0] * n.x + m[0, 1] * n.y + m[0, 2] * n.z), (float)(m[1, 0] * n.x + m[1, 1] * n.y + m[1, 2] * n.z), (float)(m[2, 0] * n.x + m[2, 1] * n.y + m[2, 2] * n.z));
var halfLength = Uniform(random, 0.05f, 0.15f);
var acceleration = Uniform(random, 4.0f, 12.0f);
// n at the midpoint close enough to zero that the acceleration alone carries it across.
if (length(gradient) < 1.5f || abs(n.w) > 0.4f * length(gradient) * acceleration * halfLength * halfLength * 0.5f)
continue;
var c = -(float)Math.Sign(n.w) * acceleration * normalize(gradient);
var span = new Span(new SpanPath { midpoint = q, velocity = None, speed = acceleration * halfLength, acceleration = acceleration, halfLength = halfLength }, c);
var fold = stage.Ridged(span.Path).fold;
var departure = 0.0;
for (var s = 0; s < 64; s++)
{
var value = (double)snoise(stage.Point(span.At(halfLength * (2.0f * s / 63 - 1.0f)), 0));
departure = Math.Max(departure, Math.Abs(Math.Sign(n.w) * value - Math.Abs(value)));
}
if (departure > 0.0)
crossings++;
Assert.True(departure <= fold * (1.0 + 1.0e-3), $"q {q}, acceleration {acceleration:R}, half-length {halfLength:R}: the branch departure is {departure:R}, above fold {fold:R}");
if (departure > 0.0)
tightest = Math.Min(tightest, departure / fold);
cases++;
}
Assert.True(crossings >= 100, $"only {crossings} of {cases} paths crossed zero by acceleration alone");
output.WriteLine($"IV-REPORT span fold, acceleration only ({(transformed ? "transformed" : "identity")}): {crossings} of {cases} paths cross zero, tightest departure / fold {tightest:R}");
}
/// <summary>
/// The fold's velocity term under a stretch: the velocity runs along the stretched axis (M scales it by
/// 1/stretch), the midpoint is where n reaches zero inside the span, so the departure is about twice
/// |grad n . w| halfLength and a velocity stretched the wrong way (or not at all) puts fold under it.
/// </summary>
[Theory]
[InlineData(0.25f)]
[InlineData(0.5f)]
[InlineData(4.0f)]
public void SpanFoldCoversTheStretchedVelocity(float stretch)
{
var random = new System.Random(0x57E1C + (int)(stretch * 8.0f));
var cases = 0;
var crossings = 0;
var tightest = double.MaxValue;
for (var trial = 0; trial < 400_000 && cases < 60; trial++)
{
var axis = UnitVector(random);
var stage = new Stage(new SpanStage { frequency = Uniform(random, 0.5f, 2.0f), stretch = stretch, axis = axis, shift = Vector(random, -4.0f, 4.0f), offset = Vector(random, -4.0f, 4.0f), octaveStep = None }, 1, 2.0f, 0.5f);
var q = Vector(random, -30.0f, 30.0f);
var velocity = axis * Uniform(random, 0.5f, 2.0f);
var halfLength = 0.05f;
var path = new SpanPath { midpoint = q, velocity = velocity, speed = length(velocity), acceleration = 0.0f, halfLength = halfLength };
var n = snoise_grad(stage.Point(q, 0));
var m = stage.Matrix();
var w = new double[3];
for (var i = 0; i < 3; i++)
w[i] = m[i, 0] * velocity.x + m[i, 1] * velocity.y + m[i, 2] * velocity.z;
var along = Math.Abs(n.x * w[0] + n.y * w[1] + n.z * w[2]);
if (along < 2.0 || Math.Abs(n.w) > 0.5 * along * halfLength)
continue;
var fold = stage.Ridged(path).fold;
var departure = 0.0;
for (var s = 0; s < 64; s++)
{
var t = halfLength * (2.0f * s / 63 - 1.0f);
var value = (double)snoise(stage.Point(q + velocity * t, 0));
departure = Math.Max(departure, Math.Abs(Math.Sign(n.w) * value - Math.Abs(value)));
}
if (departure > 0.0)
crossings++;
Assert.True(departure <= fold * (1.0 + 1.0e-3), $"stage {stage}, q {q}: the branch departure is {departure:R}, above fold {fold:R}");
if (departure > 0.0)
tightest = Math.Min(tightest, departure / fold);
cases++;
}
Assert.True(crossings >= 30, $"only {crossings} of {cases} paths crossed zero along the stretched velocity");
output.WriteLine($"IV-REPORT span fold, stretched velocity {stretch}: {crossings} of {cases} paths cross zero, tightest departure / fold {tightest:R}");
}
/// <summary>
/// Both bounds are sums over the octaves, weighted by |gain|^o: for a stage of several octaves,
/// curvature and fold equal the weighted sum of the bounds of each octave taken as a stage of its own,
/// to within the rounding of forming the octave's coordinate another way. A weight misplaced (gain
/// instead of |gain|, a division for a product), a dropped octave or a mis-scaled octave fails.
/// </summary>
[Fact]
public void SpanBoundsAreSumsOverTheOctaves()
{
var random = new System.Random(0x0C7A5);
var compared = 0;
for (var k = 0; k < 80; k++)
{
var stage = Stage.DrawSmooth(random, 300.0f, 90.0f, 5, 2);
var span = Span.Draw(random, Vector(random, -30.0f, 30.0f), new[] { 0.02f, 0.1f, 0.5f }[k % 3], Uniform(random, 0.0f, 2.0f));
var whole = stage.Ridged(span.Path);
double curvature = 0.0, fold = 0.0;
for (var o = 0; o < stage.Octaves; o++)
{
var one = stage.OneOctave(o).Ridged(span.Path);
var weight = Math.Abs(stage.Amplitude(o));
curvature += weight * one.curvature;
fold += weight * one.fold;
}
Assert.InRange(whole.curvature, curvature * (1.0 - 5.0e-3) - 1.0e-4, curvature * (1.0 + 5.0e-3) + 1.0e-4);
Assert.InRange(whole.fold, fold * (1.0 - 5.0e-3) - 1.0e-4, fold * (1.0 + 5.0e-3) + 1.0e-4);
if (fold > 1.0e-3)
compared++;
}
Assert.True(compared >= 20, $"only {compared} of 80 stages had a fold to compare");
}
}

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

File diff suppressed because it is too large Load diff

17
packages/cultmath/tools/compile-glsl.ps1
View file

@ -29,12 +29,21 @@ if ($Phacelle) {
if ($Cycles) {
$library += Get-Content -LiteralPath (Join-Path $repoRoot "shaders\CultMath.Cycles.glsl") -Raw
}
# A parameter of a struct type the library declares is passed as that struct's constructor with each
# field's zero, read from its declaration (a field may be a struct itself); every other type is its zero.
$structs = @{}
foreach ($declaration in [regex]::Matches($library, '(?s)struct (\w+)\s*\{(.*?)\};')) {
$structs[$declaration.Groups[1].Value] = @([regex]::Matches($declaration.Groups[2].Value, '(?m)^\s*(\w+)\s+\w+\s*;') | ForEach-Object { $_.Groups[1].Value })
}
function Get-Zero([string]$type) {
if ($structs.ContainsKey($type)) { return "$type($((@($structs[$type] | ForEach-Object { Get-Zero $_ })) -join ', '))" }
if ($type -eq 'int' -or $type -like 'ivec*') { return "$type(0)" }
if ($type -eq 'uint' -or $type -like 'uvec*') { return "$type(0u)" }
return "$type(0.0)"
}
$calls = foreach ($match in [regex]::Matches($library, '(?m)^\w+\s+(cultmath_\w+)\s*\(([^)]*)\)')) {
$arguments = foreach ($parameter in ($match.Groups[2].Value -split ',')) {
$type = ($parameter.Trim() -split '\s+')[0]
if ($type -eq 'int' -or $type -like 'ivec*') { "$type(0)" }
elseif ($type -eq 'uint' -or $type -like 'uvec*') { "$type(0u)" }
else { "$type(0.0)" }
Get-Zero (($parameter.Trim() -split '\s+')[0])
}
" $($match.Groups[1].Value)($($arguments -join ', '));"
}

42
packages/cultmath/tools/compile-hlsl.ps1 Normal file
View file

@ -0,0 +1,42 @@
param(
# The HLSL files to compile, each as a lib_6_3 library with its own directory as the include path.
# Default: both copies of CultMath.hlsl (the canonical one and the Unity package's byte copy).
# tools\compile-hlsl.ps1 -ShaderPath <a file holding 'float f( { return 0; }'> is the negative control.
[string[]]$ShaderPath = @(),
# The dxc to run. Default: the x64 install get-dxc.ps1 makes; it is never searched for, so a dxc
# that cannot launch (one built for another architecture) is an error, not a skipped check.
[string]$DxcPath = ""
)
$ErrorActionPreference = "Stop"
$repoRoot = Resolve-Path (Join-Path $PSScriptRoot "..")
if (-not $DxcPath) {
$DxcPath = Join-Path $repoRoot ".tools\dxc\current\bin\x64\dxc.exe"
}
if (-not $ShaderPath) {
$ShaderPath = @(
(Join-Path $repoRoot "shaders\CultMath.hlsl"),
(Join-Path $repoRoot "unity\org.gamecult.cultmath\Shaders\CultMath.hlsl"))
}
if (-not (Test-Path -LiteralPath $DxcPath)) {
throw "dxc not found (parameter DxcPath). Run tools\get-dxc.ps1 first or pass -DxcPath."
}
$outputDirectory = Join-Path $repoRoot ".tools\hlsl"
New-Item -ItemType Directory -Force -Path $outputDirectory | Out-Null
$index = 0
foreach ($shader in $ShaderPath) {
$index++
if (-not (Test-Path -LiteralPath $shader -PathType Leaf)) {
throw "shader not found (parameter ShaderPath, entry $index)."
}
$shaderFile = Get-Item -LiteralPath $shader
$output = Join-Path $outputDirectory "$index-$($shaderFile.BaseName).dxil"
& $DxcPath -T lib_6_3 -HV 2021 -WX -I $shaderFile.DirectoryName $shaderFile.FullName -Fo $output
if ($LASTEXITCODE -ne 0) {
throw "dxc failed with exit code $LASTEXITCODE (parameter ShaderPath, entry $index)."
}
Write-Host "compiled shader $index of $(@($ShaderPath).Count) (lib_6_3, -HV 2021, -WX)"
}

7
packages/cultmath/tools/compile-hlsl.ps1.meta Normal file
View file

@ -0,0 +1,7 @@
fileFormatVersion: 2
guid: 3f6c1e7a9b2d4c58a0e4d7b1c95f2a86
DefaultImporter:
externalObjects: {}
userData:
assetBundleName:
assetBundleVariant:

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

@ -14,7 +14,7 @@ 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;
static const float CULTMATH_SNOISE_THIRD = 535.8009;
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); }

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

@ -302,10 +302,10 @@ float3x3 cultmath_snoise_hessian(float3 value)
// 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);
if (b == 0.0)
return t1 != t0 ? clamp(a, lo, hi) * (t1 - t0) : 0.0;
float direction = t1 < t0 ? -1.0 : 1.0;
float start = min(t0, t1);
float stop = max(t0, t1);
@ -315,7 +315,9 @@ float cultmath_integral_clamp_linear(float a, float b, float lo, float hi, float
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);
float held = enter > start ? before * (enter - start) : 0.0;
float kept = stop > leave ? after * (stop - leave) : 0.0;
float total = held + (a + b * (0.5 * (enter + leave))) * (leave - enter) + kept;
return direction * total;
}
@ -603,6 +605,119 @@ CultCellular cultmath_cellular(float3 p)
#include "CultMath.Cycles.hlsl"
#include "CultMath.Affine.hlsl"
// Span stages (invariant 8): the SpanPath, SpanStage and SpanNoise of math.Span.cs, field for field in
// C#'s order, and cultmath_span_stage, its octave stage, statement for statement so the text is a bit-for-bit
// mirror. They follow the includes because the bounds read CULTMATH_SNOISE_LIPSCHITZ, _HESSIAN and _THIRD.
// A shader caller builds a SpanStage itself (SpanStage(1.0, 1.0, vec3(0.0), vec3(0.0), vec3(0.0), vec3(0.0))
// is the identity of math.SpanStage.Identity, which has no shader form); see math.span_stage's comments for
// the span, the transform, the bounds and the allowances.
struct SpanPath
{
float3 midpoint;
float3 velocity;
float speed;
float acceleration;
float halfLength;
};
struct SpanStage
{
float frequency;
float stretch;
float3 axis;
float3 shift;
float3 offset;
float3 octaveStep;
};
struct SpanNoise
{
float4 field;
float curvature;
float fold;
};
// The octave stage behind cultmath_span_fbm and cultmath_span_ridged, as math.span_stage is behind span_fbm
// and span_ridged. The Frobenius norm of the Hessian is written out in C#'s order (GLSL has no _mRC member
// access); the C# local named step is octaveShift here (it would shadow GLSL's step()).
SpanNoise cultmath_span_stage(bool ridged, SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain)
{
float3 q = path.midpoint;
float halfLength = path.halfLength;
octaves = clamp(octaves, 0, 16);
bool stretched = dot(stage.axis, stage.axis) > 0.0 && stage.stretch != 1.0;
float3 direction = stretched ? normalize(stage.axis) : float3(0.0, 0.0, 0.0);
float squeeze = stretched ? 1.0 / stage.stretch - 1.0 : 0.0;
float norm = abs(stage.frequency) * (stretched ? max(1.0, 1.0 / stage.stretch) : 1.0);
float3 mq = stretched ? stage.frequency * (q + squeeze * dot(direction, q) * direction) : stage.frequency * q;
float3 mv = stretched ? stage.frequency * (path.velocity + squeeze * dot(direction, path.velocity) * direction) : stage.frequency * path.velocity;
float amplitude = 1.0;
float scale = 1.0;
float value = 0.0;
float3 gradient = float3(0.0, 0.0, 0.0);
float curvature = 0.0;
float fold = 0.0;
for (int i = 0; i < octaves; i++)
{
float3 octaveShift = (float)i * stage.octaveStep;
float3 x = scale * (mq + stage.shift) + (stage.offset + octaveShift);
float4 n = cultmath_snoise_grad(x);
float3 g = float3(n.x, n.y, n.z);
float3 mg = stretched ? stage.frequency * (g + squeeze * dot(direction, g) * direction) : stage.frequency * g;
if (ridged)
{
float s = sign(n.w);
value += amplitude * (1.0 - abs(n.w));
gradient += -amplitude * scale * s * mg;
}
else
{
value += amplitude * n.w;
gradient += amplitude * scale * mg;
}
float weight = abs(amplitude);
float factor = abs(scale) * norm;
float nu = factor * path.speed;
float rounding = 9.5367431640625e-7 * (abs(scale) * (cultmath_csum(abs(mq)) + cultmath_csum(abs(stage.shift))) + cultmath_csum(abs(stage.offset)) + cultmath_csum(abs(octaveShift)) + 1.7320508 * nu * halfLength);
float valueSlack = 6.103515625e-5 + CULTMATH_SNOISE_LIPSCHITZ * rounding;
float gradientSlack = 2.44140625e-4 + CULTMATH_SNOISE_HESSIAN * rounding;
float hessianSlack = 9.765625e-4 + CULTMATH_SNOISE_THIRD * rounding;
float3x3 hm = cultmath_snoise_hessian(x);
float frobenius = sqrt(hm[0][0] * hm[0][0] + hm[0][1] * hm[0][1] + hm[0][2] * hm[0][2] + hm[1][0] * hm[1][0] + hm[1][1] * hm[1][1] + hm[1][2] * hm[1][2]
+ hm[2][0] * hm[2][0] + hm[2][1] * hm[2][1] + hm[2][2] * hm[2][2]);
float hessianBound = min(CULTMATH_SNOISE_HESSIAN, frobenius + hessianSlack + CULTMATH_SNOISE_THIRD * nu * halfLength);
curvature += weight * factor * factor * hessianBound;
if (ridged)
{
float3 w = scale * mv;
float reach = hessianBound * nu * halfLength;
float gl = min(CULTMATH_SNOISE_LIPSCHITZ, length(g) + gradientSlack + reach);
float travel = (abs(dot(g, w)) + gradientSlack * length(w)) * halfLength + (hessianBound * nu * length(w) + gl * factor * path.acceleration) * (0.5 * halfLength * halfLength) + 2.0 * valueSlack;
fold += 2.0 * weight * max(travel - abs(n.w), 0.0);
}
scale *= lacunarity;
amplitude *= gain;
}
SpanNoise result;
result.field = float4(gradient, value);
result.curvature = curvature;
result.fold = fold;
return result;
}
SpanNoise cultmath_span_fbm(SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain)
{
return cultmath_span_stage(false, path, stage, octaves, lacunarity, gain);
}
SpanNoise cultmath_span_ridged(SpanPath path, SpanStage stage, int octaves, float lacunarity, float gain)
{
return cultmath_span_stage(true, path, stage, octaves, lacunarity, gain);
}
float cultmath_value_noise(float2 position)
{
float2 cell = floor(position);

29
src/GameCult.Caching.MessagePack/CultMathJson.cs
View file

@ -14,7 +14,7 @@ namespace GameCult.Caching.MessagePack;
/// The canonical System.Text.Json shape of every public CultMath value type, the JSON twin of
/// <see cref="CultMathResolver"/>. Vectors, quaternions and colors are objects keyed by their
/// field names (<c>{"x":1,"y":2}</c>), matrices are arrays of row objects, <c>rect</c> is
/// <c>{"min":..,"max":..}</c>, <c>CultCellular</c> is <c>{"nearest":..,"edge":..,"id":..}</c>, and <c>CultPhasor</c> is <c>{"cos":..,"sin":..}</c>.
/// <c>{"min":..,"max":..}</c>, <c>CultCellular</c> is <c>{"nearest":..,"edge":..,"id":..}</c>, <c>CultPhasor</c> is <c>{"cos":..,"sin":..}</c>, <c>SpanNoise</c> is <c>{"field":..,"curvature":..,"fold":..}</c>, <c>SpanPath</c> is <c>{"midpoint":..,"velocity":..,"speed":..,"acceleration":..,"halfLength":..}</c>, and <c>SpanStage</c> is <c>{"frequency":..,"stretch":..,"axis":..,"shift":..,"offset":..,"octaveStep":..}</c>.
/// Property names are fixed: <c>PropertyNamingPolicy</c> does not rename them, and
/// <c>PropertyNameCaseInsensitive</c> makes reading case-insensitive.
/// <para>
@ -481,7 +481,7 @@ public static class CultMathJson
Add<Random, uint>(UIntScalar.Instance, new Fields<uint>(("state", Val<uint>())),
a => new Random { state = a[0] }, (v, a) => a[0] = v.state),
// Compositions: matrices are arrays of row objects, rect is {min, max}, CultCellular is {nearest, edge, id}, CultPhasor is {cos, sin}.
// Compositions: matrices are arrays of row objects, rect is {min, max}, CultCellular is {nearest, edge, id}, CultPhasor is {cos, sin}, SpanNoise is {field, curvature, fold}, SpanPath is {midpoint, velocity, speed, acceleration, halfLength}, SpanStage is {frequency, stretch, axis, shift, offset, octaveStep}.
Add<float2x2, float>(f, new Rows<float>(2, Xy<float>()),
a => new float2x2(a[0], a[1], a[2], a[3]),
(v, a) => { a[0] = v[0].x; a[1] = v[0].y; a[2] = v[1].x; a[3] = v[1].y; }),
@ -514,6 +514,31 @@ public static class CultMathJson
a[0] = v.cos.x; a[1] = v.cos.y; a[2] = v.cos.z; a[3] = v.cos.w;
a[4] = v.sin.x; a[5] = v.sin.y; a[6] = v.sin.z; a[7] = v.sin.w;
}),
Add<SpanNoise, float>(f, new Fields<float>(("field", Xyzw<float>()), ("curvature", Val<float>()), ("fold", Val<float>())),
a => new SpanNoise(new float4(a[0], a[1], a[2], a[3]), a[4], a[5]),
(v, a) =>
{
a[0] = v.field.x; a[1] = v.field.y; a[2] = v.field.z; a[3] = v.field.w;
a[4] = v.curvature; a[5] = v.fold;
}),
Add<SpanPath, float>(f, new Fields<float>(("midpoint", Xyz<float>()), ("velocity", Xyz<float>()), ("speed", Val<float>()), ("acceleration", Val<float>()), ("halfLength", Val<float>())),
a => new SpanPath { midpoint = new float3(a[0], a[1], a[2]), velocity = new float3(a[3], a[4], a[5]), speed = a[6], acceleration = a[7], halfLength = a[8] },
(v, a) =>
{
a[0] = v.midpoint.x; a[1] = v.midpoint.y; a[2] = v.midpoint.z;
a[3] = v.velocity.x; a[4] = v.velocity.y; a[5] = v.velocity.z;
a[6] = v.speed; a[7] = v.acceleration; a[8] = v.halfLength;
}),
Add<SpanStage, float>(f, new Fields<float>(("frequency", Val<float>()), ("stretch", Val<float>()), ("axis", Xyz<float>()), ("shift", Xyz<float>()), ("offset", Xyz<float>()), ("octaveStep", Xyz<float>())),
a => new SpanStage { frequency = a[0], stretch = a[1], axis = new float3(a[2], a[3], a[4]), shift = new float3(a[5], a[6], a[7]), offset = new float3(a[8], a[9], a[10]), octaveStep = new float3(a[11], a[12], a[13]) },
(v, a) =>
{
a[0] = v.frequency; a[1] = v.stretch;
a[2] = v.axis.x; a[3] = v.axis.y; a[4] = v.axis.z;
a[5] = v.shift.x; a[6] = v.shift.y; a[7] = v.shift.z;
a[8] = v.offset.x; a[9] = v.offset.y; a[10] = v.offset.z;
a[11] = v.octaveStep.x; a[12] = v.octaveStep.y; a[13] = v.octaveStep.z;
}),
};
var map = new Dictionary<Type, JsonConverter>();

55
src/GameCult.Caching.MessagePack/CultMathResolver.cs
View file

@ -180,7 +180,7 @@ public sealed class CultMathResolver : IFormatterResolver
return value;
}),
// Compositions: matrices are rows, rect is [min, max], CultCellular is [nearest, edge, id], CultPhasor is [cos, sin].
// Compositions: matrices are rows, rect is [min, max], CultCellular is [nearest, edge, id], CultPhasor is [cos, sin], SpanNoise is [field, curvature, fold], SpanPath is [midpoint, velocity, speed, acceleration, halfLength], SpanStage is [frequency, stretch, axis, shift, offset, octaveStep].
Shape<float2x2>((ref MessagePackWriter w, float2x2 v) => { w.WriteArrayHeader(2); PutFloat2(ref w, v[0]); PutFloat2(ref w, v[1]); },
(ref MessagePackReader r) =>
{
@ -236,6 +236,59 @@ public sealed class CultMathResolver : IFormatterResolver
}
return value;
}),
Shape<SpanNoise>((ref MessagePackWriter w, SpanNoise v) => { w.WriteArrayHeader(3); PutFloat4(ref w, v.field); w.Write(v.curvature); w.Write(v.fold); },
(ref MessagePackReader r) =>
{
var value = default(SpanNoise);
for (int i = 0, n = r.ReadArrayHeader(); i < n; i++)
{
switch (i)
{
case 0: value.field = GetFloat4(ref r); break;
case 1: value.curvature = r.ReadSingle(); break;
case 2: value.fold = r.ReadSingle(); break;
default: r.Skip(); break;
}
}
return value;
}),
Shape<SpanPath>((ref MessagePackWriter w, SpanPath v) => { w.WriteArrayHeader(5); PutFloat3(ref w, v.midpoint); PutFloat3(ref w, v.velocity); w.Write(v.speed); w.Write(v.acceleration); w.Write(v.halfLength); },
(ref MessagePackReader r) =>
{
var value = default(SpanPath);
for (int i = 0, n = r.ReadArrayHeader(); i < n; i++)
{
switch (i)
{
case 0: value.midpoint = GetFloat3(ref r); break;
case 1: value.velocity = GetFloat3(ref r); break;
case 2: value.speed = r.ReadSingle(); break;
case 3: value.acceleration = r.ReadSingle(); break;
case 4: value.halfLength = r.ReadSingle(); break;
default: r.Skip(); break;
}
}
return value;
}),
Shape<SpanStage>((ref MessagePackWriter w, SpanStage v) => { w.WriteArrayHeader(6); w.Write(v.frequency); w.Write(v.stretch); PutFloat3(ref w, v.axis); PutFloat3(ref w, v.shift); PutFloat3(ref w, v.offset); PutFloat3(ref w, v.octaveStep); },
(ref MessagePackReader r) =>
{
var value = default(SpanStage);
for (int i = 0, n = r.ReadArrayHeader(); i < n; i++)
{
switch (i)
{
case 0: value.frequency = r.ReadSingle(); break;
case 1: value.stretch = r.ReadSingle(); break;
case 2: value.axis = GetFloat3(ref r); break;
case 3: value.shift = GetFloat3(ref r); break;
case 4: value.offset = GetFloat3(ref r); break;
case 5: value.octaveStep = GetFloat3(ref r); break;
default: r.Skip(); break;
}
}
return value;
}),
};
var map = new Dictionary<Type, object>();

9
tests/GameCult.Caching.Tests/CultMathJsonTests.cs
View file

@ -45,6 +45,12 @@ namespace GameCult.Caching.Tests
"""{"nearest":{"x":1,"y":2,"z":3,"w":4},"edge":{"x":5,"y":6,"z":7,"w":8},"id":0.25}""");
Json(new CultPhasor(new float4(1f, 2f, 3f, 4f), new float4(5f, 6f, 7f, 8f))).Should().Be(
"""{"cos":{"x":1,"y":2,"z":3,"w":4},"sin":{"x":5,"y":6,"z":7,"w":8}}""");
Json(new SpanNoise(new float4(1f, 2f, 3f, 4f), 0.25f, 0.5f)).Should().Be(
"""{"field":{"x":1,"y":2,"z":3,"w":4},"curvature":0.25,"fold":0.5}""");
Json(new SpanPath { midpoint = new float3(1f, 2f, 3f), velocity = new float3(4f, 5f, 6f), speed = 7f, acceleration = 8f, halfLength = 9f }).Should().Be(
"""{"midpoint":{"x":1,"y":2,"z":3},"velocity":{"x":4,"y":5,"z":6},"speed":7,"acceleration":8,"halfLength":9}""");
Json(new SpanStage { frequency = 1f, stretch = 2f, axis = new float3(3f, 4f, 5f), shift = new float3(6f, 7f, 8f), offset = new float3(9f, 10f, 11f), octaveStep = new float3(12f, 13f, 14f) }).Should().Be(
"""{"frequency":1,"stretch":2,"axis":{"x":3,"y":4,"z":5},"shift":{"x":6,"y":7,"z":8},"offset":{"x":9,"y":10,"z":11},"octaveStep":{"x":12,"y":13,"z":14}}""");
}
[Test]
@ -69,6 +75,9 @@ namespace GameCult.Caching.Tests
RoundTrip(new rect(-1f, -2f, 5f, 7f));
RoundTrip(new CultCellular(new float4(1, 2, 3, 4), new float4(5, 6, 7, 8), 0.25f));
RoundTrip(new CultPhasor(new float4(1, 2, 3, 4), new float4(5, 6, 7, 8)));
RoundTrip(new SpanNoise(new float4(1, 2, 3, 4), 0.25f, 0.5f));
RoundTrip(new SpanPath { midpoint = new float3(1f, 2f, 3f), velocity = new float3(4f, 5f, 6f), speed = 7f, acceleration = 8f, halfLength = 9f });
RoundTrip(new SpanStage { frequency = 1f, stretch = 2f, axis = new float3(3f, 4f, 5f), shift = new float3(6f, 7f, 8f), offset = new float3(9f, 10f, 11f), octaveStep = new float3(12f, 13f, 14f) });
}
[Test]

21
tests/GameCult.Caching.Tests/CultMathSerializationTests.cs
View file

@ -102,6 +102,9 @@ namespace GameCult.Caching.Tests
RoundTrip(new rect(-1f, -2f, 5f, 7f));
RoundTrip(new CultCellular(new float4(1, 2, 3, 4), new float4(5, 6, 7, 8), 0.25f));
RoundTrip(new CultPhasor(new float4(1, 2, 3, 4), new float4(5, 6, 7, 8)));
RoundTrip(new SpanNoise(new float4(1, 2, 3, 4), 0.25f, 0.5f));
RoundTrip(new SpanPath { midpoint = new float3(1f, 2f, 3f), velocity = new float3(4f, 5f, 6f), speed = 7f, acceleration = 8f, halfLength = 9f });
RoundTrip(new SpanStage { frequency = 1f, stretch = 2f, axis = new float3(3f, 4f, 5f), shift = new float3(6f, 7f, 8f), offset = new float3(9f, 10f, 11f), octaveStep = new float3(12f, 13f, 14f) });
}
[Test]
@ -137,6 +140,24 @@ namespace GameCult.Caching.Tests
"92" +
"94CA3F800000CA40000000CA40400000CA40800000" +
"94CA40A00000CA40C00000CA40E00000CA41000000");
Hex(new SpanNoise(new float4(1f, 2f, 3f, 4f), 0.25f, 0.5f)).Should().Be(
"93" +
"94CA3F800000CA40000000CA40400000CA40800000" +
"CA3E800000" +
"CA3F000000");
Hex(new SpanPath { midpoint = new float3(1f, 2f, 3f), velocity = new float3(4f, 5f, 6f), speed = 7f, acceleration = 8f, halfLength = 9f }).Should().Be(
"95" +
"93CA3F800000CA40000000CA40400000" +
"93CA40800000CA40A00000CA40C00000" +
"CA40E00000CA41000000CA41100000");
Hex(new SpanStage { frequency = 1f, stretch = 2f, axis = new float3(3f, 4f, 5f), shift = new float3(6f, 7f, 8f), offset = new float3(9f, 10f, 11f), octaveStep = new float3(12f, 13f, 14f) }).Should().Be(
"96" +
"CA3F800000" +
"CA40000000" +
"93CA40400000CA40800000CA40A00000" +
"93CA40C00000CA40E00000CA41000000" +
"93CA41100000CA41200000CA41300000" +
"93CA41400000CA41500000CA41600000");
}
[Test]