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-affine-forms: affine forms for noise bounds, in C#, HLSL and GLSL

Browse source 
Reduced affine arithmetic (centre, one shared symbol, one error term) for CultMath:
af_* ops, af_snoise / af_fbm with gradient-plus-curvature forms, af_frustum_*, the
2D iv_snoise_ball, and the centred warp bound. Mirrored bit for bit in
CultMath.Affine.hlsl (and the Unity copy), lowered into CultMath.glsl, with fixture
families that check enclosure on device, including small-reach and C#-centre legs.
Scenario (c) is a deep well an order of magnitude past the fog band; affine beats the
centred interval there (3.57x / 3.46x vs 3.20x / 3.12x), asserted, not skipped.
Soul held cut-affine-forms (s1-s3) and cut-affine-shaders (s1-s3).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Meta Vens 2026-10-04 12:57:11 +02:00
commit 156782da0e
21 changed files with 7497 additions and 114 deletions
Download patch file Download diff file Expand all files Collapse all files

3
packages/cultmath/README.md
View file

@ -34,7 +34,8 @@ shader semantics.
parity. HLSL already owns `float2`, `float3`, and `float4`; the include
exposes `cultmath_*` functions for shared semantics HLSL does not provide.
It includes `shaders/CultMath.Phacelle.hlsl` (MPL-2.0) and
`shaders/CultMath.Interval.hlsl` from the same directory, so ship all three.
`shaders/CultMath.Interval.hlsl` and `shaders/CultMath.Affine.hlsl` from the same
directory, so ship all four.
- interval arithmetic for culling empty space with a proof: `iv_*` functions
over `float2(lo, hi)` (`iv_add`, `iv_mul`, `iv_exp`, `iv_smoothstep`, ...)
and `iv_snoise_ball`/`iv_fbm_ball`, which bound `snoise` and `fbm_grad` over a

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

