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

cultmath: span stages take SpanPath and SpanStage and widen curvature and fold by float rounding allowances

Browse source 
span_fbm/span_ridged/span_stage lose the 14-argument signatures for SpanPath, SpanStage,
octaves, lacunarity, gain. Per octave S_o and e_n, e_g, e_h (af_snoise's derivation plus the
Hessian's 2^-10) widen the local Hessian tier, the gradient bound and the fold travel. Tests
compare against the exact double kernel in coordinate bands, add the pinned rank-one span,
widen the stage parameter ranges, and state the lab-form mapping. Not built or run yet.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Meta Vens 2026-10-15 20:49:30 +02:00
commit 0a1b49a231
9 changed files with 492 additions and 216 deletions
Download patch file Download diff file Expand all files Collapse all files

5
packages/cultmath/README.md
View file

@ -82,8 +82,9 @@ 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)`; and, in C# for now (shaders follow), the octave stages `span_fbm` and
`span_ridged`, which return a `SpanNoise` (`field` as `float4(gradient, value)`,
lo, hi)`; and, in C# for now (shaders follow), the octave stages `span_fbm(path, stage,
octaves, lacunarity, gain)` and `span_ridged`, over a `SpanPath` and a `SpanStage` (`SpanStage.Identity`
is the plain fbm), 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,

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

@ -809,34 +809,62 @@ configuration bound. A change to the kernels must re-run the measurements.
### Span stages
`span_fbm` and `span_ridged` (`math.Span.cs`, C# only so far) are the octave stages
of a warp chain, each returning one `SpanNoise { float4 field; float curvature; float
fold; }`. The span is the path `q(t)`, `|t| <= halfLength`, through `q` with the given
`velocity` at its midpoint, `|q'| <= speed` and `|q''| <= acceleration` along it.
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, and `stretch` must be positive and finite (unguarded, like `lacunarity`).
`field` is the octave sum and its gradient in `q` (`f_o M grad n_o`); at frequency 1,
stretch 1 and zero shift, offset and octaveStep it equals `fbm_grad` and
`ridged_grad` bit for bit. `span_ridged` sums `gain^o (1 - |n_o|)`, `ridged_grad`'s
convention; the Span Lab's `1 - 2|n|` is twice this minus the gain sum, which the
consumer forms.
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.
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 +
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| halfLength + (Hl nu |w| + Gl f_o
|M| acceleration) halfLength^2 / 2`, `w = f_o M velocity` and `Gl = min(SNOISE_LIPSCHITZ,
|grad n_mid| + 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. Rounding is not widened: the
Frobenius norm exceeds the operator norm except on a rank-one Hessian, and the pins
carry a 10% margin. `SpanStageTests` encloses the Hessian norm and the branch
departure at 64 points of seeded spans, in the local and the global regime.
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), the error is at most 0.057 of
`e_n`, 0.046 of `e_g` and 0.075 of `e_h` (the last at the smallest band, where the
absolute Hessian error is 2.3e-5 and the absolute term carries it). 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

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

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