@ -222,6 +222,18 @@ lowering that inflates the ball fails too. Only those two bounds fail the
check; the largest ulp distance is reported as for any ulp-bounded family.
`GoldenFixtureMatchesCSharp` applies the same check to C#'s own results.
Every `af_*` function is an ulp-bounded family, and `af_snoise`, `af_fbm` and `af_frustum_ball` also carry a two-sided `check`
of their own, applied the same way: the device's own `snoise` or fBm at 75 sampled points of the case's set must lie in the device's
form with no tolerance (a device whose gradient error exceeds the allowance fails), its width must not pass the Lipschitz form's
allowance, and the device's `af_frustum_ball` radius must enclose the sampled ray points and stay within its rounding widening.
The reference is the device's own noise, which is right for culling soundness but would pass a device whose noise is wrong
everywhere, so two appended families (`af_snoise` and `af_fbm`, 256 cases each) add the centre leg: the fixture's own `x0`
(C#'s value at the centre) must lie in the device form's `[x0 - e, x0 + e]`, for every `af_snoise` and `af_fbm` family. They
also draw reach `|axis| + radius` from 1e-6 to 1e-3, the short slices near a camera, where the float32 allowance is most of `e`:
the main families' reach is at least about 1e-3, where the second-order remainder is far larger and a device that drops the
allowance still passes. `GoldenFixtureMatchesCSharp` requires each of the two to hold cases that fail on C#'s own results once
the allowance is removed.
`GoldenFixtureMatchesCSharp` evaluates C# on the committed fixture's own
arguments. Arguments are drawn only when the fixture is regenerated
(`CULTMATH_WRITE_GLSL=1`), because the draw itself calls `MathF.Pow`, whose
@ -586,9 +598,9 @@ over 2,000 balls x 64 points are the defence. A change to the `snoise` kernel
must re-run `MeasureLipschitz` and re-pin `L` and the witness.
What the ball bounds guarantee: enclosure of `snoise` and `fbm_grad(...).w` over
the ball, in float32. What they do not: anything about tapes, pruning, choice
tracking or affine forms; those are later steps of
`docs/cultmath-tape-target.md` (CultLib root). How loose they are is a number:
the ball, in float32. What they do not: anything about tapes, pruning or choice tracking. Correlation
between the operands of one expression is what "Affine forms" below adds.
How loose the ball bounds are is a number:
`NoiseBoundTests.TightnessReport` prints the mean interval width over the
sampled range: 5.08 for `iv_snoise_ball` at `r = 0.05`, 1.96 at 0.25 and 1.28
at 1.0; 6.64 and 2.82 for 4-octave `iv_fbm_ball` at 0.05 and 0.25.
@ -616,30 +628,43 @@ printed. The flow is not counted. At `N = 8`, in Aetheria's units on its
- (a) Height fog, camera above the safety band, rays level and up: the whole
grid is one envelope probe per tile, against 256 envelope evaluations per ray,
16384x. Every cell is oracle-empty and the pre-pass proves them all.
- (b) Inside the fog: 1.15x, and 0.99x with the warp. Most cells are fog, and
- (b) Inside the fog: 1.16x, and 1.11x with the warp. Most cells are fog, and
the pre-pass probes each one.
- (c) The zone bowl and four wells, the camera gazing across. With no warp:
2.74x, oracle ceiling 7.26x, efficiency 0.378, cull fraction 0.952, probe
overhead 0.961. With the warp `D = 60`: 1.39x, oracle ceiling 7.17x,
efficiency 0.194, cull fraction 0.842, probe overhead 0.993. With the warp the
pre-pass cuts envelope evaluations from 177.46 to 33.49 per ray, but the
oracle march costs 0.194 of the tile march, so about four fifths of the tile
march's cost goes to cells no ray samples nonzero: the warp bound swamps the
noise ball, so the ball cannot prove those cells empty. That gap belongs to
the tape target's affine forms, not to a better pre-pass.
- (d) r1's uniform slab, which has no envelope: 1.23x, and 1.48x with no warp.
- (c) The deep well (ruling `wells-scale-order-of-magnitude`): the zone bowl
plus one well at its centre, `PowerPulse` exponent 2, scale 2220, depth 600.
The band `G = 60` is where only the noise decides density (`s` within
`F -+ 1.5 A` with the fade at 1); depth and half-depth radius are each
`10 G`. The camera's `xz` is uniform within 300 of the centre, at
`y = F - h(xz) / 2`, so it looks across the bowl from inside it. Both rows run
on the same seeds. With no warp, the centred interval takes 3.20x (oracle
ceiling 5.24x, efficiency 0.611, cull fraction 0.988, probe overhead 0.940)
and the affine bound 3.57x (efficiency 0.682, cull 0.993, overhead 0.905). With
the warp `D = 60`: 3.12x (ceiling 5.31x, efficiency 0.588, cull 0.987, overhead
0.929) against 3.46x (efficiency 0.651, cull 0.993, overhead 0.880), so affine
is 1.12x and 1.11x cheaper. The range mask (the oracle cells plus every cell
where any of 5 x 5 slopes x 9 depths has density above zero) gives a range
ceiling of 0.859 (0.848 with the warp), the cost a bound that proves exactly
the range-empty cells would reach against the oracle's; the efficiencies
against it are 0.794 and 0.712 (warp 0), 0.768 and 0.694 (warp `D`).
- (c-band) The shallow wells, four mass-derived wells in the zone bowl with
the camera gazing across, whose geometry is band-dominated, printed only.
Centred interval 3.01x (efficiency 0.414, cull 0.960, overhead 0.958), affine
3.95x (0.544, 0.980, 0.901); with the warp 2.75x (0.383, 0.954, 0.944) against
3.54x (0.494, 0.976, 0.865).
- (d) r1's uniform slab, which has no envelope: 1.38x with the warp 0.05, and
1.67x with no warp.
- (e) The void at the shipped point (ruling `operator-shipped-void`): `Rh = 198`,
`S = 50`, `Rc = Rh + S = 248`, `e = -ln 1.5 / ln(1 - (Rh / Rc)^2)`, the
shipped camera, 200 tiles. The fixed-step reference takes 36.34 steps per
pixel. The same cells masked by the pre-pass take 20.13, with identical
transmittance. The footprint-aware march takes 13.54: 2.68x fewer. `snoise`
per pixel goes from 32.98 to 26.69, and combined cost from 40.25 to 29.52 with
pixel. The same cells masked by the pre-pass take 20.08, with identical
transmittance. The footprint-aware march takes 13.50: 2.69x fewer. `snoise`
per pixel goes from 32.98 to 26.70, and combined cost from 40.25 to 29.54 with
the probes, 1.36x. Every one of the 200 tiles proves a body. Before the body
both marches are measured against a converged one (midpoint quadrature at
`h = 0.5` over every cell, trusting no pre-pass): the relative optical-depth
error at the body start is at most 0.0608 for the fixed-step march and 0.08639
for the footprint march, 1.42x. Over the grid of hollow radius, camera offset
and ramp width the steps ratio runs from 1.57x (`Rh = 100`, low, `S = 150`) to
and ramp width the steps ratio runs from 1.59x (`Rh = 100`, low, `S = 150`) to
27.45x (`Rh = 1000`, centred, `S = 10`), the cost ratio from 1.26x to 3.89x,
and the depth-error ratio from 0.85 to 4.72 (`Rh = 200`, centred, `S = 150`,
where the errors are 0.01504 and 0.003188). LOD-far takes 14.19x fewer steps;
@ -669,24 +694,144 @@ interval argument, so the composed interval encloses both fields.
The contracts, at `N = 8`, the lower of warp 0 and `D`:
- (a) at least 2x combined cost;
- (c) (ruling `wells-overhead-plus-floor`): probe overhead, the tile march over
the tile march plus its probes, at least 0.90; and, hard, the pre-pass proves
at least half of the oracle-empty cells empty, so a pre-pass that skips
nothing fails. The oracle ceiling, the efficiency against it and the cull
fraction are printed, not asserted;
- (c) on the deep well (rulings `affine-wells-contract-cheaper-than-interval`
and `wells-overhead-plus-floor`): the affine bound's combined cost ratio is
above the centred interval's at warp 0 and at warp `D`; and, hard, the
pre-pass proves at least half of the oracle-empty cells empty under both
bounds, so a pre-pass that skips nothing fails. The probe overhead, the
efficiencies against the oracle and the range ceiling, the ceilings and the
cull fraction are printed, not asserted; (c-band) is printed only;
- (e) at the shipped void: the masked fixed-step march agrees with the reference
exactly; the footprint march's optical depth at the body start is off the
converged march's by at most twice the fixed-step march's error; and it takes
at most half the reference's steps per pixel.
This run meets every contract: (a) 16384x, (c) probe overhead 0.961 and cull
fraction 0.842, (e) 2.68x fewer steps and a depth error 1.42x the fixed-step
march's. The fields, the grids and the 0.2 weight are not tuned toward the
contracts.
This run meets every contract: (a) 16384x, (c) affine 3.57x and 3.46x against
the centred interval's 3.20x and 3.12x, with a cull fraction of 0.987 at the
lowest, (e) 2.69x fewer steps and a depth error 1.42x the fixed-step march's.
The fields, the grids, the 0.2 envelope weight and the 2.0 gradient weight are
not tuned toward the contracts. An affine probe is one `snoise_grad`, counted at
2.0 `snoise` evaluations (an estimate, printed beside the count).
Consumers: the gamecult.org ground shader (through the GLSL lowering) and
Aetheria's nebula raymarch at its CultMath pin bump.
## Affine forms
`math.Affine.cs` carries reduced affine forms with one shared symbol: a
`float3(x0, a, e)` is the set `x0 + a eps + e delta` with `eps` and `delta` in
`[-1, 1]` and `e >= 0` (Messine's AF1, as Gamito and Maddock use it for gradient
noise, in one register as an interval is a `float2`). `eps` is shared by every
form over one region (for a screen-tile slice, the depth parameter, `z = z_m +
h eps`), so `af_sub(x, x)` has `a = 0`; `delta` is private to the form, and
every term that cannot be tracked is condensed into `e`. A slice is a
`float4(centre, radius)` with a `float3` axis. No op takes a compound type, a
reference parameter or a swizzle write, so each mirrors into the HLSL common
subset.
The ops: `af_point`, `af_symbol`, `af_from_iv`, `af_range` (the bridge back to
`iv_*`), `af_add`, `af_sub`, `af_neg`, `af_scale`, `af_add_iv`, `af_mul` (the
product `|a_x a_y|` and the cross-error terms condensed into `e`, so there is
no `af_sqr`), `af_snoise` and `af_fbm`, and the slice pair `af_frustum_axis` and
`af_frustum_ball`.
Rounding. Every op returns a form whose set at each `eps` contains the float32
value of its function at every point of its operands' sets at that `eps`. Each
op folds its own float32 rounding into `e` as `S 2^-20`, `S` the sum of the
magnitudes its roundings are relative to (named in the op's comment): at most
`8u S` is needed with `u = 2^-24`, and the factor of two over that is the
allowance for a device whose `+` and `*` are faithful rather than correctly
rounded. `AffineTests.EveryOpEnclosesItsPoints` evaluates the form's set in
double and the function in float32 with no tolerance. Products and `af_from_iv` also
add `2^-126`, the most an underflow to a subnormal can lose (`TinyOperandsEncloseTheirPoints`).
`af_from_iv` halves its endpoints before it sums and differences them, and `af_range`
scales its terms before it sums them, so endpoints near the float maximum do not overflow
in the middle of the op; a form that overflows anyway gives `af_range` infinite endpoints,
never NaN (`FromIntervalAndRangeHoldAtTheExtremes`). `af_range` finds a NaN endpoint by its bit pattern, because
`lo == lo` is a compare dxc marks `fast` and deletes under its default flags. `AbsorptionsArePinned` and
`SnoiseFormCarriesItsDocumentedWidth` pin each absorption and `af_snoise`'s float32 allowance both
ways: the enclosure tests alone cannot see them on a CPU, where float32 is exact.
`af_snoise(centre, axis, radius)` encloses `snoise` over `centre + axis eps + w`,
`|w| <= radius`, from one `snoise_grad`. It returns the narrower of the
Lipschitz form `(n, 0, L R)` and the centred form `(n, g . axis, |g| radius +
M R^2 / 2)` (`R = |axis| + radius`, `L = SNOISE_LIPSCHITZ`, `M =
SNOISE_HESSIAN`), so its total width is never more than the Lipschitz form's.
The centred form is Taylor's theorem, which is why it beats the ball when the
radius is small and the axis carries the variation. `af_fbm` compounds octaves
as `iv_fbm_ball` does. The float32 evaluation allowance is derived in
`af_snoise`'s comment.
`SNOISE_HESSIAN = 56.050385` is 1.10 times the refined maximum of the spectral
norm of the Hessian of `snoise(float3)`, measured by
`NoiseBoundTests.MeasureHessian` over 1e6 seeded points (central differences of
`snoise_grad`) and refined by ascent from the 1e4 largest: sampled 50.1366,
refined 50.9549, witness (-245.19885, -53.01323, -91.47363).
`HessianConstantPinsSampledCurvature` fails on a 1% move.
`SNOISE2_LIPSCHITZ = 8.117002` is the same for the gradient of `snoise(float2)`
(`MeasureLipschitz2`: sampled 7.2575, refined 7.3791, witness (-67.12454,
-129.53296), `|snoise2|` at most 0.99964), and it backs `iv_snoise_ball(float2,
float)`. Both are empirical with a margin, not proofs. A random ball cannot see a
curvature constant 20% low, so `SNOISE_HESSIAN` is defended by
`HessianConstantPinsSampledCurvature` and `SnoiseFormEnclosesTheCurvatureWitness`
(`af_snoise` at the point where the curvature is realised, along 4,000 directions: a
constant of 40 misses there). The margin over the function itself: a scan of the whole
289^3 hash period found a realised maximum of 50.58, and bounding each corner's gradient
independently over every hash configuration gives 52.22, so 56.05 sits 7% above the
configuration bound. A change to the kernels must re-run the measurements.
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`:
every ray point at depth `z` is `centre + axis eps(z) + w`. Its derivation
and the clamp residual are in the comment. The caller rotates, translates and
scales the slice into the noise's frame, and owns the rounding of that.
The centred warp. A flow-warped sample reads `snoise` at `p + flow(p) shift`. The
caller does not enlarge the ball by the whole displacement `|flow| shift`; it
moves the centre by the flow at the centre, `c + flow(c) shift`, and widens the
radius by the flow's variation across the ball, `r + |shift| rho` with `rho =
FlowAmplitude (SNOISE2_LIPSCHITZ FlowFrequency r + a) sqrt 2` (the clamp to the unit
disc is nonexpansive), where `a = 2^-14 + L2 2^-20 (|q|_1 + r FlowFrequency)` is the
float32 allowance of the flow's two `snoise(float2)` reads (the point `q` plus its offset rounds at the
magnitude of `q`, and a float32 flow at two points of a ball differs by more than the exact
function can). CultMath owns neither the flow nor the envelope: the warp's
centre and variation are the caller's, and a consumer adds them to the slice's
`warpVariation` or to the form's radius. With two phases that cross-fade, each
phase takes its own centre; one ball serves when the phases' centres coincide
(warp 0). `NoiseBoundTests.CentredWarpStaysEnclosed` checks 2,000 balls, radius 1e-3 to
500, and fails when `rho` is dropped or halved; `AffineWarpStaysEnclosed` checks
3,000 tile slices of `FogField.Slice`, which returns the final `af_snoise` arguments
(each phase's moved centre and grown radius, the detail's unwarped set, the shared axis),
so the affine bound wires nothing itself: the reach over which the flow is varied is
`|axis| + radius`, each phase carries its own growth, and the gate reads each argument's
reach. It fails when the reach is the ball's radius or the larger of the two, a phase's
growth is dropped or taken from the other phase, the axis or the detail radius is scaled,
or `rho` is halved. `AffineCompositionTightensItsRanges` sees the shared-symbol
composition itself: the bound is exactly the intersection of the affine arithmetic's piece
and the interval composition over the forms' ranges, and the affine piece is the narrower
of the two in at least 35% of the open slices (the piece alone is wider in total).
The deep well is the representative Aetheria case: the camera inside a sun's
bowl, looking across. `DeepWellIsTenBandsDeepAndWide` computes the band `G` from
`FogField`'s own floor offset, amplitude and noise range (60) and requires the
depth (633 at the centre against the rim) and the half-depth radius to be at
least `10 G`, every drawn camera sitting inside the well below the rim; a well of
depth 150 fails. The shallow wells stay as the printed shallow-well, band-dominated row (c-band).
The test field's affine bound takes depth as `af_symbol(z_m, h)`, builds `s`
affine in it with the height interval, the footprint and the slope-depth
product in `e` through `af_add_iv`, takes the noise through `af_frustum_*` and
`af_snoise` per phase and for the detail octave, fades by `af_from_iv` and
`af_mul`, reads the range with `af_range`, and intersects with the centred
interval's bound (both sound). `FogField.Bound` makes no direct `snoise` call.
The cull-map render (`DeepWellCullMap`, explicit) writes a PNG when
`CULTMATH_WRITE_CULLMAP` names a directory: one deep-well camera at warp `D`,
960 x 540 in 8 x 8 tiles, with panels for dense cells per ray under the centred
interval, under the affine bound and under the oracle mask (one logarithmic
scale) and the transmittance. Its run: 4.23, 3.60 and 2.02 cells per ray on
average, 44 the most, mean transmittance 0.278.
## Rules
- Keep type names HLSL-shaped: `float2`, `float3`, `float4`, `float3x3`, and `math`.

4
packages/cultmath/scripts/build-unity-package.ps1
View file

@ -84,7 +84,9 @@ $requiredTemplateFiles = @(
"Shaders\CultMath.Phacelle.hlsl",
"Shaders\CultMath.Phacelle.hlsl.meta",
"Shaders\CultMath.Interval.hlsl",
"Shaders\CultMath.Interval.hlsl.meta"
"Shaders\CultMath.Interval.hlsl.meta",
"Shaders\CultMath.Affine.hlsl",
"Shaders\CultMath.Affine.hlsl.meta"
)
foreach ($relativePath in $requiredTemplateFiles) {
if (-not (Test-Path -LiteralPath (Join-Path $templateRoot $relativePath) -PathType Leaf)) {

132
packages/cultmath/shaders/CultMath.Affine.hlsl Normal file
View file

@ -0,0 +1,132 @@
#ifndef CULTMATH_AFFINE_HLSL
#define CULTMATH_AFFINE_HLSL
// Reduced affine forms with one shared symbol over float3(x0, a, e): the set { x0 + a eps + e delta :
// delta in [-1, 1] } at the region's shared eps in [-1, 1], e >= 0. Every cultmath_af_* function mirrors
// math.Affine.cs bit for bit; see that file and docs/design.md ("Affine forms") for each op's enclosure
// and the float32 rounding allowance it folds into e. Included by CultMath.hlsl after
// CultMath.Interval.hlsl, whose CULTMATH_SNOISE_LIPSCHITZ it shares, and after its own
// definition of cultmath_snoise_grad.
// An upper bound on the spectral norm of the Hessian of cultmath_snoise(float3); provenance in
// math.Affine.cs (SNOISE_HESSIAN).
static const float CULTMATH_SNOISE_HESSIAN = 56.050385;
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); }
float3 cultmath_af_from_iv(float2 interval)
{
float x0 = interval.x * 0.5 + interval.y * 0.5;
float e = interval.y * 0.5 - interval.x * 0.5;
return float3(x0, 0.0, e + (abs(x0) + e) * 9.5367431640625e-7 + 1.17549435e-38);
}
// An endpoint that would be NaN is the infinity on its side. NaN is found by its bit pattern: dxc marks
// compares fast, so a test that needs NaN to make `lo == lo` false is deleted under its default flags.
float2 cultmath_af_range(float3 x)
{
float r = abs(x.y) + x.z;
float w = abs(x.x) * 9.5367431640625e-7 + r * 9.5367431640625e-7;
float lo = x.x - r - w;
float hi = x.x + r + w;
return float2((asuint(lo) & 0x7fffffffu) > 0x7f800000u ? asfloat(0xff800000u) : lo,
(asuint(hi) & 0x7fffffffu) > 0x7f800000u ? asfloat(0x7f800000u) : hi);
}
float3 cultmath_af_add(float3 x, float3 y)
{
float x0 = x.x + y.x;
float a = x.y + y.y;
float e = x.z + y.z;
return float3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7);
}
float3 cultmath_af_sub(float3 x, float3 y)
{
float x0 = x.x - y.x;
float a = x.y - y.y;
float e = x.z + y.z;
return float3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7);
}
float3 cultmath_af_neg(float3 x) { return float3(-x.x, -x.y, x.z); }
float3 cultmath_af_scale(float3 x, float s)
{
float x0 = x.x * s;
float a = x.y * s;
float e = x.z * abs(s);
return float3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7 + 1.17549435e-38);
}
float3 cultmath_af_add_iv(float3 x, float2 interval) { return cultmath_af_add(x, cultmath_af_from_iv(interval)); }
float3 cultmath_af_mul(float3 x, float3 y)
{
float x0 = x.x * y.x;
float p = x.x * y.y;
float q = y.x * x.y;
float a = p + q;
float e = abs(x.x) * y.z + abs(y.x) * x.z + (abs(x.y) + x.z) * (abs(y.y) + y.z);
return float3(x0, a, e + (abs(x0) + abs(p) + abs(q) + e) * 9.5367431640625e-7 + 1.17549435e-38);
}
// The narrower of the Lipschitz form and the centred (Taylor) form over every p = centre + axis eps + w,
// |w| <= radius, from one cultmath_snoise_grad at the centre; the allowance for the float32 evaluation is
// in math.Affine.cs.
float3 cultmath_af_snoise(float3 centre, float3 axis, float radius)
{
float4 g = cultmath_snoise_grad(centre);
float3 gradient = float3(g.x, g.y, g.z);
float reach = length(axis) + radius;
float lipschitz = CULTMATH_SNOISE_LIPSCHITZ * reach;
float a = dot(gradient, axis);
float remainder = length(gradient) * radius + 0.5 * CULTMATH_SNOISE_HESSIAN * reach * reach;
float allowance = 6.103515625e-5 + 2.44140625e-4 * reach
+ CULTMATH_SNOISE_LIPSCHITZ * 9.5367431640625e-7 * (abs(centre.x) + abs(centre.y) + abs(centre.z) + reach);
bool centred = abs(a) + remainder < lipschitz;
float slope = centred ? a : 0.0;
float e = (centred ? remainder : lipschitz) + allowance;
return float3(g.w, slope, e + (abs(g.w) + abs(slope) + e) * 9.5367431640625e-7);
}
// Octave sum of cultmath_af_snoise, compounding frequency and amplitude as the value-and-gradient fBm does.
float3 cultmath_af_fbm(float3 centre, float3 axis, float radius, int octaves, float lacunarity, float gain)
{
octaves = clamp(octaves, 0, 16);
float amplitude = 1.0;
float frequency = 1.0;
float3 sum = float3(0.0, 0.0, 0.0);
for (int i = 0; i < octaves; i++)
{
sum = cultmath_af_add(sum, cultmath_af_scale(cultmath_af_snoise(centre * frequency, axis * frequency, radius * abs(frequency)), amplitude));
frequency *= lacunarity;
amplitude *= gain;
}
return sum;
}
// The axis of a screen tile's slice over the depth segment [z0, z1]; camera frame of
// cultmath_iv_frustum_ball.
float3 cultmath_af_frustum_axis(float2 centreSlope, float z0, float z1)
{
float h = (z1 - z0) * 0.5;
return float3(centreSlope.x * h, centreSlope.y * h, h);
}
// The affine slice of a screen tile over [z0, z1]: centre xyz, radius w. Derivation in math.Affine.cs.
float4 cultmath_af_frustum_ball(float2 centreSlope, float z0, float z1, float footprintPerDepth, float warpVariation)
{
float zm = (z0 + z1) * 0.5;
float h = (z1 - z0) * 0.5;
float cx = centreSlope.x * zm;
float cy = centreSlope.y * zm;
float spread = (abs(centreSlope.x) + abs(centreSlope.y) + 1.0) * h;
float radius = z1 * footprintPerDepth + warpVariation;
radius += (abs(cx) + abs(cy) + abs(zm) + spread + radius) * 9.5367431640625e-7;
return float4(cx, cy, zm, radius);
}
#endif

7
packages/cultmath/shaders/CultMath.Affine.hlsl.meta Normal file
View file

@ -0,0 +1,7 @@
fileFormatVersion: 2
guid: 26f28fce38e14eed8b9b5dbd758e4ffc
ShaderIncludeImporter:
externalObjects: {}
userData:
assetBundleName:
assetBundleVariant:

11
packages/cultmath/shaders/CultMath.Interval.hlsl
View file

@ -9,6 +9,9 @@
// Lipschitz constant of cultmath_snoise(float3); provenance in math.Interval.cs (SNOISE_LIPSCHITZ).
static const float CULTMATH_SNOISE_LIPSCHITZ = 7.9640074;
// Lipschitz constant of cultmath_snoise(float2); provenance in math.Interval.cs (SNOISE2_LIPSCHITZ).
static const float CULTMATH_SNOISE2_LIPSCHITZ = 8.117002;
float2 cultmath_iv_point(float value) { return float2(value, value); }
float2 cultmath_iv_add(float2 a, float2 b) { return float2(a.x + b.x, a.y + b.y); }
float2 cultmath_iv_sub(float2 a, float2 b) { return float2(a.x - b.y, a.y - b.x); }
@ -87,6 +90,14 @@ float2 cultmath_iv_snoise_ball(float3 centre, float radius)
return float2(max(n - e, -1.0), min(n + e, 1.0));
}
// The 2D overload: [n - L2 r, n + L2 r] intersected with [-1, 1], n = cultmath_snoise(centre). One snoise, value only.
float2 cultmath_iv_snoise_ball(float2 centre, float radius)
{
float n = cultmath_snoise(centre);
float e = CULTMATH_SNOISE2_LIPSCHITZ * radius;
return float2(max(n - e, -1.0), min(n + e, 1.0));
}
// Octave sum of cultmath_iv_snoise_ball, compounding frequency and amplitude as the value-and-gradient fBm does.
float2 cultmath_iv_fbm_ball(float3 centre, float radius, int octaves, float lacunarity, float gain)
{

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

@ -1,4 +1,4 @@
// Generated from shaders/CultMath.hlsl and CultMath.Interval.hlsl by GlslLowering (packages/cultmath/tests/CultMath.Tests/GlslLowering.cs). Do not edit; regenerate with CULTMATH_WRITE_GLSL=1 dotnet test --filter GlslMirrorTests. MIT. cultmath_phacelle is in CultMath.Phacelle.glsl, to be concatenated after this file by consumers that call it.
// Generated from shaders/CultMath.hlsl and CultMath.Interval.hlsl and CultMath.Affine.hlsl by GlslLowering (packages/cultmath/tests/CultMath.Tests/GlslLowering.cs). Do not edit; regenerate with CULTMATH_WRITE_GLSL=1 dotnet test --filter GlslMirrorTests. MIT. cultmath_phacelle is in CultMath.Phacelle.glsl, to be concatenated after this file by consumers that call it.
#ifndef CULTMATH_GLSL
#define CULTMATH_GLSL
// CultMath HLSL mirror.
@ -518,6 +518,9 @@ CultCellular cultmath_cellular(vec3 p)
// Lipschitz constant of cultmath_snoise(vec3); provenance in math.Interval.cs (SNOISE_LIPSCHITZ).
const float CULTMATH_SNOISE_LIPSCHITZ = 7.9640074;
// Lipschitz constant of cultmath_snoise(vec2); provenance in math.Interval.cs (SNOISE2_LIPSCHITZ).
const float CULTMATH_SNOISE2_LIPSCHITZ = 8.117002;
vec2 cultmath_iv_point(float value) { return vec2(value, value); }
vec2 cultmath_iv_add(vec2 a, vec2 b) { return vec2(a.x + b.x, a.y + b.y); }
vec2 cultmath_iv_sub(vec2 a, vec2 b) { return vec2(a.x - b.y, a.y - b.x); }
@ -596,6 +599,14 @@ vec2 cultmath_iv_snoise_ball(vec3 centre, float radius)
return vec2(max(n - e, -1.0), min(n + e, 1.0));
}
// The 2D overload: [n - L2 r, n + L2 r] intersected with [-1, 1], n = cultmath_snoise(centre). One snoise, value only.
vec2 cultmath_iv_snoise_ball(vec2 centre, float radius)
{
float n = cultmath_snoise(centre);
float e = CULTMATH_SNOISE2_LIPSCHITZ * radius;
return vec2(max(n - e, -1.0), min(n + e, 1.0));
}
// Octave sum of cultmath_iv_snoise_ball, compounding frequency and amplitude as the value-and-gradient fBm does.
vec2 cultmath_iv_fbm_ball(vec3 centre, float radius, int octaves, float lacunarity, float gain)
{
@ -631,6 +642,139 @@ vec4 cultmath_iv_frustum_ball(vec2 centreSlope, float z0, float z1, float footpr
// Reduced affine forms with one shared symbol over vec3(x0, a, e): the set { x0 + a eps + e delta :
// delta in [-1, 1] } at the region's shared eps in [-1, 1], e >= 0. Every cultmath_af_* function mirrors
// math.Affine.cs bit for bit; see that file and docs/design.md ("Affine forms") for each op's enclosure
// and the float32 rounding allowance it folds into e. Included by CultMath.hlsl after
// CultMath.Interval.hlsl, whose CULTMATH_SNOISE_LIPSCHITZ it shares, and after its own
// definition of cultmath_snoise_grad.
// An upper bound on the spectral norm of the Hessian of cultmath_snoise(vec3); provenance in
// math.Affine.cs (SNOISE_HESSIAN).
const float CULTMATH_SNOISE_HESSIAN = 56.050385;
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); }
vec3 cultmath_af_from_iv(vec2 interval)
{
float x0 = interval.x * 0.5 + interval.y * 0.5;
float e = interval.y * 0.5 - interval.x * 0.5;
return vec3(x0, 0.0, e + (abs(x0) + e) * 9.5367431640625e-7 + 1.17549435e-38);
}
// An endpoint that would be NaN is the infinity on its side. NaN is found by its bit pattern: dxc marks
// compares fast, so a test that needs NaN to make `lo == lo` false is deleted under its default flags.
vec2 cultmath_af_range(vec3 x)
{
float r = abs(x.y) + x.z;
float w = abs(x.x) * 9.5367431640625e-7 + r * 9.5367431640625e-7;
float lo = x.x - r - w;
float hi = x.x + r + w;
return vec2((floatBitsToUint(lo) & 0x7fffffffu) > 0x7f800000u ? uintBitsToFloat(0xff800000u) : lo,
(floatBitsToUint(hi) & 0x7fffffffu) > 0x7f800000u ? uintBitsToFloat(0x7f800000u) : hi);
}
vec3 cultmath_af_add(vec3 x, vec3 y)
{
float x0 = x.x + y.x;
float a = x.y + y.y;
float e = x.z + y.z;
return vec3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7);
}
vec3 cultmath_af_sub(vec3 x, vec3 y)
{
float x0 = x.x - y.x;
float a = x.y - y.y;
float e = x.z + y.z;
return vec3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7);
}
vec3 cultmath_af_neg(vec3 x) { return vec3(-x.x, -x.y, x.z); }
vec3 cultmath_af_scale(vec3 x, float s)
{
float x0 = x.x * s;
float a = x.y * s;
float e = x.z * abs(s);
return vec3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7 + 1.17549435e-38);
}
vec3 cultmath_af_add_iv(vec3 x, vec2 interval) { return cultmath_af_add(x, cultmath_af_from_iv(interval)); }
vec3 cultmath_af_mul(vec3 x, vec3 y)
{
float x0 = x.x * y.x;
float p = x.x * y.y;
float q = y.x * x.y;
float a = p + q;
float e = abs(x.x) * y.z + abs(y.x) * x.z + (abs(x.y) + x.z) * (abs(y.y) + y.z);
return vec3(x0, a, e + (abs(x0) + abs(p) + abs(q) + e) * 9.5367431640625e-7 + 1.17549435e-38);
}
// The narrower of the Lipschitz form and the centred (Taylor) form over every p = centre + axis eps + w,
// |w| <= radius, from one cultmath_snoise_grad at the centre; the allowance for the float32 evaluation is
// in math.Affine.cs.
vec3 cultmath_af_snoise(vec3 centre, vec3 axis, float radius)
{
vec4 g = cultmath_snoise_grad(centre);
vec3 gradient = vec3(g.x, g.y, g.z);
float reach = length(axis) + radius;
float lipschitz = CULTMATH_SNOISE_LIPSCHITZ * reach;
float a = dot(gradient, axis);
float remainder = length(gradient) * radius + 0.5 * CULTMATH_SNOISE_HESSIAN * reach * reach;
float allowance = 6.103515625e-5 + 2.44140625e-4 * reach
+ CULTMATH_SNOISE_LIPSCHITZ * 9.5367431640625e-7 * (abs(centre.x) + abs(centre.y) + abs(centre.z) + reach);
bool centred = abs(a) + remainder < lipschitz;
float slope = centred ? a : 0.0;
float e = (centred ? remainder : lipschitz) + allowance;
return vec3(g.w, slope, e + (abs(g.w) + abs(slope) + e) * 9.5367431640625e-7);
}
// Octave sum of cultmath_af_snoise, compounding frequency and amplitude as the value-and-gradient fBm does.
vec3 cultmath_af_fbm(vec3 centre, vec3 axis, float radius, int octaves, float lacunarity, float gain)
{
octaves = clamp(octaves, 0, 16);
float amplitude = 1.0;
float frequency = 1.0;
vec3 sum = vec3(0.0, 0.0, 0.0);
for (int i = 0; i < octaves; i++)
{
sum = cultmath_af_add(sum, cultmath_af_scale(cultmath_af_snoise(centre * frequency, axis * frequency, radius * abs(frequency)), amplitude));
frequency *= lacunarity;
amplitude *= gain;
}
return sum;
}
// The axis of a screen tile's slice over the depth segment [z0, z1]; camera frame of
// cultmath_iv_frustum_ball.
vec3 cultmath_af_frustum_axis(vec2 centreSlope, float z0, float z1)
{
float h = (z1 - z0) * 0.5;
return vec3(centreSlope.x * h, centreSlope.y * h, h);
}
// The affine slice of a screen tile over [z0, z1]: centre xyz, radius w. Derivation in math.Affine.cs.
vec4 cultmath_af_frustum_ball(vec2 centreSlope, float z0, float z1, float footprintPerDepth, float warpVariation)
{
float zm = (z0 + z1) * 0.5;
float h = (z1 - z0) * 0.5;
float cx = centreSlope.x * zm;
float cy = centreSlope.y * zm;
float spread = (abs(centreSlope.x) + abs(centreSlope.y) + 1.0) * h;
float radius = z1 * footprintPerDepth + warpVariation;
radius += (abs(cx) + abs(cy) + abs(zm) + spread + radius) * 9.5367431640625e-7;
return vec4(cx, cy, zm, radius);
}
float cultmath_value_noise(vec2 position)
{
vec2 cell = floor(position);

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

@ -509,6 +509,7 @@ CultCellular cultmath_cellular(float3 p)
#include "CultMath.Phacelle.hlsl"
#include "CultMath.Interval.hlsl"
#include "CultMath.Affine.hlsl"
float cultmath_value_noise(float2 position)
{

241
packages/cultmath/src/CultMath/math.Affine.cs Normal file
View file

@ -0,0 +1,241 @@
namespace CultMath;
public static partial class math
{
// Reduced affine forms with one shared symbol (design.md, "Affine forms"): float3(x0, a, e) is the
// set { x0 + a eps + e delta : delta in [-1, 1] } at the region's shared eps in [-1, 1], with
// e >= 0. Every form over one region shares that eps (for a screen-tile slice, the depth
// parameter), and delta is private to the form. This is Messine's AF1 as Gamito and Maddock use it
// for gradient noise (t0 + t1 e1 + t2 e2), in one register, as an interval is a float2.
//
// Every af_* op returns a form whose set, at each eps, contains f evaluated pointwise in float32 at
// every point of its operands' sets at that same eps. Rounding rule: each op folds its own float32
// rounding into e as (S) 2^-20, where S is the sum of the magnitudes the op's roundings are
// relative to (named in each comment). With u = 2^-24, every op below needs at most 8u S (its
// component roundings, the pointwise result's own rounding, and the roundings of computing e), so
// 2^-20 = 16u covers it twice over; the factor of two is the allowance for a device whose + and *
// are faithful (within one ulp) rather than correctly rounded, and fma contraction only removes
// roundings. Products also add 2^-126, the most any underflow to a subnormal can lose. The
// enclosure tests evaluate the form's set in double and f in float32 with no tolerance
// (AffineTests.EveryOpEnclosesItsPoints). No op takes a compound type or a reference parameter: each is a
// plain function of float values, so it mirrors into the HLSL common subset.
// An upper bound on the spectral norm of the Hessian of snoise(float3). Provenance:
// NoiseBoundTests.MeasureHessian (slow, explicit) takes the largest Hessian norm (central
// differences of snoise_grad, largest absolute eigenvalue) over 1e6 seeded points, refines the 1e4
// largest by ascent, and multiplies the refined maximum by 1.10;
// NoiseBoundTests.HessianConstantPinsSampledCurvature pins this value against a fast re-measure and
// the ascent from MeasureHessian's witness. Empirical with a margin, not a proof. Random balls
// cannot see a constant 20% low, so the defence is those two curvature tests plus
// NoiseBoundTests.SnoiseFormEnclosesTheCurvatureWitness, where af_snoise misses at the point the
// curvature is realised if the constant is under it. The margin over the function: a scan of the
// whole 289^3 hash period found a realised maximum of 50.58, and bounding each corner's gradient
// independently over every hash configuration gives 52.22 (Soul pass s1). A change to the snoise
// kernel must re-run MeasureHessian and re-pin it.
public const float SNOISE_HESSIAN = 56.050385f;
/// <summary>The point form (x, 0, 0): the value x at every eps. Exact.</summary>
public static float3 af_point(float value) => new(value, 0.0f, 0.0f);
/// <summary>
/// The shared symbol itself, (x0, a, 0): the value x0 + a eps. It is how a region's parameter enters
/// (for a tile slice, depth z = z_m + h eps), so it is exact by definition: eps is defined through it.
/// </summary>
public static float3 af_symbol(float x0, float a) => new(x0, a, 0.0f);
/// <summary>
/// Every value of the interval [lo, hi], uncorrelated with eps: (lo / 2 + hi / 2, 0, hi / 2 - lo / 2),
/// the halves taken before the sum and the difference so that endpoints near the float maximum do not
/// overflow in them. Rounding: x0 and e each round once, and a value of the interval is within e of x0
/// only up to those; S = |x0| + e. Halving a subnormal endpoint loses at most 2^-150, and the sum and
/// difference of subnormals are exact, so 2^-126 covers the underflow (the e of an interval spanning
/// one subnormal step would otherwise round to 0). |x0| + e is at most the larger endpoint
/// in magnitude, so it does not overflow. If e + its absorption exceeds the float maximum it is
/// infinity, which still encloses.
/// </summary>
public static float3 af_from_iv(float2 interval)
{
var x0 = interval.x * 0.5f + interval.y * 0.5f;
var e = interval.y * 0.5f - interval.x * 0.5f;
return new float3(x0, 0.0f, e + (abs(x0) + e) * 9.5367431640625e-7f + 1.17549435e-38f);
}
/// <summary>
/// The interval of every value of the form over every eps: x0 -+ (|a| + e), the bridge to iv_*. The
/// endpoints round inward by at most u (|x0| + |a| + e) each, so each is moved outward by
/// (|x0| + |a| + e) 2^-20, formed term by term so that the sum never passes the float maximum before it
/// is scaled. A form with an infinite component (an overflowed op) gives infinite
/// endpoints, and an endpoint that would be NaN (infinity minus infinity, or a NaN component) is
/// the infinity on its side, so the result is never NaN and still encloses. NaN is found by its bit
/// pattern, not by <c>lo == lo</c>, which a compiler that assumes no NaNs (dxc's default) deletes.
/// </summary>
public static float2 af_range(float3 x)
{
var r = abs(x.y) + x.z;
var w = abs(x.x) * 9.5367431640625e-7f + r * 9.5367431640625e-7f;
var lo = x.x - r - w;
var hi = x.x + r + w;
return new float2((asuint(lo) & 0x7fffffffu) > 0x7f800000u ? float.NegativeInfinity : lo,
(asuint(hi) & 0x7fffffffu) > 0x7f800000u ? float.PositiveInfinity : hi);
}
/// <summary>
/// Encloses x + y: (x0 + y0, a_x + a_y, e_x + e_y). The shared symbol keeps the two forms
/// correlated, so af_sub(x, x) has a = 0. Rounding: S = |x0| + |a| + e of the result.
/// </summary>
public static float3 af_add(float3 x, float3 y)
{
var x0 = x.x + y.x;
var a = x.y + y.y;
var e = x.z + y.z;
return new float3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7f);
}
/// <summary>Encloses x - y: (x0 - y0, a_x - a_y, e_x + e_y). Rounding as af_add.</summary>
public static float3 af_sub(float3 x, float3 y)
{
var x0 = x.x - y.x;
var a = x.y - y.y;
var e = x.z + y.z;
return new float3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7f);
}
/// <summary>Encloses -x: (-x0, -a, e). Exact.</summary>
public static float3 af_neg(float3 x) => new(-x.x, -x.y, x.z);
/// <summary>
/// Encloses x * s: (x0 s, a s, e |s|). Rounding: S = |x0| + |a| + e of the result, plus 2^-126
/// for underflow.
/// </summary>
public static float3 af_scale(float3 x, float s)
{
var x0 = x.x * s;
var a = x.y * s;
var e = x.z * abs(s);
return new float3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7f + 1.17549435e-38f);
}
/// <summary>Encloses x + v for v in the interval, uncorrelated with eps: af_add(x, af_from_iv(interval)).</summary>
public static float3 af_add_iv(float3 x, float2 interval) => af_add(x, af_from_iv(interval));
/// <summary>
/// Encloses x * y, AF1's product: (x0 y0, x0 a_y + y0 a_x, |x0| e_y + |y0| e_x + (|a_x| + e_x)(|a_y| + e_y)).
/// The last term condenses the eps^2 term a_x a_y and every cross term with an error symbol into e,
/// so a form multiplied by itself is wider than its square: there is no af_sqr until a consumer
/// names one. Rounding: a can cancel, so S = |x0| + |x0 a_y| + |y0 a_x| + e, plus 2^-126 for
/// underflow.
/// </summary>
public static float3 af_mul(float3 x, float3 y)
{
var x0 = x.x * y.x;
var p = x.x * y.y;
var q = y.x * x.y;
var a = p + q;
var e = abs(x.x) * y.z + abs(y.x) * x.z + (abs(x.y) + x.z) * (abs(y.y) + y.z);
return new float3(x0, a, e + (abs(x0) + abs(p) + abs(q) + e) * 9.5367431640625e-7f + 1.17549435e-38f);
}
/// <summary>
/// Encloses snoise(p) for every p = centre + axis eps + w with |w| &lt;= radius, from one snoise_grad at
/// the centre (value n, gradient g). With R = |axis| + radius it returns the narrower of two sound
/// forms: the Lipschitz form (n, 0, L R), L = SNOISE_LIPSCHITZ, and the centred form
/// (n, g . axis, |g| radius + M R^2 / 2), M = SNOISE_HESSIAN. The centred form is Taylor's theorem:
/// snoise(c + d) = n + g . d + r with |r| &lt;= M |d|^2 / 2, d = axis eps + w, g . w &lt;= |g| radius.
/// radius must be at least 0.
///
/// Both forms carry an allowance for the float32 evaluation of snoise, which is not the exact
/// function: 2^-14 + 2^-12 R + L 2^-20 (|centre|_1 + R). The skew and the corner offsets round the
/// input point by at most about 2u |p|_1 before any kernel term is formed (the hashing is exact on
/// integers below 2^24), so a float32 snoise of p is the exact one at a point within 2u |p|_1, which
/// moves the value by at most L 2u |p|_1; for both p and the centre, with the caller's scaling of the
/// point by a frequency (one rounding) and a factor of two for a GPU, 16u = 2^-20. The kernel's
/// absolute rounding is at most 105 x four corners x terms of magnitude below 0.5^4 x 1.8 over about
/// 30 operations, under 2^-15, doubled to 2^-14; the gradient's error times |d| &lt;= R is under 2^-13 R,
/// doubled to 2^-12 R. AffineTests.SnoiseFormEnclosesBall checks the result with no tolerance on the
/// CPU, and the fixture's enclosure check measures it on a GPU. Rounding of the form itself: S =
/// |n| + |a| + e. The Lipschitz form is chosen when its half-width L R is at most the centred form's
/// |a| + |g| radius + M R^2 / 2, so the total width is never more than the Lipschitz form's.
/// </summary>
public static float3 af_snoise(float3 centre, float3 axis, float radius)
{
var g = snoise_grad(centre);
var gradient = new float3(g.x, g.y, g.z);
var reach = length(axis) + radius;
var lipschitz = SNOISE_LIPSCHITZ * reach;
var a = dot(gradient, axis);
var remainder = length(gradient) * radius + 0.5f * SNOISE_HESSIAN * reach * reach;
var allowance = 6.103515625e-5f + 2.44140625e-4f * reach
+ SNOISE_LIPSCHITZ * 9.5367431640625e-7f * (abs(centre.x) + abs(centre.y) + abs(centre.z) + reach);
var centred = abs(a) + remainder < lipschitz;
var slope = centred ? a : 0.0f;
var e = (centred ? remainder : lipschitz) + allowance;
return new float3(g.w, slope, e + (abs(g.w) + abs(slope) + e) * 9.5367431640625e-7f);
}
/// <summary>
/// Encloses the fBm value, sum of a_i snoise(p f_i) over the octaves, for every p = centre + axis eps
/// + w with |w| &lt;= radius: the octave sum of af_snoise(centre f_i, axis f_i, radius |f_i|) scaled by
/// a_i, with frequency f_i and amplitude a_i compounded exactly as fbm_grad compounds them (its .w is
/// what AffineTests.SnoiseFormEnclosesBall encloses). The rounding of p f_i, inside fbm_grad, and of
/// the scaled centre and axis here, moves each point by at most 3u (|centre f_i|_1 + R f_i), inside
/// af_snoise's allowance for its input point. octaves is clamped to [0, 16] as there.
/// </summary>
public static float3 af_fbm(float3 centre, float3 axis, float radius, int octaves, float lacunarity, float gain)
{
octaves = clamp(octaves, 0, 16);
var amplitude = 1.0f;
var frequency = 1.0f;
var sum = new float3(0.0f, 0.0f, 0.0f);
for (var i = 0; i < octaves; i++)
{
sum = af_add(sum, af_scale(af_snoise(centre * frequency, axis * frequency, radius * abs(frequency)), amplitude));
frequency *= lacunarity;
amplitude *= gain;
}
return sum;
}
/// <summary>
/// The axis of a screen tile's slice over the depth segment [z0, z1], for af_frustum_ball's centre:
/// (m_c, 1) h with h = (z1 - z0) / 2, in the camera frame of iv_frustum_ball (a pinhole at the
/// origin looking down +z, the ray of slope m the points (m z, z)). Its z component is h itself, the
/// symbol's coefficient for depth: z = z_m + h eps. Its rounding is af_frustum_ball's.
/// </summary>
public static float3 af_frustum_axis(float2 centreSlope, float z0, float z1)
{
var h = (z1 - z0) * 0.5f;
return new float3(centreSlope.x * h, centreSlope.y * h, h);
}
/// <summary>
/// The affine slice of a screen tile over [z0, z1]: a centre (xyz) and radius (w) such that every point
/// of every ray of the tile at depth z in [z0, z1] is centre + axis eps(z) + w', |w'| &lt;= radius, with
/// axis = af_frustum_axis(centreSlope, z0, z1), z_m = (z0 + z1) / 2, h = (z1 - z0) / 2 (both float32,
/// as computed here) and eps(z) = clamp((z - z_m) / h, -1, 1) (0 where h is 0). The centre is
/// (m_c z_m, z_m) and the radius z1 footprintPerDepth + warpVariation, where footprintPerDepth is
/// iv_frustum_ball's and warpVariation is what the caller's warp can add once its value at the centre
/// is applied to the centre (design.md, "Affine forms", the centred warp). Derivation: a point of
/// slope m_c + e (|e| &lt;= footprintPerDepth) at depth z is (m_c z_m, z_m) + (m_c, 1)(z - z_m) +
/// (e z, 0), and (m_c, 1)(z - z_m) = axis eps(z) + (m_c, 1) r(z), where the clamp residual
/// |r(z)| = |z - z_m - h eps(z)| is at most u (|z_m| + h), one rounding each of z_m and h. A caller
/// taking depth as af_symbol(z_m, h) adds the same residual, 2^-23 (|z_m| + h) covers it, to that
/// form's e. Rounding: the centre's products and the axis's products round once each (u |c|_1 and
/// u |axis|_1 against |eps| &lt;= 1), the residual is at most u (|c|_1 + |axis|_1) since
/// |(m_c, 1)| &lt;= |m_c|_1 + 1, and the radius rounds at most four times; 2^-20 (|c|_1 + |axis|_1 +
/// radius) covers it twice over (AffineTests.FrustumSliceEnclosesEveryRayPoint, extreme families
/// included). The caller rotates and translates the centre into world space, rotates the axis, and
/// scales all three by its noise frequency; the rounding of that transform is the caller's.
/// </summary>
public static float4 af_frustum_ball(float2 centreSlope, float z0, float z1, float footprintPerDepth, float warpVariation)
{
var zm = (z0 + z1) * 0.5f;
var h = (z1 - z0) * 0.5f;
var cx = centreSlope.x * zm;
var cy = centreSlope.y * zm;
var spread = (abs(centreSlope.x) + abs(centreSlope.y) + 1.0f) * h;
var radius = z1 * footprintPerDepth + warpVariation;
radius += (abs(cx) + abs(cy) + abs(zm) + spread + radius) * 9.5367431640625e-7f;
return new float4(cx, cy, zm, radius);
}
}

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

@ -0,0 +1,2 @@
fileFormatVersion: 2
guid: 86b87ffb6c8a4d72b067a64f9f8a06ce

24
packages/cultmath/src/CultMath/math.Interval.cs
View file

@ -20,6 +20,14 @@ public static partial class math
// snoise kernel must re-run MeasureLipschitz and re-pin this constant.
public const float SNOISE_LIPSCHITZ = 7.9640074f;
// Lipschitz constant of snoise(float2) in its input, with SNOISE_LIPSCHITZ's provenance rule:
// NoiseBoundTests.MeasureLipschitz2 (slow, explicit) takes the largest gradient length (central
// differences over exact float steps) over 1e6 seeded points, refines the 1e4 largest by ascent and
// multiplies by 1.10; NoiseBoundTests.Lipschitz2ConstantPins pins it against a fast re-measure and
// the ascent from MeasureLipschitz2's witness. A change to the 2D snoise kernel must re-run
// MeasureLipschitz2 and re-pin this constant.
public const float SNOISE2_LIPSCHITZ = 8.117002f;
/// <summary>The point interval [x, x].</summary>
public static float2 iv_point(float value) => new(value, value);
@ -142,6 +150,22 @@ public static partial class math
return new(max(n - e, -1.0f), min(n + e, 1.0f));
}
/// <summary>
/// Encloses snoise(x) for every 2D x within distance radius of centre: [n - L2 r, n + L2 r] intersected
/// with [-1, 1], n = snoise(centre), L2 = SNOISE2_LIPSCHITZ. radius must be at least 0. One snoise
/// evaluation, value only. Like iv_snoise_ball(float3) it is exact for the point it is handed and
/// carries no allowance for the float32 evaluation: a consumer whose 2D flow is evaluated at points
/// it rounds itself (a shifted, scaled or projected point) adds 2^-14 + L2 2^-20 (|centre|_1 + radius)
/// to the width, the convention af_snoise carries for the 3D point
/// (NoiseBoundTests.FlowVariationCarriesItsAllowance, CentredWarpStaysEnclosed).
/// </summary>
public static float2 iv_snoise_ball(float2 centre, float radius)
{
var n = snoise(centre);
var e = SNOISE2_LIPSCHITZ * radius;
return new(max(n - e, -1.0f), min(n + e, 1.0f));
}
/// <summary>
/// Encloses the fBm value, sum of a_i * snoise(x * f_i) over the octaves, for every x within distance
/// radius of centre: the octave sum of iv_snoise_ball(centre * f_i, radius * |f_i|) scaled by a_i,

699
packages/cultmath/tests/CultMath.Tests/AffineTests.cs Normal file
View file

@ -0,0 +1,699 @@
using CultMath;
using Xunit;
using static CultMath.math;
namespace CultMath.Tests;
/// <summary>
/// The affine forms of math.Affine.cs (design.md, "Affine forms"). A form (x0, a, e) is the set
/// x0 + a eps + e delta at the region's shared eps; every check evaluates that set in double and the
/// pointwise function in float32, with no tolerance. Operand values are float32 values drawn inside
/// their operands' sets at one shared eps, endpoints (eps and delta at -1 and 1) included.
/// </summary>
public sealed class AffineTests
{
private readonly ITestOutputHelper output;
public AffineTests(ITestOutputHelper output) => this.output = output;
private static float Uniform(System.Random random, float lo, float hi) => lo + (hi - lo) * random.NextSingle();
private static float LogUniform(System.Random random, float lo, float hi) =>
MathF.Exp(Uniform(random, MathF.Log(lo), MathF.Log(hi)));
private static float Signed(System.Random random, float magnitude) => random.Next(2) == 0 ? -magnitude : magnitude;
// Whether value lies in the form's set at eps, in double.
private static bool InForm(float3 form, float eps, double value) =>
Math.Abs(value - (form.x + (double)form.y * eps)) <= form.z;
// A form with its centre in [-20, 20] and |a| and e log-uniform over [1e-5, 20]; one in eight has
// a = 0 and one in eight e = 0.
private static float3 DrawForm(System.Random random)
{
var x0 = Uniform(random, -20.0f, 20.0f);
var a = random.Next(8) == 0 ? 0.0f : Signed(random, LogUniform(random, 1.0e-5f, 20.0f));
var e = random.Next(8) == 0 ? 0.0f : LogUniform(random, 1.0e-5f, 20.0f);
return new float3(x0, a, e);
}
private static float2 DrawInterval(System.Random random)
{
var lo = Uniform(random, -20.0f, 20.0f);
return new float2(lo, lo + (random.Next(8) == 0 ? 0.0f : LogUniform(random, 1.0e-5f, 40.0f)));
}
// The tiny family: magnitudes from 1e-21 to 1e-17, products near and under 2^-126 = 1.2e-38, where a
// float32 result rounds to a subnormal and the absorption S 2^-20 is smaller than that loss.
private static float3 DrawTiny(System.Random random) =>
new(random.Next(8) == 0 ? 0.0f : Signed(random, LogUniform(random, 1.0e-21f, 1.0e-17f)),
random.Next(4) == 0 ? 0.0f : Signed(random, LogUniform(random, 1.0e-21f, 1.0e-17f)),
random.Next(4) == 0 ? 0.0f : LogUniform(random, 1.0e-21f, 1.0e-17f));
private static float2 DrawTinyInterval(System.Random random)
{
var lo = Signed(random, LogUniform(random, 1.0e-21f, 1.0e-17f));
return new float2(lo, lo + (random.Next(8) == 0 ? 0.0f : LogUniform(random, 1.0e-21f, 1.0e-17f)));
}
// A float32 value in the form's set at eps, delta at an endpoint when asked; NaN when rounding leaves
// no float of the set near the drawn point (a set narrower than an ulp).
private static float Member(System.Random random, float3 form, float eps, bool endpoint)
{
var centre = form.x + (double)form.y * eps;
var delta = endpoint ? (random.Next(2) == 0 ? -1.0 : 1.0) : Uniform(random, -1.0f, 1.0f);
var value = (float)(centre + form.z * delta);
for (var k = 0; k < 4 && !InForm(form, eps, value); k++)
value = value > centre ? MathF.BitDecrement(value) : MathF.BitIncrement(value);
return InForm(form, eps, value) ? value : float.NaN;
}
private static float MemberOf(System.Random random, float2 interval, bool endpoint) =>
endpoint ? (random.Next(2) == 0 ? interval.x : interval.y) : Math.Clamp(Uniform(random, interval.x, interval.y), interval.x, interval.y);
private readonly record struct Operands(float3 X, float3 Y, float2 Iv, float S);
private readonly record struct Values(float X, float Y, float V, float S);
// One op under the harness: Build gives the result's membership test (eps, value) from the operand
// forms, F the pointwise float32 function of operand values.
private sealed record Op(string Name, Func<Operands, Func<float, double, bool>> Build, Func<Values, float> F);
private static Func<float, double, bool> Form(float3 form) => (eps, value) => InForm(form, eps, value);
private static Func<float, double, bool> Range(float2 range) => (_, value) => value >= range.x && value <= range.y;
private static readonly Op[] Ops =
{
new("af_point", o => Form(af_point(o.X.x)), v => v.S),
new("af_symbol", o => Form(af_symbol(o.X.x, o.X.y)), v => v.S),
new("af_from_iv", o => Form(af_from_iv(o.Iv)), v => v.V),
new("af_range", o => Range(af_range(o.X)), v => v.X),
new("af_add", o => Form(af_add(o.X, o.Y)), v => v.X + v.Y),
new("af_sub", o => Form(af_sub(o.X, o.Y)), v => v.X - v.Y),
new("af_neg", o => Form(af_neg(o.X)), v => -v.X),
new("af_scale", o => Form(af_scale(o.X, o.S)), v => v.X * v.S),
new("af_add_iv", o => Form(af_add_iv(o.X, o.Iv)), v => v.X + v.V),
new("af_mul", o => Form(af_mul(o.X, o.Y)), v => v.X * v.Y),
};
// Runs one op over `forms` seeded operand draws x 16 (eps, delta) draws; returns the misses and up to
// four witnesses. The first four draws per form put eps at an endpoint, the first eight delta.
// af_point and af_symbol are checked on their defining value: the point x0 itself, and x0 + a eps in
// double (the symbol defines eps, so it has no float32 evaluation to enclose).
private static (int Misses, List<string> Witness) Harness(Op op, int seed, int forms = 10_000, Func<System.Random, float3>? draw = null, bool tiny = false)
{
draw ??= tiny ? DrawTiny : DrawForm;
var random = new System.Random(seed);
var misses = 0;
var witness = new List<string>();
for (var t = 0; t < forms; t++)
{
var operands = new Operands(draw(random), draw(random), tiny ? DrawTinyInterval(random) : DrawInterval(random),
random.Next(8) == 0 ? 0.0f : tiny ? Signed(random, LogUniform(random, 1.0e-32f, 1.0e-19f)) : Uniform(random, -4.0f, 4.0f));
var contains = op.Build(operands);
for (var d = 0; d < 16; d++)
{
var eps = d < 4 ? (d % 2 == 0 ? -1.0f : 1.0f) : Uniform(random, -1.0f, 1.0f);
var endpoint = d < 8;
var x = Member(random, operands.X, eps, endpoint);
var y = Member(random, operands.Y, eps, endpoint);
if (float.IsNaN(x) || float.IsNaN(y))
continue;
var v = MemberOf(random, operands.Iv, endpoint);
double f = op.Name switch
{
"af_point" => operands.X.x,
"af_symbol" => operands.X.x + (double)operands.X.y * eps,
_ => op.F(new Values(x, y, v, operands.S)),
};
if (!contains(eps, f))
{
misses++;
if (witness.Count < 4)
witness.Add($"{op.Name}: X {operands.X}, Y {operands.Y}, iv {operands.Iv}, s {operands.S:R}; eps {eps:R}, values ({x:R}, {y:R}, {v:R}) give {f:R}");
}
}
}
return (misses, witness);
}
[Fact]
public void EveryOpEnclosesItsPoints()
{
var failures = new List<string>();
foreach (var op in Ops)
{
var (misses, witness) = Harness(op, 0xAF00);
if (misses > 0)
failures.Add($"{op.Name}: {misses} misses" + Environment.NewLine + string.Join(Environment.NewLine, witness));
}
Assert.True(failures.Count == 0, string.Join(Environment.NewLine, failures));
}
private static float3 Shrink(float3 form) => new(form.x, form.y, form.z > 0.0f ? MathF.BitDecrement(form.z) : 0.0f);
// Each op with its rounding absorption removed: the propagated error alone, computed in float32.
private static readonly Op[] Unabsorbed =
{
new("af_from_iv", o => Form(new float3((o.Iv.x + o.Iv.y) * 0.5f, 0.0f, (o.Iv.y - o.Iv.x) * 0.5f)), v => v.V),
new("af_range", o => Range(new float2(o.X.x - (MathF.Abs(o.X.y) + o.X.z), o.X.x + (MathF.Abs(o.X.y) + o.X.z))), v => v.X),
new("af_add", o => Form(new float3(o.X.x + o.Y.x, o.X.y + o.Y.y, o.X.z + o.Y.z)), v => v.X + v.Y),
new("af_sub", o => Form(new float3(o.X.x - o.Y.x, o.X.y - o.Y.y, o.X.z + o.Y.z)), v => v.X - v.Y),
new("af_scale", o => Form(new float3(o.X.x * o.S, o.X.y * o.S, o.X.z * MathF.Abs(o.S))), v => v.X * v.S),
new("af_add_iv", o => Form(new float3(o.X.x + (o.Iv.x + o.Iv.y) * 0.5f, o.X.y, o.X.z + (o.Iv.y - o.Iv.x) * 0.5f)), v => v.X + v.V),
new("af_mul", o => Form(new float3(o.X.x * o.Y.x, o.X.x * o.Y.y + o.Y.x * o.X.y,
MathF.Abs(o.X.x) * o.Y.z + MathF.Abs(o.Y.x) * o.X.z + (MathF.Abs(o.X.y) + o.X.z) * (MathF.Abs(o.Y.y) + o.Y.z))), v => v.X * v.Y),
};
// One-ulp mutants that survive, each with the reason. A rounding absorption bounds the worst case
// over every operand value, and 2^-20 S is 16u S, at least eight ulps of the larger of |x0|, |a| and
// e; one ulp less of e leaves the rest of that margin, which no sampled rounding can cross.
private static readonly Dictionary<string, string> OneUlpSurvivors = new()
{
["af_from_iv"] = "the absorption (|x0| + e) 2^-20 is at least eight ulps of e over the one rounding of x0 and e it covers",
["af_range"] = "each endpoint moves outward by (|x0| + |a| + e) 2^-20, eight ulps past the one rounding it covers",
["af_add"] = "the absorption 2^-20 S is eight ulps of S's largest term past the three roundings it covers",
["af_sub"] = "as af_add",
["af_scale"] = "as af_add",
["af_add_iv"] = "two absorptions, af_from_iv's and af_add's",
["af_mul"] = "the absorption 2^-20 (|x0| + |x0 a_y| + |y0 a_x| + e) is eight ulps past the products it covers",
};
/// <summary>
/// The rounding absorptions are load-bearing: with each op's absorption removed (Unabsorbed), the
/// harness finds a miss for every op. Each op's e shrunk by one ulp is run too; it survives for every
/// op, because the absorption is a worst-case bound many ulps wide, and each survivor is named with
/// its reason in OneUlpSurvivors. A one-ulp mutant that starts failing, or one that is not named,
/// fails this test, so the list cannot go stale.
/// </summary>
[Fact]
public void EveryRoundingTermIsCaught()
{
var problems = new List<string>();
foreach (var mutant in Unabsorbed)
{
var (misses, _) = Harness(mutant, 0xAF01);
output.WriteLine($"AF-REPORT {mutant.Name} without absorption: {misses} misses");
if (misses == 0)
problems.Add($"{mutant.Name}: removing the rounding absorption is not caught");
}
foreach (var op in Ops.Where(o => o.Name is not ("af_point" or "af_symbol" or "af_neg")))
{
var shrunk = op.Name == "af_range"
? op with { Build = o => { var r = af_range(o.X); return Range(new float2(MathF.BitIncrement(r.x), MathF.BitDecrement(r.y))); } }
: op with { Build = o => Form(Shrink(BuildForm(op.Name, o))) };
var (misses, _) = Harness(shrunk, 0xAF02);
output.WriteLine($"AF-REPORT {op.Name} with e one ulp short: {misses} misses");
var named = OneUlpSurvivors.ContainsKey(op.Name);
if (misses == 0 && !named)
problems.Add($"{op.Name}: e one ulp short survives and is not named");
if (misses > 0 && named)
problems.Add($"{op.Name}: e one ulp short is now caught; remove it from OneUlpSurvivors");
}
Assert.True(problems.Count == 0, string.Join(Environment.NewLine, problems));
}
private static float3 BuildForm(string name, Operands o) => name switch
{
"af_from_iv" => af_from_iv(o.Iv),
"af_add" => af_add(o.X, o.Y),
"af_sub" => af_sub(o.X, o.Y),
"af_scale" => af_scale(o.X, o.S),
"af_add_iv" => af_add_iv(o.X, o.Iv),
"af_mul" => af_mul(o.X, o.Y),
_ => throw new ArgumentException(name),
};
/// <summary>
/// The shared symbol correlates: for x = af_symbol(x0, a), af_sub(x, x) and af_add(x, af_neg(x)) have
/// a = 0 and e at most their own rounding, and af_sub(af_scale(x, 2), af_add(x, x)) has a = 0. A
/// mutant of af_add, af_sub or af_scale that condensed a into e would leave e of order |a|.
/// </summary>
[Fact]
public void SharedSymbolCorrelates()
{
var random = new System.Random(0xAF03);
for (var t = 0; t < 10_000; t++)
{
var x = af_symbol(Uniform(random, -20.0f, 20.0f), Signed(random, LogUniform(random, 1.0e-5f, 20.0f)));
foreach (var (name, form) in new[]
{
("af_sub(x, x)", af_sub(x, x)),
("af_add(x, af_neg(x))", af_add(x, af_neg(x))),
("af_sub(af_scale(x, 2), af_add(x, x))", af_sub(af_scale(x, 2.0f), af_add(x, x))),
})
{
var rounding = (MathF.Abs(x.x) + MathF.Abs(x.y)) * 1.0e-5f;
Assert.True(form.y == 0.0f && form.z <= rounding, $"{name} for x = {x} is {form}: a must be 0 and e at most {rounding:R}");
}
}
}
// Forms whose range straddles zero (|x0| below |a| + e) or is one-signed (|x0| above it).
private static float3 DrawSigned(System.Random random, bool straddle)
{
var a = Signed(random, LogUniform(random, 1.0e-3f, 10.0f));
var e = LogUniform(random, 1.0e-3f, 10.0f);
var reach = MathF.Abs(a) + e;
return new float3(Signed(random, straddle ? Uniform(random, 0.0f, 0.9f * reach) : reach * Uniform(random, 1.1f, 4.0f)), a, e);
}
/// <summary>
/// af_mul encloses the product over forms that straddle zero and forms that are one-signed, every
/// sign combination. The eps^2 term a_x a_y is what the condensed e carries for eps at +-1: a mutant
/// dropping |a_x a_y| from e misses there.
/// </summary>
[Fact]
public void MulEnclosesEverySignCase()
{
var mul = Ops.Single(o => o.Name == "af_mul");
var dropped = mul with
{
Build = o =>
{
var r = af_mul(o.X, o.Y);
return Form(new float3(r.x, r.y, r.z - MathF.Abs(o.X.y * o.Y.y)));
},
};
foreach (var (sx, sy) in new[] { (true, true), (true, false), (false, true), (false, false) })
{
var seed = 0xAF04 + (sx ? 1 : 0) + (sy ? 2 : 0);
var flip = 0;
float3 Draw(System.Random r) => DrawSigned(r, (flip++ % 2 == 0) ? sx : sy);
var (misses, witness) = Harness(mul, seed, 5_000, Draw);
Assert.True(misses == 0, $"af_mul, straddling ({sx}, {sy}): {misses} misses" + Environment.NewLine + string.Join(Environment.NewLine, witness));
flip = 0;
var (caught, _) = Harness(dropped, seed, 5_000, Draw);
Assert.True(caught > 0, $"af_mul without |a_x a_y|, straddling ({sx}, {sy}), is not caught");
}
}
// af_snoise's float32 allowance, as its doc comment states it, for the width comparison below.
private static double Allowance(float3 centre, double reach) =>
Math.ScaleB(1.0, -14) + Math.ScaleB(reach, -12)
+ SNOISE_LIPSCHITZ * Math.ScaleB(Math.Abs(centre.x) + Math.Abs(centre.y) + Math.Abs(centre.z) + reach, -20);
private static float3 UnitVector(System.Random random)
{
float3 v;
do
{
v = new float3(Uniform(random, -1.0f, 1.0f), Uniform(random, -1.0f, 1.0f), Uniform(random, -1.0f, 1.0f));
}
while (dot(v, v) > 1.0f || dot(v, v) < 1.0e-4f);
return normalize(v);
}
// A float32 point of the set centre + axis eps + w, |w| <= radius, checked in double; the first eight
// w are on the surface. Where the float32 sum rounds a point out of the set, w is pulled in by the
// overshoot and 2^-21 (|centre|_1 + |axis|_1) and the sum retried; a point that still rounds out is
// reported as null and skipped.
private static float3? PointOf(System.Random random, float3 centre, float3 axis, float radius, float eps, int index)
{
var w = UnitVector(random) * (index < 8 ? radius : radius * MathF.Cbrt(random.NextSingle()));
for (var k = 0; k < 6; k++)
{
var p = centre + axis * eps + w;
var (dx, dy, dz) = (p.x - (centre.x + (double)axis.x * eps), p.y - (centre.y + (double)axis.y * eps), p.z - (centre.z + (double)axis.z * eps));
var distance = Math.Sqrt(dx * dx + dy * dy + dz * dz);
if (distance <= radius)
return p;
w *= (float)(Math.Max(0.0, 2.0 * radius - distance - Math.ScaleB(Math.Abs(centre.x) + Math.Abs(centre.y) + Math.Abs(centre.z) + Math.Abs(axis.x) + Math.Abs(axis.y) + Math.Abs(axis.z), -21)) / radius);
}
return null;
}
/// <summary>
/// af_snoise and af_fbm over 2,000 (centre, axis, radius) each: centres in [-20, 20]^3, radius
/// log-uniform in [1e-4, 2], |axis| uniform in [0, 2] (an eighth exactly 0). 64 points
/// centre + axis eps + w per case, eps at -1 and 1 for the first eight; snoise (fbm_grad(...).w) in
/// float32 lies in the form at eps, no tolerance. Dropping the M R^2 / 2 remainder from a centred
/// form misses. af_snoise's total width is never more than the Lipschitz form's, L R plus the
/// allowance.
/// </summary>
[Fact]
public void SnoiseFormEnclosesBall()
{
var random = new System.Random(0xAF10);
var (misses, remainderMisses, skipped, centredCases) = (0, 0, 0, 0);
var witness = new List<string>();
for (var b = 0; b < 4000; b++)
{
var fbm = b >= 2000;
var centre = new float3(Uniform(random, -20.0f, 20.0f), Uniform(random, -20.0f, 20.0f), Uniform(random, -20.0f, 20.0f));
var axis = random.Next(8) == 0 ? new float3(0.0f, 0.0f, 0.0f) : UnitVector(random) * Uniform(random, 0.0f, 2.0f);
var radius = LogUniform(random, 1.0e-4f, 2.0f);
var (octaves, lacunarity, gain) = (random.Next(0, 7), Uniform(random, 1.5f, 2.5f), Uniform(random, 0.3f, 0.7f));
var form = fbm ? af_fbm(centre, axis, radius, octaves, lacunarity, gain) : af_snoise(centre, axis, radius);
var reach = (double)length(axis) + radius;
var remainder = 0.5 * SNOISE_HESSIAN * reach * reach;
var mutant = new float3(form.x, form.y, (float)(form.z - remainder));
if (!fbm)
{
var lipschitz = (SNOISE_LIPSCHITZ * reach + Allowance(centre, reach) + Math.ScaleB(Math.Abs(form.x), -20)) * (1.0 + Math.ScaleB(1.0, -18));
Assert.True(Math.Abs(form.y) + form.z <= lipschitz, $"af_snoise({centre}, {axis}, {radius:R}) = {form} is wider than the Lipschitz form's half-width {lipschitz:R}");
centredCases += form.y != 0.0f ? 1 : 0;
}
for (var i = 0; i < 64; i++)
{
var eps = i < 8 ? (i % 2 == 0 ? -1.0f : 1.0f) : Uniform(random, -1.0f, 1.0f);
var p = PointOf(random, centre, axis, radius, eps, i);
if (p is not { } point)
{
skipped++;
continue;
}
var value = fbm ? fbm_grad(point, octaves, lacunarity, gain).w : snoise(point);
if (!InForm(form, eps, value))
{
misses++;
if (witness.Count < 4)
witness.Add($"{(fbm ? $"af_fbm(..., {octaves}, {lacunarity:R}, {gain:R})" : "af_snoise")}({centre}, {axis}, {radius:R}) = {form} misses {value:R} at eps {eps:R}, {point}");
}
if (!fbm && form.y != 0.0f && !InForm(mutant, eps, value))
remainderMisses++;
}
}
output.WriteLine($"AF-REPORT af_snoise/af_fbm enclosure: {misses} misses, {skipped} points skipped (rounded out of the set), {centredCases} of 2000 af_snoise cases centred, {remainderMisses} misses with the remainder dropped");
Assert.True(misses == 0, string.Join(Environment.NewLine, witness));
Assert.True(skipped < 64 * 4000 / 100, $"{skipped} points rounded out of their sets");
Assert.True(remainderMisses > 0, "dropping M R^2 / 2 from the centred form is not caught");
}
/// <summary>
/// At radius 1e-3 and |axis| = 1e-3, af_snoise takes the centred form in every one of 2,000 seeded
/// cases (a is nonzero), and its mean total width is under half of the Lipschitz form's 2 L R. A
/// mutant that always returns the Lipschitz form has a = 0.
/// </summary>
[Fact]
public void SnoiseChoosesTheCentredFormWhereNarrower()
{
var random = new System.Random(0xAF11);
var sum = 0.0;
for (var b = 0; b < 2000; b++)
{
var centre = new float3(Uniform(random, -20.0f, 20.0f), Uniform(random, -20.0f, 20.0f), Uniform(random, -20.0f, 20.0f));
var axis = UnitVector(random) * 1.0e-3f;
var form = af_snoise(centre, axis, 1.0e-3f);
var lipschitz = 2.0 * SNOISE_LIPSCHITZ * 2.0e-3;
var width = 2.0 * (Math.Abs(form.y) + form.z);
Assert.True(form.y != 0.0f && width < lipschitz, $"af_snoise({centre}, {axis}, 1e-3) = {form}: not the centred form, or no narrower than 2 L R = {lipschitz:R}");
sum += width / lipschitz;
}
output.WriteLine($"AF-REPORT af_snoise at radius 1e-3: mean width over 2 L R {sum / 2000:F3}");
Assert.True(sum / 2000 < 0.5, $"mean width over 2 L R is {sum / 2000:F3}, not under half");
}
/// <summary>
/// af_frustum_ball and af_frustum_axis over iv_frustum_ball's families (NoiseBoundTests.DrawBallCase:
/// the r2 domain and the degenerate, thin, wide-slope, far, from-the-camera, tiny-focal and huge-warp
/// extremes), 2,000 tiles each, warp passed as the warp variation: every exact ray point (m z, z),
/// formed in double from the float32 slope and depth at the footprint's corners at both depth ends
/// (and 56 more for the r2 domain), moved by a full-length outward flow of length warp, is within
/// radius of centre + axis eps(z), eps(z) = clamp((z - z_m) / h, -1, 1) with z_m and h the centre's z
/// and the axis's z (0 when h is 0). The radius is pinned to the derivation both ways:
/// r &lt;= w &lt;= r + 2^-19 (r + |c|_1 + |axis|_1), r = z1 footprintPerDepth + warp in double.
/// </summary>
[Fact]
public void FrustumSliceEnclosesEveryRayPoint()
{
var random = new System.Random(0xAF12);
string[] families = { "r2 domain", "degenerate", "thin", "wide slope", "far", "from the camera", "tiny focal", "huge warp" };
var failures = new Dictionary<string, int>();
var witness = new List<string>();
foreach (var family in families)
{
failures[family] = 0;
for (var t = 0; t < 2000; t++)
{
var (n, focal, slope, z0, z1, warp) = NoiseBoundTests.DrawBallCase(random, family);
var half = n / (2.0f * focal);
var footprint = n / (MathF.Sqrt(2.0f) * focal);
var axis = af_frustum_axis(slope, z0, z1);
var ball = af_frustum_ball(slope, z0, z1, footprint, warp);
var (cx, cy, cz, w) = ((double)ball.x, (double)ball.y, (double)ball.z, (double)ball.w);
var r = z1 * (double)footprint + warp;
var c1 = Math.Abs(cx) + Math.Abs(cy) + Math.Abs(cz);
var a1 = Math.Abs((double)axis.x) + Math.Abs((double)axis.y) + Math.Abs((double)axis.z);
if (!(r <= w && w <= r + Math.ScaleB(r + c1 + a1, -19)))
{
failures[family]++;
if (witness.Count < 8)
witness.Add($"{family} pin: af_frustum_ball({slope}, {z0:R}, {z1:R}, {footprint:R}, {warp:R}) = {ball}, derivation {r:R}");
}
// The axis is (m_c, 1) h exactly, and the radius carries the documented widening, the
// rounding of the centre, the axis (spread) and the radius at 2^-20 each, to 15% (the
// radius rounds too).
var h = (z1 - z0) * 0.5f;
var spread = ((double)Math.Abs(slope.x) + Math.Abs(slope.y) + 1.0) * h;
var widening = Math.ScaleB(c1 + spread + r, -20);
if (!(axis.z == h && axis.x == slope.x * h && axis.y == slope.y * h && Math.Abs(w - r - widening) <= 0.15 * widening + Math.ScaleB(w, -22)))
{
failures[family]++;
if (witness.Count < 8)
witness.Add($"{family} axis or widening: af_frustum_ball({slope}, {z0:R}, {z1:R}, {footprint:R}, {warp:R}) = {ball}, axis {axis}, widening {w - r:R} against {widening:R}");
}
for (var i = 0; i < (family == "r2 domain" ? 64 : 8); i++)
{
var offset = i < 8
? new float2((i & 1) == 0 ? -half : half, (i & 2) == 0 ? -half : half)
: new float2(Uniform(random, -half, half), Uniform(random, -half, half));
var z = i < 8 ? ((i & 4) == 0 ? z0 : z1) : Uniform(random, z0, z1);
var m = slope + offset;
var eps = axis.z == 0.0f ? 0.0 : Math.Clamp((z - cz) / axis.z, -1.0, 1.0);
var (qx, qy, qz) = (cx + axis.x * eps, cy + axis.y * eps, cz + axis.z * eps);
var (dx, dy, dz) = (m.x * (double)z - qx, m.y * (double)z - qy, z - qz);
var away = Math.Sqrt(dx * dx + dy * dy + dz * dz);
var distance = away + warp;
if (!(distance <= w))
{
failures[family]++;
if (witness.Count < 8)
witness.Add($"{family}: af_frustum_ball({slope}, {z0:R}, {z1:R}, {footprint:R}, {warp:R}) = {ball}, axis {axis} misses slope {m}, depth {z:R} by {(distance - w) / w:R} r");
}
}
}
}
Assert.True(failures.Values.Sum() == 0, string.Join(", ", families.Select(f => $"{f} {failures[f]}")) + Environment.NewLine + string.Join(Environment.NewLine, witness));
}
/// <summary>
/// Underflow. Operands of magnitude 1e-21 to 1e-17 put products near and under 2^-126, where float32
/// rounds to a subnormal and the relative absorption S 2^-20 is smaller than the loss: every op
/// still encloses, no tolerance. af_mul and af_scale carry the 2^-126 term; a mutant that subtracts
/// it misses here.
/// </summary>
[Fact]
public void TinyOperandsEncloseTheirPoints()
{
var failures = new List<string>();
foreach (var op in Ops)
{
var (misses, witness) = Harness(op, 0xAF05, 10_000, tiny: true);
if (misses > 0)
failures.Add($"{op.Name}: {misses} misses" + Environment.NewLine + string.Join(Environment.NewLine, witness));
}
Assert.True(failures.Count == 0, string.Join(Environment.NewLine, failures));
}
private const float MinNormal = 1.17549435e-38f;
private static bool Near(double observed, double expected, double relative = 1.0e-5) =>
Math.Abs(observed - expected) <= relative * Math.Abs(expected);
/// <summary>
/// The rounding absorption of each op, pinned both ways. With no propagated error in e, af_add,
/// af_sub, af_scale, af_mul and af_from_iv return an e equal to their documented absorption (and
/// 2^-126 where documented) to a part in 1e5 of its double value, and af_range moves each endpoint
/// outward by (|x0| + |a| + e) 2^-20 to 30% (the endpoints round once). A factor, a sign or an
/// operator changed in any of them fails here; the one-ulp changes of the worst-case absorption are
/// the OneUlpSurvivors set EveryRoundingTermIsCaught names.
/// </summary>
[Fact]
public void AbsorptionsArePinned()
{
var u20 = Math.ScaleB(1.0, -20);
var random = new System.Random(0xAF06);
for (var t = 0; t < 2_000; t++)
{
var (x0, y0) = (Uniform(random, -20.0f, 20.0f), Uniform(random, -20.0f, 20.0f));
var (ax, ay) = (Signed(random, LogUniform(random, 1.0e-3f, 20.0f)), Signed(random, LogUniform(random, 1.0e-3f, 20.0f)));
var scale = Signed(random, LogUniform(random, 1.0e-2f, 4.0f));
var x = new float3(x0, ax, 0.0f);
var y = new float3(y0, ay, 0.0f);
var sum = af_add(x, y);
Assert.True(Near(sum.z, (Math.Abs((double)sum.x) + Math.Abs((double)sum.y)) * u20), $"af_add({x}, {y}) = {sum}: e is not (|x0| + |a|) 2^-20");
var difference = af_sub(x, y);
Assert.True(Near(difference.z, (Math.Abs((double)difference.x) + Math.Abs((double)difference.y)) * u20), $"af_sub({x}, {y}) = {difference}: e is not (|x0| + |a|) 2^-20");
var scaled = af_scale(x, scale);
Assert.True(Near(scaled.z, (Math.Abs((double)scaled.x) + Math.Abs((double)scaled.y)) * u20 + MinNormal), $"af_scale({x}, {scale:R}) = {scaled}: e is not (|x0| + |a|) 2^-20 + 2^-126");
// a_x = 0 leaves the product's propagated error 0 and its a = x0 a_y.
var px = new float3(x0, 0.0f, 0.0f);
var product = af_mul(px, y);
var terms = Math.Abs((double)product.x) + Math.Abs((double)x0 * ay);
Assert.True(Near(product.z, terms * u20 + MinNormal), $"af_mul({px}, {y}) = {product}: e is not (|x0 y0| + |x0 a_y|) 2^-20 + 2^-126");
var iv = DrawInterval(random);
var form = af_from_iv(iv);
var (mid, half) = (((double)iv.x + iv.y) / 2.0, ((double)iv.y - iv.x) / 2.0);
var expected = half + (Math.Abs(mid) + half) * u20 + MinNormal;
Assert.True(Math.Abs(form.z - expected) <= 1.0e-5 * (Math.Abs(mid) + half) * u20 + 1.0e-6 * half + 1.0e-12 * expected,
$"af_from_iv({iv}) = {form}: e is not half-width + (|mid| + half-width) 2^-20 + 2^-126 = {expected:R}");
var wide = new float3(Uniform(random, -20.0f, 20.0f), Signed(random, LogUniform(random, 1.0e-3f, 20.0f)), LogUniform(random, 1.0e-3f, 20.0f));
var range = af_range(wide);
var reach = Math.Abs((double)wide.y) + wide.z;
var absorption = (Math.Abs((double)wide.x) + reach) * u20;
var outward = ((double)range.y - range.x) / 2.0 - reach;
Assert.True(Math.Abs(outward - absorption) <= 0.3 * absorption, $"af_range({wide}) = {range}: endpoints are {outward:R} past x0 -+ (|a| + e), not (|x0| + |a| + e) 2^-20 = {absorption:R}");
}
}
/// <summary>
/// af_from_iv and af_range at the float32 extremes. Subnormal endpoints (including [0, the least
/// subnormal], whose halved width underflows to 0) and endpoints near the float maximum (where lo +
/// hi, hi - lo or the absorption would overflow) give a form whose set holds both endpoints and the
/// midpoint, and an af_range that holds them with no NaN endpoint. A form with an infinite or NaN
/// component, as an overflowed op returns, gives an af_range with no NaN endpoint.
/// </summary>
[Fact]
public void FromIntervalAndRangeHoldAtTheExtremes()
{
var least = BitConverter.Int32BitsToSingle(1);
var max = float.MaxValue;
var random = new System.Random(0xAF07);
var specials = new[] { 0.0f, least, 2 * least, 100 * least, MinNormal - least, MinNormal, 3.0f * MinNormal, 1.0e-30f, 1.0f, 1.0e20f, 1.0e30f, max / 4, max / 2, max * 0.75f, max };
var misses = new List<string>();
void Check(float lo, float hi)
{
var form = af_from_iv(new float2(lo, hi));
var range = af_range(form);
if (float.IsNaN(form.x) || float.IsNaN(form.z) || float.IsNaN(range.x) || float.IsNaN(range.y))
misses.Add($"af_from_iv([{lo:R}, {hi:R}]) = {form}, af_range = {range}: NaN");
foreach (var value in new[] { lo, hi, (float)(((double)lo + hi) / 2.0) })
{
if (!InForm(form, 0.0f, value))
misses.Add($"af_from_iv([{lo:R}, {hi:R}]) = {form} misses {value:R}");
if (!(range.x <= value && value <= range.y))
misses.Add($"af_range(af_from_iv([{lo:R}, {hi:R}])) = {range} misses {value:R}");
}
// Tight where the answer is representable: an endpoint is no farther than 2^-18 (|lo| + |hi|)
// from the interval's, unless that bound passes the float maximum, where infinity is the answer.
var slack = Math.ScaleB(Math.Abs((double)lo) + Math.Abs((double)hi), -18) + 1.0e-37;
if ((double)lo - slack > -(double)max && range.x < (double)lo - slack)
misses.Add($"af_range(af_from_iv([{lo:R}, {hi:R}])).x = {range.x:R}, more than 2^-18 (|lo| + |hi|) under lo");
if ((double)hi + slack < (double)max && range.y > (double)hi + slack)
misses.Add($"af_range(af_from_iv([{lo:R}, {hi:R}])).y = {range.y:R}, more than 2^-18 (|lo| + |hi|) over hi");
}
foreach (var a in specials)
foreach (var b in specials)
{
Check(-a, b);
Check(-Math.Max(a, b), -Math.Min(a, b));
Check(Math.Min(a, b), Math.Max(a, b));
}
for (var t = 0; t < 200_000; t++)
{
var sign = random.Next(2) == 0 ? 1 : -1;
var lo = sign * BitConverter.Int32BitsToSingle(random.Next(0, random.Next(2) == 0 ? 1 << 20 : 0x00800000));
var width = BitConverter.Int32BitsToSingle(random.Next(0, random.Next(2) == 0 ? 1 << 12 : 1 << 22));
Check(lo, lo + width);
}
var odd = new[] { float.PositiveInfinity, float.NegativeInfinity, float.NaN, 0.0f, 1.0f, max };
foreach (var x0 in odd)
foreach (var a in odd)
foreach (var e in odd)
{
var range = af_range(new float3(x0, a, e));
// The range encloses x0 -+ (|a| + e) wherever that is a number, and is the whole line on the
// side where it is not (infinity minus infinity, a NaN component): the infinity on that side.
var spread = Math.Abs((double)a) + e;
var enclosesLow = double.IsNaN(x0 - spread) ? range.x == float.NegativeInfinity : range.x <= x0 - spread;
var enclosesHigh = double.IsNaN(x0 + spread) ? range.y == float.PositiveInfinity : range.y >= x0 + spread;
if (!enclosesLow || !enclosesHigh)
misses.Add($"af_range(({x0}, {a}, {e})) = {range}");
}
Assert.True(misses.Count == 0, $"{misses.Count} misses:" + Environment.NewLine + string.Join(Environment.NewLine, misses.Take(8)));
}
/// <summary>
/// af_snoise's width, pinned both ways. Over 6,000 seeded (centre, axis, radius) with reach R from 1e-3
/// to 3 (centres in [-20, 20]^3, a third out to [-1000, 1000]^3, a third of the axes and a third of the
/// radii 0), the form is the documented rule evaluated in double from snoise_grad: x0 = n; the centred
/// form when |g . axis| + |g| radius + M R^2 / 2 is under L R (checked away from that boundary), else
/// the Lipschitz form; e = that half-width + the float32 allowance (2^-14 + 2^-12 R + L 2^-20 (|centre|_1 + R))
/// plus the form's own rounding absorption (|n| + |a| + e) 2^-20; each of x0, a and e within a part in
/// 1e6 of it. The enclosure tests pass on an allowance or an absorption of any size on the CPU, where
/// float32 is exact; this is the test that sees them.
/// </summary>
[Fact]
public void SnoiseFormCarriesItsDocumentedWidth()
{
var u20 = Math.ScaleB(1.0, -20);
double lipschitz = SNOISE_LIPSCHITZ, hessian = SNOISE_HESSIAN;
var random = new System.Random(0xAF13);
var (centred, wide) = (0, 0);
for (var t = 0; t < 6_000; t++)
{
var spread = t % 3 == 2 ? 1000.0f : 20.0f;
var centre = new float3(Uniform(random, -spread, spread), Uniform(random, -spread, spread), Uniform(random, -spread, spread));
var axis = random.Next(3) == 0 ? new float3(0.0f, 0.0f, 0.0f) : UnitVector(random) * LogUniform(random, 1.0e-3f, 3.0f);
var radius = random.Next(3) == 0 ? 0.0f : LogUniform(random, 1.0e-3f, 3.0f);
var g = snoise_grad(centre);
var reach = (double)length(axis) + radius;
var lip = lipschitz * reach;
var a = (double)g.x * axis.x + (double)g.y * axis.y + (double)g.z * axis.z;
var remainder = Math.Sqrt((double)g.x * g.x + (double)g.y * g.y + (double)g.z * g.z) * radius + 0.5 * hessian * reach * reach;
var allowance = Math.ScaleB(1.0, -14) + Math.ScaleB(reach, -12)
+ lipschitz * u20 * (Math.Abs((double)centre.x) + Math.Abs((double)centre.y) + Math.Abs((double)centre.z) + reach);
double Width(bool useCentred, out double slope)
{
slope = useCentred ? a : 0.0;
var e = (useCentred ? remainder : lip) + allowance;
return e + (Math.Abs((double)g.w) + Math.Abs(slope) + e) * u20;
}
var zCentred = Width(true, out var slopeCentred);
var zLipschitz = Width(false, out _);
var form = af_snoise(centre, axis, radius);
var observedCentred = Math.Abs(form.z - zCentred) <= Math.Abs(form.z - zLipschitz);
var expectedCentred = Math.Abs(a) + remainder < lip;
if (Math.Abs(Math.Abs(a) + remainder - lip) > 1.0e-3 * lip)
Assert.True(observedCentred == expectedCentred || zCentred == zLipschitz, $"af_snoise({centre}, {axis}, {radius:R}) = {form} took the {(observedCentred ? "centred" : "Lipschitz")} form; the rule takes the {(expectedCentred ? "centred" : "Lipschitz")} one (|a| + remainder = {Math.Abs(a) + remainder:R}, L R = {lip:R})");
var z = observedCentred ? zCentred : zLipschitz;
var slope = observedCentred ? slopeCentred : 0.0;
Assert.True(form.x == g.w, $"af_snoise({centre}, {axis}, {radius:R}).x0 = {form.x:R}, snoise = {g.w:R}");
Assert.True(Math.Abs(form.y - slope) <= 1.0e-5 * Math.Abs((double)length(new float3(g.x, g.y, g.z)) * length(axis)) + 1.0e-12, $"af_snoise({centre}, {axis}, {radius:R}).a = {form.y:R}, g . axis = {slope:R}");
Assert.True(Math.Abs(form.z - z) <= 1.0e-6 * z, $"af_snoise({centre}, {axis}, {radius:R}) = {form}: e is {form.z:R}, the documented rule gives {z:R}");
centred += observedCentred ? 1 : 0;
wide += !observedCentred ? 1 : 0;
}
output.WriteLine($"AF-REPORT af_snoise width: {centred} centred cases, {wide} Lipschitz cases");
Assert.True(centred > 1000 && wide > 1000, $"the cases do not exercise both forms: {centred} centred, {wide} Lipschitz");
}
}

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

@ -163,7 +163,7 @@ public sealed class GlslMirrorTests
}
if (want["check"] is not null)
AssertFrustumBall(want["cases"]!.AsArray());
output.WriteLine($"{name}: {Families[f].Verify!(want["cases"]!.AsArray())}");
}
committed.Remove("functions");
@ -188,8 +188,9 @@ public sealed class GlslMirrorTests
// FrustumBall applied to C#'s own results: every case's ball encloses the exact one and is no wider
// than its rounding widening allows.
private static void AssertFrustumBall(JsonArray cases)
private static string AssertFrustumBall(JsonArray cases)
{
double low = double.PositiveInfinity, high = double.PositiveInfinity;
foreach (var entry in cases)
{
var sides = ((string)entry!).Split(" -> ");
@ -203,7 +204,148 @@ public sealed class GlslMirrorTests
var c1 = Math.Abs(mx * zm) + Math.Abs(my * zm) + Math.Abs(zm);
Assert.True(b[3] >= reach, $"iv_frustum_ball case {entry}: radius short of the exact ball by {reach - b[3]}");
Assert.True(b[3] <= reach + Math.ScaleB(radius + c1, -19), $"iv_frustum_ball case {entry}: radius past the exact ball by {b[3] - reach}, more than 2^-19 (r + |c|_1)");
low = Math.Min(low, (b[3] - reach) / radius);
high = Math.Min(high, (reach + Math.ScaleB(radius + c1, -19) - b[3]) / Math.ScaleB(radius + c1, -19));
}
return $"least slack: {low:E3} of r above the exact ball, {high:F3} of 2^-19 (r + |c|_1) below its widening";
}
private static readonly double[] CheckEps = { -1.0, -0.5, 0.0, 0.5, 1.0 };
// The origin, the six axis directions and the eight cube corners, normalised.
private static readonly (double X, double Y, double Z)[] CheckDirections = CheckDirectionList().ToArray();
private static IEnumerable<(double X, double Y, double Z)> CheckDirectionList()
{
yield return (0.0, 0.0, 0.0);
foreach (var s in new[] { -1.0, 1.0 })
{
yield return (s, 0.0, 0.0);
yield return (0.0, s, 0.0);
yield return (0.0, 0.0, s);
}
var corner = 1.0 / Math.Sqrt(3.0);
foreach (var sx in new[] { -corner, corner })
foreach (var sy in new[] { -corner, corner })
foreach (var sz in new[] { -corner, corner })
yield return (sx, sy, sz);
}
// SampledForm applied to C#'s own results, with C#'s own snoise or fBm as the reference value. A form
// must enclose its reference at every sample, and be no wider than the Lipschitz form's allowance.
// allowance is the part of e that is the float32 allowance; a family that allowanceDecides must have a
// sample that falls outside the form once it is removed, so the family can tell a device whose error
// passes the allowance from one that does not.
private static string AssertSampledForm(JsonArray cases, string name, Func<string[], float3, float> reference, Func<string[], double[], double> widthBound,
Func<string[], double[], double> allowance, bool allowanceDecides)
{
double least = double.PositiveInfinity, leastRelative = double.PositiveInfinity, tightest = double.PositiveInfinity;
var decided = 0;
foreach (var entry in cases)
{
var sides = ((string)entry!).Split(" -> ");
var tokens = sides[0].Split(' ');
double[] a = Floats(string.Join(' ', tokens.Take(7))), form = Floats(sides[1]);
double x0 = form[0], slope = form[1], e = form[2];
var withoutAllowance = e - allowance(tokens, a);
var allowanceNeeded = false;
foreach (var eps in CheckEps)
foreach (var d in CheckDirections)
{
var p = new float3((float)(a[0] + a[3] * eps + a[6] * d.X), (float)(a[1] + a[4] * eps + a[6] * d.Y), (float)(a[2] + a[5] * eps + a[6] * d.Z));
double value = reference(tokens, p), middle = x0 + slope * eps;
var inside = Math.Min(value - (middle - e), middle + e - value);
Assert.True(inside >= 0.0, $"{name} case {entry}: the reference {value:R} at eps {eps} and direction {d} is outside the form's range [{middle - e:R}, {middle + e:R}] by {-inside:R}");
least = Math.Min(least, inside);
leastRelative = Math.Min(leastRelative, inside / e);
allowanceNeeded |= Math.Abs(value - middle) > withoutAllowance;
}
if (allowanceNeeded)
decided++;
var bound = widthBound(tokens, a);
Assert.True(Math.Abs(slope) + e <= bound + Math.ScaleB(1.0 + bound, -16), $"{name} case {entry}: the width {Math.Abs(slope) + e:R} is past B + 2^-16 (1 + B) = {bound + Math.ScaleB(1.0 + bound, -16):R}");
tightest = Math.Min(tightest, (bound + Math.ScaleB(1.0 + bound, -16) - Math.Abs(slope) - e) / bound);
}
Assert.True(!allowanceDecides || decided > 0, $"{name}: no sample falls outside the form once the float32 allowance is removed, so the family cannot fail a device whose error passes it");
return $"least enclosure slack {least:E3} ({leastRelative:E3} of e), least tightness slack {tightest:E3} of B, {decided} of {cases.Count} cases need the allowance";
}
// The float32 allowance of one octave of af_snoise: 2^-14 + 2^-12 R + L 2^-20 (|c|_1 + R).
private static double SnoiseAllowance(double centre1, double reach) =>
Math.ScaleB(1.0, -14) + Math.ScaleB(reach, -12) + math.SNOISE_LIPSCHITZ * Math.ScaleB(centre1 + reach, -20);
// The Lipschitz form's width bound for one octave of the sampled-form check.
private static double SnoiseWidthBound(double centre1, double reach) => math.SNOISE_LIPSCHITZ * reach + SnoiseAllowance(centre1, reach);
private static double Length3(double[] a, int first) => Math.Sqrt(a[first] * a[first] + a[first + 1] * a[first + 1] + a[first + 2] * a[first + 2]);
private static double SnoiseOctave(Func<double, double, double> octave, double[] a, double frequency) =>
octave((Math.Abs(a[0]) + Math.Abs(a[1]) + Math.Abs(a[2])) * frequency, (Length3(a, 3) + a[6]) * frequency);
// The amplitude-weighted sum of a per-octave quantity over the fBm's octaves.
private static double FbmOctaves(string[] t, double[] a, Func<double, double, double> octave)
{
var (octaves, lacunarity, gain) = ((int)Convert.ToUInt32(t[7], 16), Floats(t[8])[0], Floats(t[9])[0]);
double sum = 0.0, amplitude = 1.0, frequency = 1.0;
for (var i = 0; i < Math.Clamp(octaves, 0, 16); i++)
{
sum += amplitude * SnoiseOctave(octave, a, frequency);
amplitude *= gain;
frequency *= lacunarity;
}
return sum;
}
private static string AssertSnoiseForm(JsonArray cases, bool allowanceDecides = false) =>
AssertSampledForm(cases, "af_snoise", (_, p) => math.snoise(p),
(_, a) => SnoiseOctave(SnoiseWidthBound, a, 1.0), (_, a) => SnoiseOctave(SnoiseAllowance, a, 1.0), allowanceDecides);
private static string AssertFbmForm(JsonArray cases, bool allowanceDecides = false) =>
AssertSampledForm(cases, "af_fbm",
(t, p) => math.fbm_grad(p, (int)Convert.ToUInt32(t[7], 16), (float)Floats(t[8])[0], (float)Floats(t[9])[0]).w,
(t, a) => FbmOctaves(t, a, SnoiseWidthBound), (t, a) => FbmOctaves(t, a, SnoiseAllowance), allowanceDecides);
// AffineFrustumBall applied to C#'s own results.
private static string AssertAffineFrustumBall(JsonArray cases)
{
double enclosure = double.PositiveInfinity, lower = double.PositiveInfinity, upper = double.PositiveInfinity;
var diagonal = new[] { (0.0, 0.0), (1.0, 0.0), (-1.0, 0.0), (0.0, 1.0), (0.0, -1.0), (Math.Sqrt(0.5), Math.Sqrt(0.5)), (Math.Sqrt(0.5), -Math.Sqrt(0.5)), (-Math.Sqrt(0.5), Math.Sqrt(0.5)), (-Math.Sqrt(0.5), -Math.Sqrt(0.5)) };
foreach (var entry in cases)
{
var sides = ((string)entry!).Split(" -> ");
var a = Floats(sides[0]);
var b = Floats(sides[1]);
double mx = a[0], my = a[1], z0 = a[2], z1 = a[3], fp = a[4], warp = a[5];
double zm = (z0 + z1) / 2.0, h = (z1 - z0) / 2.0;
var r0 = z1 * fp + warp;
var c1 = Math.Abs(mx * zm) + Math.Abs(my * zm) + Math.Abs(zm);
var spread = (Math.Abs(mx) + Math.Abs(my) + 1.0) * h;
double dx = b[0] - mx * zm, dy = b[1] - my * zm, dz = b[2] - zm;
var d = Math.Sqrt(dx * dx + dy * dy + dz * dz);
foreach (var z in new[] { z0, (3.0 * z0 + z1) / 4.0, zm, (z0 + 3.0 * z1) / 4.0, z1 })
{
var eps = h > 0.0 ? Math.Clamp((z - zm) / h, -1.0, 1.0) : 0.0;
foreach (var (ex, ey) in diagonal)
{
double px = (mx + fp * ex) * z - (b[0] + mx * h * eps), py = (my + fp * ey) * z - (b[1] + my * h * eps), pz = z - (b[2] + h * eps);
var slack = b[3] - (Math.Sqrt(px * px + py * py + pz * pz) + warp);
Assert.True(slack >= 0.0, $"af_frustum_ball case {entry}: the ray point at depth {z:R} and offset ({ex}, {ey}) is outside the ball by {-slack:R}");
enclosure = Math.Min(enclosure, slack / b[3]);
}
}
Assert.True(b[3] >= r0 + d, $"af_frustum_ball case {entry}: radius short of R0 + d by {r0 + d - b[3]}");
Assert.True(b[3] <= r0 + d + Math.ScaleB(r0 + c1 + spread, -19), $"af_frustum_ball case {entry}: radius past R0 + d by {b[3] - r0 - d}, more than 2^-19 (R0 + |c0|_1 + spread)");
lower = Math.Min(lower, (b[3] - r0 - d) / r0);
upper = Math.Min(upper, (r0 + d + Math.ScaleB(r0 + c1 + spread, -19) - b[3]) / Math.ScaleB(r0 + c1 + spread, -19));
}
return $"least enclosure slack {enclosure:E3} of the radius, least slack above R0 + d {lower:E3} of R0, below its widening {upper:F3} of 2^-19 (R0 + |c0|_1 + spread)";
}
// A committed case's arguments, decoded into the C# method's parameter types in Hex's field order.
@ -249,7 +391,7 @@ public sealed class GlslMirrorTests
// phacelle's sums of exp, sin and cos cancel; Windows differs from Linux by up to 3.6e-7 of it.
private const double ScaledPlatformBound = 9.5367431640625e-7;
private sealed record Family(string Glsl, string Tolerance, Func<System.Random, object[]> Arguments, string? Check = null, Platform Platform = Platform.None);
private sealed record Family(string Glsl, string Tolerance, Func<System.Random, object[]> Arguments, string? Check = null, Platform Platform = Platform.None, Func<JsonArray, string>? Verify = null);
private static float Uniform(System.Random r, float lo, float hi) => lo + (hi - lo) * r.NextSingle();
private static float3 Point3(System.Random r) => new(Uniform(r, -50.0f, 50.0f), Uniform(r, -50.0f, 50.0f), Uniform(r, -50.0f, 50.0f));
@ -291,6 +433,76 @@ public sealed class GlslMirrorTests
"and d = |GPU centre - c| (Euclidean); every case must have r + d <= GPU radius <= r + d + 2^-19 (r + |c|_1). " +
"Report the largest ulp distance per component as for any ulp-bounded family; only the two bounds fail the check.";
// An affine form (x0, a, e): e >= 0, spanning the widths a probe's forms take.
private static float3 Form(System.Random r) => new(Uniform(r, -20.0f, 20.0f), Uniform(r, -5.0f, 5.0f), MathF.Pow(10.0f, Uniform(r, -4.0f, 1.0f)));
// A slice's axis and radius in noise space, log-uniform from 1e-3 to 0.5: below about 0.15 the
// centred form of af_snoise wins, above it the Lipschitz form does, so both are drawn.
private static float3 Axis3(System.Random r)
{
var scale = MathF.Pow(10.0f, Uniform(r, -3.0f, -0.3f));
return new float3(Uniform(r, -1.0f, 1.0f) * scale, Uniform(r, -1.0f, 1.0f) * scale, Uniform(r, -1.0f, 1.0f) * scale);
}
private static float BallRadius(System.Random r) => MathF.Pow(10.0f, Uniform(r, -3.0f, -0.3f));
// A slice of reach 1e-6 to 1e-3, the short slices and small footprints near a camera, where the float32
// allowance is most of the form's width: the centre in [-50, 50] or [-1, 1], an axis of random direction
// (zero in one case of four) and a radius, each log-uniform in 1e-6 to 1e-3.
private static object[] SmallReachSlice(System.Random r)
{
var span = r.Next(2) == 0 ? 50.0f : 1.0f;
var centre = new float3(Uniform(r, -span, span), Uniform(r, -span, span), Uniform(r, -span, span));
var axis = new float3(0.0f, 0.0f, 0.0f);
if (r.Next(4) != 0)
{
var direction = math.normalize(new float3(Uniform(r, -1.0f, 1.0f), Uniform(r, -1.0f, 1.0f), Uniform(r, -1.0f, 1.0f)));
axis = direction * MathF.Pow(10.0f, Uniform(r, -6.0f, -3.0f));
}
return new object[] { centre, axis, MathF.Pow(10.0f, Uniform(r, -6.0f, -3.0f)) };
}
private const string SampledForm =
"the arguments are (cx, cy, cz, ax, ay, az, radius, ...) and the result is the form (x0, a, e). Enclosure: for eps in " +
"{-1, -0.5, 0, 0.5, 1} and d in {0, +-ex, +-ey, +-ez, (+-1, +-1, +-1) / sqrt(3)} (15 directions), form the point " +
"p = c + axis eps + radius d in double and round it to float32; the device's own reference value at p must lie in " +
"[x0 + a eps - e, x0 + a eps + e] with no tolerance, so a device whose gradient error exceeds the allowance fails. " +
"Tightness: with L = 7.9640074 and, for each octave (one for af_snoise), R = |axis| + radius, " +
"B = L R + 2^-14 + 2^-12 R + L 2^-20 (|c|_1 + R), the width |a| + e must be at most B + 2^-16 (1 + B). " +
"Report the least slack of each side, absolutely and as a fraction of e; only the two bounds fail the check. ";
private const string SnoiseForm = SampledForm + "The reference value is cultmath_snoise(p).";
private const string FbmForm = SampledForm +
"The reference value is cultmath_fbm_grad(p, octaves, lacunarity, gain).w, and B is the sum over octaves i of gain^i B_i, " +
"where B_i is B for the centre lacunarity^i c, the axis lacunarity^i axis and the radius lacunarity^i radius.";
// The two small-reach families add the centre leg to every af_snoise and af_fbm family. The device's
// reference value is its own noise, so a device whose noise is wrong everywhere encloses its own wrong
// values; C#'s value at the centre, the form's own x0 in the fixture, is the independent reading.
private const string CentreLeg =
" Centre leg, for every cultmath_af_snoise and cultmath_af_fbm family of this fixture, the earlier ones included: " +
"the fixture's own x0 in the result (C#'s value at the centre) must lie in [x0 - e, x0 + e] of the GPU's form, with " +
"no tolerance, so a device whose noise differs from C#'s by more than its e fails.";
private const string SmallReach =
" These cases have reach |axis| + radius of 1e-6 to 1e-3, where the float32 allowance is most of e and the second order " +
"remainder is not; the enclosure above fails a device whose gradient error exceeds the allowance.";
private const string SnoiseSmallForm = SnoiseForm + SmallReach + CentreLeg;
private const string FbmSmallForm = FbmForm + SmallReach + CentreLeg;
private const string AffineFrustumBall =
"enclosure and tightness: from the arguments (mx, my, z0, z1, fp, warpVariation) compute in double zm = (z0 + z1) / 2, " +
"h = (z1 - z0) / 2, axis = (mx h, my h, h), c0 = (mx zm, my zm, zm), |c0|_1 = |mx zm| + |my zm| + |zm|, R0 = z1 fp + warpVariation, " +
"spread = (|mx| + |my| + 1) h and d = |GPU centre - c0|. Enclosure: for z in {z0, (3 z0 + z1) / 4, zm, (z0 + 3 z1) / 4, z1}, " +
"eps = clamp((z - zm) / h, -1, 1) and the nine footprint offsets e in {(0, 0), (+-1, 0), (0, +-1), (+-1, +-1) / sqrt(2)}, " +
"the point P = ((mx + fp e.x) z, (my + fp e.y) z, z) satisfies |P - (GPU centre + axis eps)| + warpVariation <= GPU radius. " +
"Tightness: R0 + d <= GPU radius <= R0 + d + 2^-19 (R0 + |c0|_1 + spread). " +
"Report the least slack of each side; only the bounds fail the check.";
// A Phacelle stripe wave vector: each component in [-6, 6], about one stripe per cell.
private static float3 Side(System.Random r) => new(Uniform(r, -6.0f, 6.0f), Uniform(r, -6.0f, 6.0f), Uniform(r, -6.0f, 6.0f));
@ -331,12 +543,38 @@ public sealed class GlslMirrorTests
var slope = new float2(Uniform(r, -1.0f, 1.0f), Uniform(r, -1.0f, 1.0f));
var z0 = Uniform(r, 0.1f, 100.0f);
return new object[] { slope, z0, z0 + Uniform(r, 0.01f, 50.0f), Uniform(r, 0.0001f, 0.05f), Uniform(r, 0.0f, 2.0f) };
}, FrustumBall),
}, FrustumBall, Verify: AssertFrustumBall),
new("vec2 cultmath_iv_exp(vec2)", "ulp-bounded, platform exp", r => new object[] { Interval(r) }, Platform: Platform.Ulp),
new("CultPhasor cultmath_phacelle(vec3, vec3, float, float)", "ulp-bounded, platform exp, sin, cos; across platforms |diff| <= 2^-20 max(|v|, 1)",
r => new object[] { Point3(r), Side(r), Uniform(r, 0.0f, 1.0f), 0.5f }, Platform: Platform.Scaled),
new("float cultmath_snoise(vec3)", "ulp-bounded, on the simplex tie set", r => new object[] { TiePoint3(r) }),
new("vec4 cultmath_snoise_grad(vec3)", "ulp-bounded, on the simplex tie set", r => new object[] { TiePoint3(r) }),
new("vec3 cultmath_af_point(float)", "ulp-bounded", r => new object[] { Uniform(r, -20.0f, 20.0f) }),
new("vec3 cultmath_af_symbol(float, float)", "ulp-bounded", r => new object[] { Uniform(r, -20.0f, 20.0f), Uniform(r, -5.0f, 5.0f) }),
new("vec3 cultmath_af_from_iv(vec2)", "ulp-bounded", r => new object[] { Interval(r) }),
new("vec2 cultmath_af_range(vec3)", "ulp-bounded", r => new object[] { Form(r) }),
new("vec3 cultmath_af_add(vec3, vec3)", "ulp-bounded", r => new object[] { Form(r), Form(r) }),
new("vec3 cultmath_af_sub(vec3, vec3)", "ulp-bounded", r => new object[] { Form(r), Form(r) }),
new("vec3 cultmath_af_neg(vec3)", "ulp-bounded", r => new object[] { Form(r) }),
new("vec3 cultmath_af_scale(vec3, float)", "ulp-bounded", r => new object[] { Form(r), Uniform(r, -4.0f, 4.0f) }),
new("vec3 cultmath_af_add_iv(vec3, vec2)", "ulp-bounded", r => new object[] { Form(r), Interval(r) }),
new("vec3 cultmath_af_mul(vec3, vec3)", "ulp-bounded", r => new object[] { Form(r), Form(r) }),
new("vec3 cultmath_af_snoise(vec3, vec3, float)", "ulp-bounded", r => new object[] { Point3(r), Axis3(r), BallRadius(r) }, SnoiseForm, Verify: cases => AssertSnoiseForm(cases)),
new("vec3 cultmath_af_fbm(vec3, vec3, float, int, float, float)", "ulp-bounded", r => new object[] { Point3(r), Axis3(r), BallRadius(r), 4, 2.0f, 0.5f }, FbmForm, Verify: cases => AssertFbmForm(cases)),
new("vec3 cultmath_af_frustum_axis(vec2, float, float)", "ulp-bounded", r =>
{
var z0 = Uniform(r, 0.1f, 100.0f);
return new object[] { new float2(Uniform(r, -1.0f, 1.0f), Uniform(r, -1.0f, 1.0f)), z0, z0 + Uniform(r, 0.01f, 50.0f) };
}),
new("vec4 cultmath_af_frustum_ball(vec2, float, float, float, float)", "ulp-bounded", r =>
{
var slope = new float2(Uniform(r, -1.0f, 1.0f), Uniform(r, -1.0f, 1.0f));
var z0 = Uniform(r, 0.1f, 100.0f);
return new object[] { slope, z0, z0 + Uniform(r, 0.01f, 50.0f), Uniform(r, 0.0001f, 0.05f), Uniform(r, 0.0f, 2.0f) };
}, AffineFrustumBall, Verify: AssertAffineFrustumBall),
new("vec2 cultmath_iv_snoise_ball(vec2, float)", "ulp-bounded", r => new object[] { new float2(Uniform(r, -50.0f, 50.0f), Uniform(r, -50.0f, 50.0f)), Uniform(r, 0.001f, 2.0f) }),
new("vec3 cultmath_af_snoise(vec3, vec3, float)", "ulp-bounded, reach 1e-6 to 1e-3", SmallReachSlice, SnoiseSmallForm, Verify: cases => AssertSnoiseForm(cases, allowanceDecides: true)),
new("vec3 cultmath_af_fbm(vec3, vec3, float, int, float, float)", "ulp-bounded, reach 1e-6 to 1e-3", r => SmallReachSlice(r).Concat(new object[] { 4, 2.0f, 0.5f }).ToArray(), FbmSmallForm, Verify: cases => AssertFbmForm(cases, allowanceDecides: true)),
};
// The fixture text, each case's arguments given by family and point index.

1
packages/cultmath/tests/CultMath.Tests/HlslMirrorTests.cs
View file

@ -66,6 +66,7 @@ public sealed class HlslMirrorTests
[InlineData("CultMath.hlsl")]
[InlineData("CultMath.Phacelle.hlsl")]
[InlineData("CultMath.Interval.hlsl")]
[InlineData("CultMath.Affine.hlsl")]
public void UnityPackageShaderIsIdenticalToTheCanonicalShader(string name)
{
var root = GetCultMathRoot();

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

@ -33,6 +33,10 @@ public sealed class HlslSourceCompatibilityTests
Assert.NotNull(ShaderFunction(shader, "cultmath_phacelle", typeof(float3), typeof(float3), typeof(float), typeof(float)));
Assert.NotNull(ShaderFunction(shader, "cultmath_iv_snoise_ball", typeof(float3), typeof(float)));
Assert.NotNull(ShaderFunction(shader, "cultmath_iv_fbm_ball", typeof(float3), typeof(float), typeof(int), typeof(float), typeof(float)));
Assert.NotNull(ShaderFunction(shader, "cultmath_iv_snoise_ball", typeof(float2), typeof(float)));
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)));
}
// HLSL functions carry no access modifier, so they compile as private instance methods.
@ -133,7 +137,7 @@ public sealed class HlslSourceCompatibilityTests
}
Assert.True(mismatches.Count == 0, $"{mismatches.Count} mismatches:{Environment.NewLine}{string.Join(Environment.NewLine, mismatches.Take(12))}");
Assert.Equal(53, compared.Count);
Assert.Equal(67, compared.Count);
}
/// <summary>
@ -165,6 +169,81 @@ public sealed class HlslSourceCompatibilityTests
}
}
/// <summary>
/// The generic comparison draws centre, axis and radius from [-100, 100], where the reach is far past the
/// 0.15 or so below which <c>cultmath_af_snoise</c>'s centred form is the narrower one, so an HLSL edit
/// of the choice between the centred and the Lipschitz form (a threshold, a reordered sum) survives it.
/// This sweep draws realistic slices (axis and radius log-uniform in 1e-3 to 0.5, the reach of a tile
/// near a camera), which land on both sides of the choice and next to it, for <c>cultmath_af_snoise</c>
/// and <c>cultmath_af_fbm</c>.
/// </summary>
[Fact]
public void AffineSnoiseMirrorMatchesCSharpBitForBitOnBothSidesOfItsBranch()
{
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 snoiseMirror = ShaderFunction(shaderType, "cultmath_af_snoise", typeof(float3), typeof(float3), typeof(float))!;
var fbmMirror = ShaderFunction(shaderType, "cultmath_af_fbm", typeof(float3), typeof(float3), typeof(float), typeof(int), typeof(float), typeof(float))!;
var random = new System.Random(0x5EED3);
float Next(float extent) => (random.NextSingle() * 2.0f - 1.0f) * extent;
float Reach() => 1.0e-3f * MathF.Pow(500.0f, random.NextSingle());
int centred = 0, lipschitz = 0;
for (var i = 0; i < 20000; i++)
{
var centre = new float3(Next(50.0f), Next(50.0f), Next(50.0f));
var axis = normalize(new float3(Next(1.0f), Next(1.0f), Next(1.0f))) * Reach();
var radius = Reach();
var expected = af_snoise(centre, axis, radius);
var actual = Values(snoiseMirror.Invoke(shader, new object[] { centre, axis, radius })!);
Assert.True(BitwiseEqual(Values(expected), actual), $"cultmath_af_snoise differs from C# at centre = {centre}, axis = {axis}, radius = {radius}");
if (expected.y != 0.0f)
centred++;
else
lipschitz++;
if (i % 8 == 0)
{
var fbmExpected = Values(af_fbm(centre, axis, radius, 4, 2.0f, 0.5f));
var fbmActual = Values(fbmMirror.Invoke(shader, new object[] { centre, axis, radius, 4, 2.0f, 0.5f })!);
Assert.True(BitwiseEqual(fbmExpected, fbmActual), $"cultmath_af_fbm differs from C# at centre = {centre}, axis = {axis}, radius = {radius}");
}
}
Assert.True(centred > 2000 && lipschitz > 2000, $"the sweep meets one branch too rarely: centred {centred}, Lipschitz {lipschitz} of 20000");
}
/// <summary>
/// The generic comparison stops at 1e7, so an HLSL edit of <c>cultmath_af_range</c>'s NaN test that
/// misreads the finite values near the float maximum (or an overflowed endpoint) survives it. Forms
/// from zero to the float maximum, infinite and NaN components, both signs, in every slot: the mirror
/// returns C#'s bits, and neither endpoint is ever NaN.
/// </summary>
[Fact]
public void AffineRangeMirrorMatchesCSharpBitForBitAtTheFloatExtremes()
{
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_af_range", typeof(float3))!;
var magnitudes = new[] { 0.0f, 1.0f, 1.0e30f, 1.0e38f, float.MaxValue / 2.0f, float.MaxValue, float.PositiveInfinity, float.NaN };
var values = magnitudes.Concat(magnitudes.Where(m => !float.IsNaN(m)).Select(m => -m)).ToArray();
foreach (var x0 in values)
foreach (var a in values)
foreach (var e in values)
{
var form = new float3(x0, a, e);
var expected = af_range(form);
var actual = Values(mirror.Invoke(shader, new object[] { form })!);
Assert.True(BitwiseEqual(Values(expected), actual), $"cultmath_af_range differs from C# at {form}");
Assert.False(float.IsNaN(expected.x) || float.IsNaN(expected.y), $"af_range({form}) is {expected}, a NaN endpoint");
}
}
/// <summary>
/// The dense sweep almost never meets a cell whose box bound is just under the 2.25 cut and whose
/// weight is still nonzero: a mirror cut of 2.245, 2.24 or 2.2 survives it. These points sit next

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

File diff suppressed because it is too large Load diff

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

File diff suppressed because it is too large Load diff

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

@ -0,0 +1,132 @@
#ifndef CULTMATH_AFFINE_HLSL
#define CULTMATH_AFFINE_HLSL
// Reduced affine forms with one shared symbol over float3(x0, a, e): the set { x0 + a eps + e delta :
// delta in [-1, 1] } at the region's shared eps in [-1, 1], e >= 0. Every cultmath_af_* function mirrors
// math.Affine.cs bit for bit; see that file and docs/design.md ("Affine forms") for each op's enclosure
// and the float32 rounding allowance it folds into e. Included by CultMath.hlsl after
// CultMath.Interval.hlsl, whose CULTMATH_SNOISE_LIPSCHITZ it shares, and after its own
// definition of cultmath_snoise_grad.
// An upper bound on the spectral norm of the Hessian of cultmath_snoise(float3); provenance in
// math.Affine.cs (SNOISE_HESSIAN).
static const float CULTMATH_SNOISE_HESSIAN = 56.050385;
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); }
float3 cultmath_af_from_iv(float2 interval)
{
float x0 = interval.x * 0.5 + interval.y * 0.5;
float e = interval.y * 0.5 - interval.x * 0.5;
return float3(x0, 0.0, e + (abs(x0) + e) * 9.5367431640625e-7 + 1.17549435e-38);
}
// An endpoint that would be NaN is the infinity on its side. NaN is found by its bit pattern: dxc marks
// compares fast, so a test that needs NaN to make `lo == lo` false is deleted under its default flags.
float2 cultmath_af_range(float3 x)
{
float r = abs(x.y) + x.z;
float w = abs(x.x) * 9.5367431640625e-7 + r * 9.5367431640625e-7;
float lo = x.x - r - w;
float hi = x.x + r + w;
return float2((asuint(lo) & 0x7fffffffu) > 0x7f800000u ? asfloat(0xff800000u) : lo,
(asuint(hi) & 0x7fffffffu) > 0x7f800000u ? asfloat(0x7f800000u) : hi);
}
float3 cultmath_af_add(float3 x, float3 y)
{
float x0 = x.x + y.x;
float a = x.y + y.y;
float e = x.z + y.z;
return float3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7);
}
float3 cultmath_af_sub(float3 x, float3 y)
{
float x0 = x.x - y.x;
float a = x.y - y.y;
float e = x.z + y.z;
return float3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7);
}
float3 cultmath_af_neg(float3 x) { return float3(-x.x, -x.y, x.z); }
float3 cultmath_af_scale(float3 x, float s)
{
float x0 = x.x * s;
float a = x.y * s;
float e = x.z * abs(s);
return float3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7 + 1.17549435e-38);
}
float3 cultmath_af_add_iv(float3 x, float2 interval) { return cultmath_af_add(x, cultmath_af_from_iv(interval)); }
float3 cultmath_af_mul(float3 x, float3 y)
{
float x0 = x.x * y.x;
float p = x.x * y.y;
float q = y.x * x.y;
float a = p + q;
float e = abs(x.x) * y.z + abs(y.x) * x.z + (abs(x.y) + x.z) * (abs(y.y) + y.z);
return float3(x0, a, e + (abs(x0) + abs(p) + abs(q) + e) * 9.5367431640625e-7 + 1.17549435e-38);
}
// The narrower of the Lipschitz form and the centred (Taylor) form over every p = centre + axis eps + w,
// |w| <= radius, from one cultmath_snoise_grad at the centre; the allowance for the float32 evaluation is
// in math.Affine.cs.
float3 cultmath_af_snoise(float3 centre, float3 axis, float radius)
{
float4 g = cultmath_snoise_grad(centre);
float3 gradient = float3(g.x, g.y, g.z);
float reach = length(axis) + radius;
float lipschitz = CULTMATH_SNOISE_LIPSCHITZ * reach;
float a = dot(gradient, axis);
float remainder = length(gradient) * radius + 0.5 * CULTMATH_SNOISE_HESSIAN * reach * reach;
float allowance = 6.103515625e-5 + 2.44140625e-4 * reach
+ CULTMATH_SNOISE_LIPSCHITZ * 9.5367431640625e-7 * (abs(centre.x) + abs(centre.y) + abs(centre.z) + reach);
bool centred = abs(a) + remainder < lipschitz;
float slope = centred ? a : 0.0;
float e = (centred ? remainder : lipschitz) + allowance;
return float3(g.w, slope, e + (abs(g.w) + abs(slope) + e) * 9.5367431640625e-7);
}
// Octave sum of cultmath_af_snoise, compounding frequency and amplitude as the value-and-gradient fBm does.
float3 cultmath_af_fbm(float3 centre, float3 axis, float radius, int octaves, float lacunarity, float gain)
{
octaves = clamp(octaves, 0, 16);
float amplitude = 1.0;
float frequency = 1.0;
float3 sum = float3(0.0, 0.0, 0.0);
for (int i = 0; i < octaves; i++)
{
sum = cultmath_af_add(sum, cultmath_af_scale(cultmath_af_snoise(centre * frequency, axis * frequency, radius * abs(frequency)), amplitude));
frequency *= lacunarity;
amplitude *= gain;
}
return sum;
}
// The axis of a screen tile's slice over the depth segment [z0, z1]; camera frame of
// cultmath_iv_frustum_ball.
float3 cultmath_af_frustum_axis(float2 centreSlope, float z0, float z1)
{
float h = (z1 - z0) * 0.5;
return float3(centreSlope.x * h, centreSlope.y * h, h);
}
// The affine slice of a screen tile over [z0, z1]: centre xyz, radius w. Derivation in math.Affine.cs.
float4 cultmath_af_frustum_ball(float2 centreSlope, float z0, float z1, float footprintPerDepth, float warpVariation)
{
float zm = (z0 + z1) * 0.5;
float h = (z1 - z0) * 0.5;
float cx = centreSlope.x * zm;
float cy = centreSlope.y * zm;
float spread = (abs(centreSlope.x) + abs(centreSlope.y) + 1.0) * h;
float radius = z1 * footprintPerDepth + warpVariation;
radius += (abs(cx) + abs(cy) + abs(zm) + spread + radius) * 9.5367431640625e-7;
return float4(cx, cy, zm, radius);
}
#endif

7
packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Affine.hlsl.meta Normal file
View file

@ -0,0 +1,7 @@
fileFormatVersion: 2
guid: 5dbcf0fdd66a4a74b22f5ca1c8e33058
ShaderIncludeImporter:
externalObjects: {}
userData:
assetBundleName:
assetBundleVariant:

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

@ -9,6 +9,9 @@
// Lipschitz constant of cultmath_snoise(float3); provenance in math.Interval.cs (SNOISE_LIPSCHITZ).
static const float CULTMATH_SNOISE_LIPSCHITZ = 7.9640074;
// Lipschitz constant of cultmath_snoise(float2); provenance in math.Interval.cs (SNOISE2_LIPSCHITZ).
static const float CULTMATH_SNOISE2_LIPSCHITZ = 8.117002;
float2 cultmath_iv_point(float value) { return float2(value, value); }
float2 cultmath_iv_add(float2 a, float2 b) { return float2(a.x + b.x, a.y + b.y); }
float2 cultmath_iv_sub(float2 a, float2 b) { return float2(a.x - b.y, a.y - b.x); }
@ -87,6 +90,14 @@ float2 cultmath_iv_snoise_ball(float3 centre, float radius)
return float2(max(n - e, -1.0), min(n + e, 1.0));
}
// The 2D overload: [n - L2 r, n + L2 r] intersected with [-1, 1], n = cultmath_snoise(centre). One snoise, value only.
float2 cultmath_iv_snoise_ball(float2 centre, float radius)
{
float n = cultmath_snoise(centre);
float e = CULTMATH_SNOISE2_LIPSCHITZ * radius;
return float2(max(n - e, -1.0), min(n + e, 1.0));
}
// Octave sum of cultmath_iv_snoise_ball, compounding frequency and amplitude as the value-and-gradient fBm does.
float2 cultmath_iv_fbm_ball(float3 centre, float radius, int octaves, float lacunarity, float gain)
{

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

@ -509,6 +509,7 @@ CultCellular cultmath_cellular(float3 p)
#include "CultMath.Phacelle.hlsl"
#include "CultMath.Interval.hlsl"
#include "CultMath.Affine.hlsl"
float cultmath_value_noise(float2 position)
{