@ -135,55 +135,78 @@ public static partial class math
// 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 frequency 1, stretch 1 and zero shift, offset and octaveStep the
// field equals fbm_grad(q, octaves, lacunarity, gain) bit for bit.
public static SpanNoise span_fbm(float3 q, float3 velocity, float speed, float acceleration, float halfLength,
float frequency, float stretch, float3 axis, float3 shift, float3 offset, float3 octaveStep, int octaves, float lacunarity, float gain) =>
span_stage(false, q, velocity, speed, acceleration, halfLength, frequency, stretch, axis, shift, offset, octaveStep, octaves, lacunarity, gain);
// 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 (the Span
// Lab's 1 - 2|n| is twice this minus the gain sum, the consumer's to form). fold bounds how far the
// sum can leave the branch it is on at the midpoint. At the identity transform the field equals
// ridged_grad bit for bit.
public static SpanNoise span_ridged(float3 q, float3 velocity, float speed, float acceleration, float halfLength,
float frequency, float stretch, float3 axis, float3 shift, float3 offset, float3 octaveStep, int octaves, float lacunarity, float gain) =>
span_stage(true, q, velocity, speed, acceleration, halfLength, frequency, stretch, axis, shift, offset, octaveStep, 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 the path q(t), |t| <= halfLength, through q at t = 0 with the given velocity there,
// |q'| <= speed and |q''| <= acceleration along it (speed at least |velocity|; all three and halfLength
// non-negative). 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.
private static float span_l1(float3 v) => abs(v.x) + abs(v.y) + abs(v.z);
// 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 + 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| + Hl nu halfLength).
// 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| halfLength + (Hl nu |w| + Gl f_o |M| acceleration) halfLength^2 / 2,
// 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.
// Rounding is not widened: the Frobenius norm exceeds the operator norm except on a rank-one
// Hessian, and the pins carry a 10% margin.
private static SpanNoise span_stage(bool ridged, float3 q, float3 velocity, float speed, float acceleration, float halfLength,
float frequency, float stretch, float3 axis, float3 shift, float3 offset, float3 octaveStep, int octaves, float lacunarity, float gain)
//
// 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. The same probe (6000 points per band, |x| from 2^-2 to 2^23) puts
// the error at most 0.057, 0.046 and 0.075 of e_n, e_g and e_h. 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(axis, axis) > 0.0f && stretch != 1.0f;
var direction = stretched ? normalize(axis) : new float3(0.0f, 0.0f, 0.0f);
var squeeze = stretched ? 1.0f / stretch - 1.0f : 0.0f;
var norm = abs(frequency) * (stretched ? max(1.0f, 1.0f / stretch) : 1.0f);
var mq = stretched ? frequency * (q + squeeze * dot(direction, q) * direction) : frequency * q;
var mv = stretched ? frequency * (velocity + squeeze * dot(direction, velocity) * direction) : frequency * velocity;
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;
@ -192,10 +215,11 @@ public static partial class math
var fold = 0.0f;
for (var i = 0; i < octaves; i++)
{
var x = scale * (mq + shift) + (offset + (float)i * octaveStep);
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 ? frequency * (g + squeeze * dot(direction, g) * direction) : frequency * g;
var mg = stretched ? stage.frequency * (g + squeeze * dot(direction, g) * direction) : stage.frequency * g;
if (ridged)
{
var s = (float)sign(n.w);
@ -210,15 +234,19 @@ public static partial class math
var weight = abs(amplitude);
var factor = abs(scale) * norm;
var nu = factor * speed;
var local = min(SNOISE_HESSIAN, span_frobenius(snoise_hessian(x)) + SNOISE_THIRD * nu * halfLength);
var nu = factor * path.speed;
var rounding = 9.5367431640625e-7f * (abs(scale) * (span_l1(mq) + span_l1(stage.shift)) + span_l1(stage.offset) + span_l1(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) + reach);
var travel = abs(dot(g, w)) * halfLength + (local * nu * length(w) + gl * factor * acceleration) * (0.5f * halfLength * 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);
}

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

@ -633,7 +633,7 @@ public sealed class NoiseBoundTests
// 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.
private static (double[][] X, double[][] P, double[] Key) Corners(double vx, double vy, double vz)
internal static (double[][] X, double[][] P, double[] Key) Corners(double vx, double vy, double vz)
{
const double cx = 1.0 / 6.0;
const double cy = 1.0 / 3.0;
@ -686,7 +686,7 @@ public sealed class NoiseBoundTests
// 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.
private static double[] HessianDouble(double vx, double vy, double vz)
internal static double[] HessianDouble(double vx, double vy, double vz)
{
var (x, p, _) = Corners(vx, vy, vz);
var hessian = new double[9];
@ -705,6 +705,23 @@ public sealed class NoiseBoundTests
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:

439
packages/cultmath/tests/CultMath.Tests/SpanStageTests.cs
View file

@ -32,48 +32,76 @@ public sealed class SpanStageTests
return normalize(v);
}
// One draw of the stage parameters: every one of them off its identity value.
private readonly record struct Stage(float Frequency, float Stretch, float3 Axis, float3 Shift, float3 Offset, float3 Step, int Octaves, float Lacunarity, float Gain)
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.
private readonly record struct Stage(SpanStage S, int Octaves, float Lacunarity, float Gain)
{
public static Stage Draw(System.Random random, int octaves) => new(
Uniform(random, 0.3f, 1.6f),
random.Next(2) == 0 ? Uniform(random, 0.2f, 0.8f) : Uniform(random, 1.25f, 5.0f),
UnitVector(random) * Uniform(random, 0.5f, 3.0f),
Vector(random, -8.0f, 8.0f),
Vector(random, -8.0f, 8.0f),
Vector(random, -3.0f, 3.0f),
octaves,
Uniform(random, 1.5f, 2.3f),
Uniform(random, 0.35f, 0.7f));
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)));
public SpanNoise Fbm(float3 q, float3 velocity, float speed, float acceleration, float halfLength) =>
span_fbm(q, velocity, speed, acceleration, halfLength, Frequency, Stretch, Axis, Shift, Offset, Step, Octaves, Lacunarity, Gain);
// 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 Ridged(float3 q, float3 velocity, float speed, float acceleration, float halfLength) =>
span_ridged(q, velocity, speed, acceleration, halfLength, Frequency, Stretch, Axis, Shift, Offset, Step, Octaves, Lacunarity, Gain);
public SpanNoise Fbm(SpanPath path) => span_fbm(path, S, Octaves, Lacunarity, Gain);
public SpanNoise Ridged(SpanPath path) => 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 = new double[] { Axis.x, Axis.y, Axis.z };
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] = Frequency * ((i == j ? 1.0 : 0.0) + (length > 0.0 ? (1.0 / Stretch - 1.0) * a[i] * a[j] / (length * length) : 0.0));
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.
public float3 Point(float3 q, int octave)
// 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 = new double[] { q.x, q.y, q.z };
var shift = new double[] { Shift.x, Shift.y, Shift.z };
var offset = new double[] { Offset.x, Offset.y, Offset.z };
var step = new double[] { Step.x, Step.y, Step.z };
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;
@ -82,28 +110,59 @@ public sealed class SpanStageTests
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) => Enumerable.Range(0, Octaves).Max(o => Magnitude(path, o));
// 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(Frequency) * Math.Max(1.0, 1.0 / 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.
private readonly record struct Span(float3 Q, float3 Velocity, float3 Acceleration, float Speed, float AccelerationBound, float HalfLength)
private 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(q, velocity, UnitVector(random) * acceleration, length(velocity) + acceleration * halfLength, acceleration, halfLength);
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) => Q + Velocity * t + Acceleration * (0.5f * t * t);
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|.
@ -121,28 +180,25 @@ public sealed class SpanStageTests
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 Hessian in q of the stage field at span point q (signed per octave for a ridged sum), by the chain
// rule from snoise_hessian: sum over octaves of weight_o f_o^2 M H_o M.
// 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.Point(q, o);
var h = snoise_hessian(x);
var weight = stage.Amplitude(o) * stage.Scale(o) * stage.Scale(o) * (ridged ? -Math.Sign(snoise(x)) : 1.0);
var hd = new double[3, 3];
for (var i = 0; i < 3; i++)
for (var j = 0; j < 3; j++)
hd[i, j] = h[i][j];
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] * hd[k, l] * m[l, j];
sum += m[i, k] * h[k * 3 + l] * m[l, j];
total[i * 3 + j] += weight * sum;
}
}
@ -168,24 +224,41 @@ public sealed class SpanStageTests
}
/// <summary>
/// At frequency 1, stretch 1 (with and without an axis) and zero shift, offset and octaveStep, the field of
/// 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 fails.
/// 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);
var none = new float3(0.0f, 0.0f, 0.0f);
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 axis = random.Next(2) == 0 ? none : UnitVector(random) * 2.0f;
var fbm = span_fbm(p, UnitVector(random), 1.0f, 0.5f, 0.25f, 1.0f, 1.0f, axis, none, none, none, octaves, lacunarity, gain);
var ridged = span_ridged(p, UnitVector(random), 1.0f, 0.5f, 0.25f, 1.0f, 1.0f, axis, none, none, none, octaves, lacunarity, gain);
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);
@ -194,10 +267,12 @@ public sealed class SpanStageTests
/// <summary>
/// The gradient is the derivative of the value through the whole transform: with every parameter off its
/// identity (frequency, stretch both ways, axis, shift, offset, octaveStep, lacunarity, gain), 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. A dropped M or lacunarity
/// factor, or M applied to the wrong side, moves it by the gradient itself.
/// 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)]
@ -207,15 +282,23 @@ public sealed class SpanStageTests
var random = new System.Random(ridged ? 0x6A1D : 0x6A1C);
var worst = 0.0;
var checkedCases = 0;
for (var k = 0; k < 400; k++)
var spans = 0;
for (var k = 0; k < 6000 && spans < 80; k++)
{
var stage = Stage.Draw(random, ridged ? random.Next(1, 3) : random.Next(1, 4));
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.3f))
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) => ridged ? stage.Ridged(p, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength) : stage.Fbm(p, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength);
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.
@ -233,18 +316,16 @@ public sealed class SpanStageTests
referenceGradient[i] += sign * stage.Amplitude(o) * stage.Scale(o) * matrix[i, j] * g[j];
}
var weight = 0.0;
for (var o = 0; o < stage.Octaves; o++)
weight += stage.Amplitude(o) * stage.Scale(o) * stage.Norm();
Assert.True(Math.Abs(referenceValue - centre.field.w) <= 5.0e-4 * (1.0 + weight), $"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]) <= 5.0e-3 * (1.0 + weight), $"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() * Math.Pow(stage.Lacunarity, stage.Octaves - 1);
var h = (float)(0.02 / fastest);
var scale = 0.0;
for (var o = 0; o < stage.Octaves; o++)
scale += stage.Amplitude(o) * stage.Scale(o) * stage.Norm();
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);
@ -276,33 +357,46 @@ public sealed class SpanStageTests
var random = new System.Random(0x57E7C4);
for (var k = 0; k < 60; k++)
{
var stage = Stage.Draw(random, 3);
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 { Stretch = 1.0f };
var noAxis = stage with { Axis = new float3(0.0f, 0.0f, 0.0f), Stretch = 3.0f };
var reference = plain.Ridged(q, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength);
var isotropic = (plain with { Axis = new float3(0.0f, 0.0f, 0.0f) }).Ridged(q, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength);
var viaNoAxis = noAxis.Ridged(q, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength);
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);
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");
var unit = stage with { Axis = new float3(0.0f, 1.0f, 0.0f) };
var scaled = stage with { Axis = new float3(0.0f, 2.0f, 0.0f) };
var a = unit.Fbm(q, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength);
var b = scaled.Fbm(q, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength);
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(q, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength);
Assert.False(SameBits(stretched.field, (unit with { Stretch = 1.0f }).Fbm(q, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength).field), "a stretch of " + unit.Stretch + " did not change the field");
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 stage parameters and per-octave bookkeeping the curvature checks share.
// 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;
private static Stage Banded(System.Random random, float offsetScale, int octaves)
{
var stage = Stage.Draw(random, octaves);
if (offsetScale > 0.0f)
stage = stage with { S = 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.HalfLength * (2.0f * s / (samples - 1) - 1.0f);
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);
@ -311,13 +405,15 @@ public sealed class SpanStageTests
}
/// <summary>
/// span-bounds-enclose, curvature: over seeded spans the operator norm of the field's q-Hessian (chain rule
/// from snoise_hessian, 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); at halfLength 0 curvature is the local tier, the Frobenius norm of the
/// midpoint Hessians scaled by amplitude, f_o^2 and |M|^2, to within float rounding; and a span through
/// MeasureHessian's witness stays enclosed in 200 directions, which the Frobenius norm at the midpoint alone
/// does not guarantee.
/// 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)]
@ -327,28 +423,34 @@ public sealed class SpanStageTests
var random = new System.Random(ridged ? 0xC0DE2 : 0xC0DE1);
var tightest = 1.0;
var checkedPoints = 0;
foreach (var halfLength in new[] { 0.002f, 0.02f, 0.2f, 1.0f, 8.0f })
for (var k = 0; k < 24; k++)
{
var stage = Stage.Draw(random, random.Next(1, 4));
var span = Span.Draw(random, Vector(random, -30.0f, 30.0f), halfLength, Uniform(random, 0.0f, 2.0f));
var result = ridged ? stage.Ridged(span.Q, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength) : stage.Fbm(span.Q, span.Velocity, span.Speed, span.AccelerationBound, span.HalfLength);
CheckCurvature(stage, ridged, span, result, 64, ref tightest, ref checkedPoints);
}
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);
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);
}
// halfLength 0: the local tier exactly.
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.Draw(random, random.Next(1, 4));
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(q, span.Velocity, span.Speed, span.AccelerationBound, 0.0f) : stage.Fbm(q, span.Velocity, span.Speed, span.AccelerationBound, 0.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]));
expected += stage.Amplitude(o) * stage.Scale(o) * stage.Scale(o) * stage.Norm() * stage.Norm() * Math.Min(SNOISE_HESSIAN, frobenius);
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));
@ -357,27 +459,50 @@ public sealed class SpanStageTests
// 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(1.0f, 1.0f, new float3(0.0f, 0.0f, 0.0f), new float3(0.0f, 0.0f, 0.0f), new float3(0.0f, 0.0f, 0.0f), new float3(0.0f, 0.0f, 0.0f), 1, 2.0f, 0.5f);
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(witness + u * distance, -u, new float3(0.0f, 0.0f, 0.0f), 1.0f, 0.0f, 0.6f);
var result = ridged ? identity.Ridged(span.Q, span.Velocity, 1.0f, 0.0f, 0.6f) : identity.Fbm(span.Q, span.Velocity, 1.0f, 0.0f, 0.6f);
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, least relative slack {tightest:R}, witness norm {through:R}");
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");
}
/// <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; fold is positive on spans that do cross zero (and the departure there is not trivially
/// zero), exactly 0 when no octave can reach zero and for span_fbm, and 0 at halfLength 0.
/// 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()
@ -385,62 +510,77 @@ public sealed class SpanStageTests
var random = new System.Random(0xF01D);
var crossing = 0;
var largest = 0.0;
var tightest = 1.0;
for (var k = 0; k < 400; k++)
{
var stage = Stage.Draw(random, random.Next(1, 4));
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 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 best = q;
var smallest = float.MaxValue;
for (var trial = 0; trial < 40; trial++)
var stage = Banded(random, offsetScale, 3);
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 candidate = Vector(random, -30.0f, 30.0f);
var value = abs(snoise(stage.Point(candidate, 0)));
if (value < smallest)
(smallest, best) = (value, candidate);
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;
}
q = best;
}
var span = Span.Draw(random, q, halfLength, Uniform(random, 0.0f, 2.0f));
var fold = stage.Ridged(q, span.Velocity, span.Speed, span.AccelerationBound, halfLength).fold;
Assert.Equal(0.0f, stage.Fbm(q, span.Velocity, span.Speed, span.AccelerationBound, halfLength).fold);
Assert.Equal(0.0f, stage.Ridged(q, span.Velocity, span.Speed, span.AccelerationBound, 0.0f).fold);
var signs = Enumerable.Range(0, stage.Octaves).Select(o => Math.Sign(snoise(stage.Point(q, o)))).ToArray();
var peak = 0.0;
for (var s = 0; s < 64; s++)
{
var t = halfLength * (2.0f * s / 63 - 1.0f);
var point = span.At(t);
var departure = 0.0;
for (var o = 0; o < stage.Octaves; o++)
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 signs = Enumerable.Range(0, stage.Octaves).Select(o => Math.Sign(stage.OneOctave(o).Fbm(PathAt(q)).field.w)).ToArray();
var peakFloat = 0.0;
var peakExact = 0.0;
for (var s = 0; s < 64; s++)
{
var n = (double)snoise(stage.Point(point, o));
departure += stage.Amplitude(o) * (signs[o] * n - Math.Abs(n));
var t = halfLength * (2.0f * s / 63 - 1.0f);
var point = span.At(t);
var departureFloat = 0.0;
var departureExact = 0.0;
for (var o = 0; o < stage.Octaves; o++)
{
var octave = stage.OneOctave(o);
var nFloat = (double)octave.Fbm(PathAt(point)).field.w;
var absFloat = 1.0 - octave.Ridged(PathAt(point)).field.w;
departureFloat += stage.Amplitude(o) * (signs[o] * nFloat - absFloat);
var x = stage.PointD(point, o);
var nExact = NoiseBoundTests.ValueDouble(x[0], x[1], x[2]);
departureExact += stage.Amplitude(o) * (signs[o] * nExact - Math.Abs(nExact));
}
peakFloat = Math.Max(peakFloat, Math.Abs(departureFloat));
peakExact = Math.Max(peakExact, Math.Abs(departureExact));
Assert.True(Math.Abs(departureFloat) <= fold * (1.0 + 1.0e-3) + 1.0e-6, $"stage {stage}, span {span}: the float branch departure at t = {t:R} is {departureFloat:R}, above fold {fold:R}");
Assert.True(Math.Abs(departureExact) <= fold * (1.0 + 1.0e-3) + 1.0e-6, $"stage {stage}, span {span}: the exact branch departure at t = {t:R} is {departureExact:R}, above fold {fold:R}");
}
peak = Math.Max(peak, Math.Abs(departure));
Assert.True(Math.Abs(departure) <= fold * (1.0 + 1.0e-3), $"stage {stage}, span {span}: the branch departure at t = {t:R} is {departure: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);
}
if (fold > 0.0)
crossing++;
largest = Math.Max(largest, peak);
if (fold > 0.0 && peak > 0.0)
tightest = Math.Min(tightest, peak / fold);
}
Assert.True(crossing >= 100, $"only {crossing} of 400 spans had a positive 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, largest departure {largest:R}, tightest departure / fold {tightest:R}");
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(1.0f, 1.0f, new float3(0.0f, 0.0f, 0.0f), new float3(0.0f, 0.0f, 0.0f), new float3(0.0f, 0.0f, 0.0f), new float3(0.0f, 0.0f, 0.0f), 1, 2.0f, 0.5f);
var calm = new Stage(SpanStage.Identity, 1, 2.0f, 0.5f);
var found = 0;
for (var k = 0; k < 400 && found < 20; k++)
{
@ -448,7 +588,7 @@ public sealed class SpanStageTests
if (abs(snoise(q)) < 0.6f)
continue;
found++;
Assert.Equal(0.0f, calm.Ridged(q, new float3(1.0f, 0.0f, 0.0f), 1.0f, 0.0f, 0.01f).fold);
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");
@ -464,8 +604,7 @@ public sealed class SpanStageTests
public void SpanFoldCoversAccelerationAlone()
{
var random = new System.Random(0xACCE1);
var identity = new Stage(1.0f, 1.0f, new float3(0.0f, 0.0f, 0.0f), new float3(0.0f, 0.0f, 0.0f), new float3(0.0f, 0.0f, 0.0f), new float3(0.0f, 0.0f, 0.0f), 1, 2.0f, 0.5f);
var still = new float3(0.0f, 0.0f, 0.0f);
var identity = new Stage(SpanStage.Identity, 1, 2.0f, 0.5f);
var cases = 0;
var crossings = 0;
var tightest = double.MaxValue;
@ -480,8 +619,8 @@ public sealed class SpanStageTests
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(q, still, c, acceleration * halfLength, acceleration, halfLength);
var fold = identity.Ridged(q, still, span.Speed, acceleration, halfLength).fold;
var span = new Span(new SpanPath { midpoint = q, velocity = None, speed = acceleration * halfLength, acceleration = acceleration, halfLength = halfLength }, c);
var fold = identity.Ridged(span.Path).fold;
var departure = 0.0;
for (var s = 0; s < 64; s++)
{