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

Merge cellular, smooth-min gradient and the struct mirror into CultMath (Asura cut 2a-i)

Browse source 
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Meta Vens 2026-09-26 02:24:47 +02:00
commit 268e0ef17c
9 changed files with 1325 additions and 12 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -100,7 +100,10 @@ needs anything else, the test fails:
no CultMath analog.
3. File-scope `static const` becomes `const`.
4. Floating literals gain `f` (`0.5` becomes `0.5f`).
5. The file body is wrapped in a class, since C# has no free functions; HLSL
5. Struct fields gain `public`. HLSL struct members have no access-modifier
concept; a C# struct's fields default to `private`, which would hide them
from the mirror test's reflection-based field comparison.
6. The file body is wrapped in a class, since C# has no free functions; HLSL
functions become private instance methods.
HLSL constructs that stay outside the target and must not appear in the mirror:
@ -120,6 +123,116 @@ compiled under `using static CultMath.math;`, so its intrinsics are C# `math`
itself: the test proves the composition in the text, not the intrinsic rules,
which only `HlslSemanticsTests` pins.
The comparison also allows exactly one struct return shape: a mirror function
whose HLSL return type is a plain struct of `float`/`floatN`/`int`/`intN`
fields is compared to its C# counterpart field by field, recursing into each
field's own components, down to bit-for-bit scalars, matching each leaf's
declaration path (field names, not just position) and declared type as well
as its value. `CultCellular` (`math.cs`, `cellular`) is the first and, for
now, only consumer. A struct return with a mismatched, reordered-by-name, or
retyped field does not get a second comparison path or a shape of its own;
it fails the same walk that already handles a bare vector return.
## Invariant 8: Value-and-Gradient Primitives
A primitive tagged invariant 8 returns its value and analytic gradient
together, laid out as `float4(gradient.xyz, value.w)`, with no value-only twin
(a value-only form would be a second path with no consumer). `smin_grad` and
`cellular` are the first primitives in this family:
- `smin_grad(float4 a, float4 b, float k)` is Inigo Quilez's quadratic-polynomial
smooth minimum ("smooth minimum",
<https://iquilezles.org/articles/smin/>), carried to value-and-gradient
form. Its blend factor `h` is affine in `b.w - a.w`, so differentiating the
value with respect to position, the terms carrying `dh/dp` cancel exactly;
the surviving gradient is `lerp(∇b, ∇a, h)`, the same `h` the value uses.
This is the analytic gradient, not an approximation, and it is continuous
across the `|a.w - b.w| = k` seam where `h` saturates to 0 or 1.
- `cellular(float3 p)` is Worley's cellular texture basis function ("A
Cellular Texture Basis Function", SIGGRAPH 1996), searched over the
jittered 5×5×5 neighbourhood of feature points selected by `pcg3d` (never
the sin-based `hash`), with cells whose position-only lower bound already
exceeds the running F2 skipped before their hash is even computed. It
returns `CultCellular { nearest, edge, id }`: `nearest` is `(∇F1, F1)`,
`edge` is `(∇F2 - ∇F1, F2 - F1)`, and `id` is the nearest cell's `pcg4d`
hash (over the cell coordinate alone, never the `pcg3d` hash that produced
its jitter) mapped to `[0, 1)`. `∇F1` and `∇F2` are each undefined exactly
at their own feature point (`F1 = 0` or `F2 = 0`); `cellular` follows the
repo's existing degenerate-normal convention
(`GameCult.Geometry.CultGeometryIsoSurface.EmitOrientedTriangle`, which
guards a zero-length normal and returns the zero vector instead of the NaN
a bare `normalize` gives there) and returns the zero vector rather than NaN
for either. `F1 = 0` is ordinary (`p` sits on a feature point); `F2 = 0`
needs two distinct cells' `pcg3d`-jittered feature points to land on the
same float32 value in every component, which is unreachable at ordinary
magnitudes but confirmed reachable (see "Precision domain" below). The
identity of the nearest and second-nearest feature point changes
discontinuously across the `F1 = F2` set (a genuine kink, not a numerical
artifact), so `∇F1` and `∇F2` are only continuous away from that set.
The 5×5×5 radius is exact, not merely safer than the classic 3×3×3 radius:
for a query at position `u` inside its own cell (`u` in `[0, 1)^3` per
component) and a candidate feature at integer cell offset `d` with jitter
`j` in `[0, 1)^3`, the per-axis displacement `(u_i - j_i) - d_i` is always
strictly greater in magnitude than `max(0, |d_i| - 1)`, because
`u_i - j_i` is always strictly inside `(-1, 1)`. So any cell with `|d_i| >=
3` on some axis is strictly farther than 2 from the query. Two real,
always-present candidates bound `F1` and `F2`: the query's own cell, whose
distance is always strictly less than `sqrt(3)` (every axis strictly under
1); and the near-side neighbour on whichever single axis has `u_i` closest
to a cell boundary. Let `m = min(u_i, 1 - u_i)` on that chosen axis, so
`m` is in `[0, 1/2]` by construction. Stepping one cell towards that
nearer boundary puts the query within `1 + m` of the far face on the
offset axis and within `1 - m` of the near face on each of the other two
axes, so that neighbour's squared distance is at most
`(1 + m)^2 + 2*(1 - m)^2 = 3 - 2m + 3m^2 <= 3` for every `m` in `[0, 1/2]`
(3 at `m = 0`, dipping to a minimum of `8/3 ~= 2.667` at `m = 1/3`, then
rising back to 2.75 at `m = 1/2`), and strictly below 3 because
the jitter never reaches exactly 0 or 1. `F2` is at most the larger of
these two real candidates (any two real candidates upper-bound the
2nd-smallest value over the full infinite candidate set), so `F1 <
sqrt(3)` and `F2 < sqrt(3)` always, comfortably under the radius-2
exclusion distance of 2, with margin `2 - sqrt(3) ~= 0.268` (Soul's
adversarial check, cut 2a-i: the worst realizable in-box `F2` is exactly
`sqrt(3)`, reached only in the unattainable limit at a cell corner). 3×3×3's
matching argument only reaches radius 1, whose excluded cells start at
distance 1, short of the `sqrt(3)` ceiling — the gap its measured failures
(9 of 400,000 points in `[-50, 50]^3`, worst error 0.088) live in.
An earlier version of this proof picked the near-side neighbour on a fixed
axis instead of whichever axis sits closest to its own boundary, and that
version is wrong (Soul, cut 2a-i): `u = (0.5, 0, 0)` with the x-neighbour's
jittered feature at `(-1, 1 - 2^-24, 1 - 2^-24)` gives a distance of about
2.06, outside the radius-2 exclusion the proof needs — `u.x = 0.5` is the
point on the cell farthest from any boundary, so `x` was the wrong axis to
fix there. Choosing the axis by proximity to its own boundary is what keeps
`m <= 1/2` and the bound at `sqrt(3)`.
`CellularAndSminGradTests` pins both the prune's bit-identical equivalence
to an unpruned 5×5×5 reference and the search's exact agreement with an
independent 7×7×7 brute-force oracle.
Precision domain: `cellular`'s exactness above is a statement about the
search radius, not about float32. Past `|p|` of roughly `2^23`, adjacent
integer cells stop being distinguishable from their neighbours at float32
precision in some directions, and distinct cells' jittered features can
round onto the same value, making `F1 = F2` ties (and, per-component,
`F2 = 0`) increasingly common from about `1e6` and confirmed at `2^24`
(measured: 4096 of 4096 sampled integer points at `2^24` hit this before
the guard above). `cellular` stays finite there; it does not stay accurate
in the sense of matching the idealized real-valued Worley field.
The lower-bound prune's bit-identity to the unpruned 5×5×5 search also has
a domain: it holds for `|p| < 2^25`. Above that, `cell + d` (an integer
cell coordinate plus a small integer offset) starts to round in float32,
so a pruned cell's own recomputed offset can land on a different neighbour
than the unpruned loop visits at that same nominal `(dx, dy, dz)`, and the
two search paths diverge (Soul, cut 2a-i: measured 180 of 1,100,000 points
differ at `|p| ~= 3.4e7`). This is a float32-precision limit on the prune's
equivalence proof, not a bug in the prune itself: below `2^25` the two
paths are proven bit-identical (`PrunedSearchIsBitIdenticalToTheUnprunedFiveCubedSearch`,
which samples `[-50, 50]^3`, is comfortably inside that domain).
Integer hashing uses the PCG hashes from Jarzynski and Olano, "Hash Functions
for GPU Rendering" (JCGT 9(3), 2020): `pcg(uint)` is O'Neill's RXS-M-XS 32/32
permutation over one LCG step, and `pcg3d`/`pcg4d` are the paper's (3 → 3) and
@ -128,9 +241,28 @@ uint multiply and add modulo 2^32 exactly as C# unchecked arithmetic does, so
they are integer-exact in HLSL and on the CPU, GPU included, unlike the
`sin`-based `hash`. CultMath has no uint vectors, so `pcg3d(int3)` and
`pcg4d(int4)` carry uint bit patterns in int components. The `float2`/`float3`/`float4`
overloads hash the IEEE-754 bits (`asuint`), so -0 and 0 hash differently.
`pcg3d(float2)` holds z at 0, which is the paper's route for (2 → N) inputs.
Seeds and other scalar uses take one output component.
overloads hash the IEEE-754 bits (`asuint`), so -0 and 0 hash differently; they
exist for callers that already hold a hash key as float bits, and `pcg3d(float2)`
holds z at 0, which is the paper's route for (2 → N) inputs. Seeds and other
scalar uses take one output component.
`cellular`'s jitter and id always hash the integer cell coordinate directly
(`pcg3d(int3(neighbor))`, `pcg4d(int4(int3(cellId), 0))`), never a cell
coordinate's float32 bit pattern. Hashing the bit pattern of an
integer-valued float is not the same permutation as hashing the integer: two
adjacent cells' bit patterns share more structure (both have the same
exponent and a mantissa that differs by one increment) than two adjacent
integers do going into `pcg3d`'s own mixing, and that residual structure
survived into the jitter. Hashing the integer cell coordinate directly is
exact and runtime-independent, and it is the HLSL mirror's own hash, so C#
and shader agree bit-for-bit. What pins this is the oracle and mirror tests,
not a uniformity measurement: reverting to hashing the float bits (or the
sin-based hash) is caught by `ProductionSearchMatchesTheBruteForceOracleExactly`,
`EveryMirrorFunctionMatchesCSharpMath`, `F1NeverExceedsF2`, and the other
tests that key off the integer-hash fixtures (cut 2a-i). `CellularAndSminGradTests`
(`JitterAndIdPassAChiSquareUniformityTest`) checks the uniformity of the hash
as `cellular` uses it: the jitter marginals and the id pass a chi-square test
at a stated significance against `pcg3d(int3(...))`.
Where HLSL is silent, CultMath keeps its own decisions and does not defer to
Unity.Mathematics: `hash` returns float, and `Random` is CultMath's own xorshift32. Engine-shaped

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

@ -313,6 +313,88 @@ int4 cultmath_pcg4d(int4 value)
int4 cultmath_pcg4d(float4 value) { return cultmath_pcg4d(int4((int)asuint(value.x), (int)asuint(value.y), (int)asuint(value.z), (int)asuint(value.w))); }
// Inigo Quilez, "smooth minimum": quadratic-polynomial smin carried to value-and-gradient form. See
// math.smin_grad's comment (math.cs) for why the gradient blend is exact, not an approximation.
float4 cultmath_smin_grad(float4 a, float4 b, float k)
{
float h = saturate(0.5 + 0.5 * (b.w - a.w) / k);
float value = lerp(b.w, a.w, h) - k * h * (1.0 - h);
float3 gradient = lerp(b.xyz, a.xyz, h);
return float4(gradient, value);
}
// Invariant 8's one struct return shape (design.md, "HLSL Target: Source Transformations"):
// compared field by field and bit for bit against CultMath.CultCellular, same field order.
struct CultCellular
{
float4 nearest;
float4 edge;
float id;
};
// Worley, "A Cellular Texture Basis Function" (SIGGRAPH 1996). See math.cellular's comment
// (math.cs) for the hashing, gradient, search-radius exactness proof, id source and
// degenerate-point conventions this mirrors exactly, including the 5x5x5 lower-bound prune.
CultCellular cultmath_cellular(float3 p)
{
float3 cell = floor(p);
float f1 = 1.0e30;
float f2 = 1.0e30;
float3 c1 = float3(0.0, 0.0, 0.0);
float3 c2 = float3(0.0, 0.0, 0.0);
float3 cellId = float3(0.0, 0.0, 0.0);
for (int dz = -2; dz <= 2; dz++)
{
float lz = max(0.0, abs((float)dz) - 1.0);
for (int dy = -2; dy <= 2; dy++)
{
float ly = max(0.0, abs((float)dy) - 1.0);
for (int dx = -2; dx <= 2; dx++)
{
float lx = max(0.0, abs((float)dx) - 1.0);
float lowerBound = sqrt(lx * lx + ly * ly + lz * lz);
if (lowerBound >= f2)
continue;
float3 neighbor = cell + float3(dx, dy, dz);
int3 hash = cultmath_pcg3d(int3(neighbor));
float3 jitter = float3(
((uint)hash.x >> 8) * (1.0 / 16777216.0),
((uint)hash.y >> 8) * (1.0 / 16777216.0),
((uint)hash.z >> 8) * (1.0 / 16777216.0));
float3 feature = neighbor + jitter;
float d = cultmath_distance(p, feature);
if (d < f1)
{
f2 = f1; c2 = c1;
f1 = d; c1 = feature; cellId = neighbor;
}
else if (d < f2)
{
f2 = d; c2 = feature;
}
}
}
}
// F1 = 0 and F2 = 0 share the same degenerate-gradient guard; see math.cellular's comment
// (math.cs) for why F2 = 0 is reachable past roughly |p| = 2^23.
float3 grad1 = f1 > 0.0 ? (p - c1) / f1 : float3(0.0, 0.0, 0.0);
float3 grad2 = f2 > 0.0 ? (p - c2) / f2 : float3(0.0, 0.0, 0.0);
// id comes from pcg4d over the winning cell's integer coordinate alone (same int-hashing rule
// as the jitter above), never from the pcg3d hash that produced its jitter (F4).
int4 idHash = cultmath_pcg4d(int4(int3(cellId), 0));
CultCellular result;
result.nearest = float4(grad1, f1);
result.edge = float4(grad2 - grad1, f2 - f1);
result.id = ((uint)idHash.w >> 8) * (1.0 / 16777216.0);
return result;
}
float cultmath_value_noise(float2 position)
{
float2 cell = floor(position);

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

@ -0,0 +1,26 @@
namespace CultMath;
/// <summary>
/// The value-and-gradient result of <see cref="math.cellular(float3)"/> (invariant 8's one struct
/// return shape). Every field follows the same <c>float4(gradient.xyz, value.w)</c> layout as a
/// bare value-and-gradient return.
/// </summary>
public struct CultCellular
{
/// <summary>(&#8711;F1, F1): the gradient and distance to the nearest feature point.</summary>
public float4 nearest;
/// <summary>(&#8711;F2 - &#8711;F1, F2 - F1): for lineae along cell borders.</summary>
public float4 edge;
/// <summary>The nearest cell's pcg4d hash, mapped to [0, 1); independent of the pcg3d hash that
/// produced that cell's jitter.</summary>
public float id;
public CultCellular(float4 nearest, float4 edge, float id)
{
this.nearest = nearest;
this.edge = edge;
this.id = id;
}
}

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

@ -0,0 +1,2 @@
fileFormatVersion: 2
guid: 76be9498b63e410891715b3ec965caf8

151
packages/cultmath/src/CultMath/math.cs
View file

@ -1,5 +1,10 @@
using System;
// Cut 2a-i fix batch 1: cellular_unit's [0,1) mapping bug (F7) is only reachable from cellular's
// private call sites; this lets CultMath.Tests probe it directly instead of hunting for a hash
// input that happens to land on the bit pattern under test.
[assembly: System.Runtime.CompilerServices.InternalsVisibleTo("CultMath.Tests")]
namespace CultMath;
public static partial class math
@ -688,6 +693,152 @@ public static partial class math
private static float3 snoise_permute(float3 value) => snoise_mod289(((value * 34.0f) + 1.0f) * value);
private static float4 snoise_permute(float4 value) => snoise_mod289(((value * 34.0f) + 1.0f) * value);
// Inigo Quilez, "smooth minimum" (https://iquilezles.org/articles/smin/): the quadratic-polynomial
// smin, carried to value-and-gradient form (invariant 8). h is affine in (b.w - a.w) inside the
// smoothing band, so differentiating value = lerp(b.w, a.w, h) - k*h*(1-h) with respect to
// position, the terms carrying dh/dp cancel exactly (dh/dp * [(a.w-b.w) - k + 2*k*h] and the
// bracket is identically zero given how h is built from a.w-b.w). What is left is exactly
// lerp(∇b, ∇a, h): the gradient blend is not an approximation, it is the analytic gradient.
// Outside the band (|a.w - b.w| >= k), h saturates to 0 or 1 and this is exactly min(a, b) with
// that input's own gradient, continuously (h and the correction term both reach the boundary at
// the same value from the smooth side, so there is no kink to exclude here the way cellular's
// F1 = F2 seam needs one).
public static float4 smin_grad(float4 a, float4 b, float k)
{
var h = saturate(0.5f + 0.5f * (b.w - a.w) / k);
var value = lerp(b.w, a.w, h) - k * h * (1.0f - h);
var gradient = lerp(new float3(b.x, b.y, b.z), new float3(a.x, a.y, a.z), h);
return new float4(gradient, value);
}
// Worley, "A Cellular Texture Basis Function" (SIGGRAPH 1996): F1/F2 and their analytic
// gradients, searched over the jittered neighbourhood of feature points. Feature points are
// selected with pcg3d over each neighbour's integer cell coordinate (design.md, "Integer
// hashing"), never the float32 bit pattern of that coordinate and never the sin-based hash.
// Hashing the integer cell coordinate directly is exact and runtime-independent, and it is what
// the HLSL mirror does, so C# and shader stay bit-for-bit aligned. What pins this is the oracle
// and mirror tests, not a uniformity measurement: reverting to hashing the float32 bit pattern
// (or the sin-based hash) is caught by ProductionSearchMatchesTheBruteForceOracleExactly,
// EveryMirrorFunctionMatchesCSharpMath, F1NeverExceedsF2, and the other tests that key off the
// integer-hash fixtures. The
// gradient of a distance field, del|p - c|, is the unit vector (p - c)/|p - c|; that is undefined exactly at
// a feature point (F1 = 0 or F2 = 0), so cellular follows the repo's existing degenerate-normal
// convention (GameCult.Geometry.CultGeometryIsoSurface.EmitOrientedTriangle: guard the zero-length
// case and return the zero vector instead of the NaN a bare normalize would produce there).
//
// The search radius is 2 cells on every axis (5x5x5), not the classic 1-cell (3x3x3) radius,
// because 3x3x3 is not exact: measured 9 of 400,000 points in [-50,50]^3 land on the wrong F2
// (worst error 0.088), and past |p| of about 2^23 the float32 jitter can put two distinct cells'
// features on the same value, which 3x3x3 has no chance of noticing either way. 5x5x5 is exact,
// not merely safer, by the following argument:
//
// Let u = p - cell be the query's position inside its own cell (u in [0,1)^3 per component), and
// let a candidate feature at integer cell offset d = (dx, dy, dz) sit at d + j for some jitter j
// in [0,1)^3. Per axis, the displacement from p to that feature is (u_i - j_i) - d_i; since
// u_i - j_i is always strictly inside (-1, 1) regardless of the actual u_i and j_i values,
// |displacement_i| is always strictly greater than max(0, |d_i| - 1). So every cell with |d_i| >=
// 3 on some axis is strictly farther than 2 from p (that axis alone already gives max(0, 3-1) =
// 2). Now bound F1 and F2 from above using two real, always-present candidates: the query's own
// cell (d = 0), whose distance is always strictly less than sqrt(3) (every axis strictly under
// 1); and the near-side neighbour on whichever single axis's u_i sits closest to a cell
// boundary. Let m = min(u_i, 1 - u_i) on that chosen axis (m is in [0, 1/2] by construction,
// since it is the smaller of the two). Stepping one cell towards that nearer boundary puts the
// query within (1 + m) of the far face on the offset axis and within (1 - m) of the near face on
// each of the other two axes, so that neighbour's squared distance is at most
// (1 + m)^2 + 2*(1 - m)^2 = 3 - 2m + 3m^2, which is <= 3 for every m in [0, 1/2] (3 at m = 0,
// dipping to a minimum of 8/3 ~= 2.667 at m = 1/3, then rising back to 2.75 at m = 1/2), and
// strictly below 3 because the jitter never reaches exactly 0 or
// 1 (u_i - j_i stays strictly inside (-1, 1), so this bound, like the own-cell one, is never
// tight). F2 is at most the larger of these two real candidates (adding candidates to a set can
// only lower or hold its 2nd-smallest value), so F1 < sqrt(3) and F2 < sqrt(3) always: roughly
// 1.732 < 2, margin 2 - sqrt(3) ~= 0.268 (Soul's adversarial check, cut 2a-i: the worst
// realizable in-box F2 is exactly sqrt(3), reached only in the unattainable limit at a cell
// corner). Every cell outside the 5x5x5 block is farther than 2, hence farther than either F1 or
// F2, so the true global F1 and F2 are always inside it. (The matching argument for 3x3x3 only
// reaches radius 1, whose excluded cells start at distance > 1, which is less than the 1.732
// ceiling — the gap 3x3x3's failures live in.)
//
// An earlier version of this proof picked the near-side neighbour on a fixed axis rather than
// whichever axis sits closest to its boundary, and that version is wrong (Soul, cut 2a-i):
// u = (0.5, 0, 0) with the x-neighbour's jittered feature at (-1, 1 - 2^-24, 1 - 2^-24) gives a
// distance of about 2.06, outside the radius-2 exclusion the proof needs. Choosing the axis by
// proximity to its own boundary is what keeps m <= 1/2 and the bound at sqrt(3); x was the wrong
// axis to fix there because u.x = 0.5 is the point on the whole cube farthest from any boundary.
//
// The loop below still visits every one of the 125 cells in a fixed order, but skips the
// pcg3d/jitter/distance work for a cell whose distance cannot possibly beat the current f2: that
// per-cell lower bound (the max(0, |d_i| - 1) term above, combined across axes) depends only on
// the integer offset, never on p or the hash, so pruning a cell the unpruned loop would also
// have rejected changes no output bit (CellularExactSearchTests pins the pruned/unpruned
// equivalence and the exactness against a 7x7x7 brute-force oracle).
public static CultCellular cellular(float3 p)
{
var cell = floor(p);
var f1 = 1.0e30f;
var f2 = 1.0e30f;
var c1 = new float3(0.0f, 0.0f, 0.0f);
var c2 = new float3(0.0f, 0.0f, 0.0f);
var cellId = new float3(0.0f, 0.0f, 0.0f);
for (var dz = -2; dz <= 2; dz++)
{
var lz = max(0.0f, abs((float)dz) - 1.0f);
for (var dy = -2; dy <= 2; dy++)
{
var ly = max(0.0f, abs((float)dy) - 1.0f);
for (var dx = -2; dx <= 2; dx++)
{
var lx = max(0.0f, abs((float)dx) - 1.0f);
var lowerBound = sqrt(lx * lx + ly * ly + lz * lz);
if (lowerBound >= f2)
continue;
var neighbor = cell + new float3(dx, dy, dz);
var hash = pcg3d(int3(neighbor));
var jitter = new float3(cellular_unit(hash.x), cellular_unit(hash.y), cellular_unit(hash.z));
var feature = neighbor + jitter;
var d = distance(p, feature);
if (d < f1)
{
f2 = f1; c2 = c1;
f1 = d; c1 = feature; cellId = neighbor;
}
else if (d < f2)
{
f2 = d; c2 = feature;
}
}
}
}
// F1 = 0 is reachable (p sits exactly on a feature point). F2 = 0 needs two distinct cells'
// pcg3d-jittered features to land on the same float32 value in every component: unreachable
// at ordinary magnitudes, but confirmed reachable past roughly |p| = 2^23 (design.md,
// "Precision domain"), so both share the same degenerate-gradient guard.
var grad1 = f1 > 0.0f ? (p - c1) / f1 : new float3(0.0f, 0.0f, 0.0f);
var grad2 = f2 > 0.0f ? (p - c2) / f2 : new float3(0.0f, 0.0f, 0.0f);
// id comes from pcg4d over the winning cell's integer coordinate alone (same int-hashing
// rule as the jitter above), never from the pcg3d hash that produced its jitter: sharing a
// hash between id and jitter would make id a deterministic function of the jitter it is
// supposed to be independent from, correlated across every cell (F4).
var idHash = pcg4d(new int4(int3(cellId), 0));
return new CultCellular(
new float4(grad1, f1),
new float4(grad2 - grad1, f2 - f1),
cellular_unit(idHash.w));
}
// Maps a hash component's uint bit pattern to [0, 1) using the top 24 bits, so the result is
// always strictly less than 1: a plain divide-by-2^32 rounds 0xFFFFFF80..0xFFFFFFFF up to
// exactly 1.0f in float32, outside the documented range. (uint >> 8) has at most 24 significant
// bits, which float32's 24-bit mantissa represents exactly, and the division is by an exact power
// of two, so both the shift-and-convert and the divide are the same IEEE-754 round-to-nearest on
// dxc and C#, bit-exact on both sides.
internal static float cellular_unit(int bits) => ((uint)bits >> 8) * (1.0f / 16777216.0f);
public static float value_noise(float2 position)
{
var cell = floor(position);

686
packages/cultmath/tests/CultMath.Tests/CellularAndSminGradTests.cs Normal file
View file

@ -0,0 +1,686 @@
using System;
using CultMath;
using Xunit;
namespace CultMath.Tests;
public sealed class CellularAndSminGradTests
{
// Central-difference tolerances for smin_grad and cellular (F5), sized separately: the two
// functions have different curvature, so the same eps does not carry the same truncation
// error, and an earlier shared 3e-4 tolerance was fragile against both (Soul, cut 2a-i: measured
// max error smin 2.74e-4, cellular nearest 1.25e-4, cellular edge 1.88e-4 on a small sample —
// and re-measured here at 5.3e-4, 2.0e-4, 2.6e-4 on a larger one — against that 3e-4 tolerance,
// sometimes under 1.1x headroom, not the "~6e-5 noise" an earlier version of this comment
// claimed). Each tolerance below is set from its own measured worst case with real headroom.
// smin_grad's value is exactly quadratic in position within a smoothing band (math.cs comment on
// smin_grad: h is affine in b.w - a.w, which is itself affine in position, so value(q) has zero
// third derivative there): central differences have no truncation error from curvature, only
// float32 rounding noise, so a larger eps directly dilutes that noise with no accuracy cost.
// Measured max error at eps=0.01f over 10,000 samples: 4.37e-5 (up from eps=0.001f's 5.3e-4,
// confirming the noise is eps-limited, not curvature-limited); this tolerance keeps about 5.7x
// headroom over that. A uniform-scale mutant on the returned gradient (a 1.01x M9/M16-style
// mutant applied here) still dies down to a 1.0005x scale (0.05%); it survives at 1.0002x.
private const float SminGradientTolerance = 2.5e-4f;
// cellular's F1/F2 are Euclidean distance fields: real curvature away from a feature point means
// central differences DO carry O(eps^2) truncation error, so (unlike smin_grad) a larger eps
// makes things worse, not better; a sweep at eps in {0.0002f, 0.0003f, 0.0005f, 0.0007f, 0.001f}
// over 8,000 samples found the committed eps=0.0005f already close to the sweet spot. Measured
// max error at that eps: nearest 2.04e-4, edge 2.57e-4; this tolerance keeps just over 5x
// headroom over the larger of the two. M9 (grad1 * 1.01) and M16 (grad2 * 1.01) both still die
// down to a 1.002x scale (0.2%); both survive at 1.001x.
private const float CellularGradientTolerance = 1.3e-3f;
// ---- smin_grad ----
[Fact]
public void SminValueNeverExceedsTheOrdinaryMinimum()
{
var random = new System.Random(0xACE1);
for (var i = 0; i < 256; i++)
{
var a = new float4(RandomUnit(random), RandomScalar(random));
var b = new float4(RandomUnit(random), RandomScalar(random));
var k = 0.05f + random.NextSingle() * 5.0f;
var result = math.smin_grad(a, b, k);
Assert.True(result.w <= math.min(a.w, b.w) + 1.0e-4f, $"a.w={a.w} b.w={b.w} k={k} smin={result.w}");
}
}
[Fact]
public void SminMatchesTheHardMinimumAndItsGradientOutsideTheBand()
{
var random = new System.Random(0xACE2);
for (var i = 0; i < 64; i++)
{
var gradA = RandomUnit(random);
var gradB = RandomUnit(random);
var k = 0.1f + random.NextSingle() * 2.0f;
// Push b comfortably above a, past the smoothing band, so a is the ordinary minimum.
var aw = random.NextSingle() * 10.0f - 5.0f;
var bw = aw + k * (1.5f + random.NextSingle());
var a = new float4(gradA, aw);
var b = new float4(gradB, bw);
var lower = math.smin_grad(a, b, k);
Assert.Equal(aw, lower.w, precision: 4);
Assert.Equal(gradA.x, lower.x, precision: 4);
Assert.Equal(gradA.y, lower.y, precision: 4);
Assert.Equal(gradA.z, lower.z, precision: 4);
// Argument order must not matter: b below a is symmetric.
var higher = math.smin_grad(b, a, k);
Assert.Equal(aw, higher.w, precision: 4);
Assert.Equal(gradA.x, higher.x, precision: 4);
Assert.Equal(gradA.y, higher.y, precision: 4);
Assert.Equal(gradA.z, higher.z, precision: 4);
}
}
[Fact]
public void SminGradientMatchesCentralDifferencesOfItsOwnValue()
{
var random = new System.Random(0xACE3);
var tested = 0;
while (tested < 200)
{
var cf = RandomUnit(random);
var cg = RandomUnit(random);
var offF = random.NextSingle() * 4.0f - 2.0f;
var offG = random.NextSingle() * 4.0f - 2.0f;
var k = 0.2f + random.NextSingle() * 2.0f;
var p = new float3(random.NextSingle() * 4.0f - 2.0f, random.NextSingle() * 4.0f - 2.0f, random.NextSingle() * 4.0f - 2.0f);
float F(float3 q) => math.dot(q, cf) + offF;
float G(float3 q) => math.dot(q, cg) + offG;
float Smin(float3 q) => math.smin_grad(new float4(cf, F(q)), new float4(cg, G(q)), k).w;
// Stay clear of the |a.w - b.w| = k seam: central differences straddling a kink would
// disagree with either side's analytic gradient (the same reason cellular excludes a
// band around its own F1 = F2 set). The band has an absolute floor, not just a
// fraction of k: eps grew to 0.01f below (SminGradientTolerance's comment), and a
// fraction-of-k-only band could be narrower than the diff can shift by eps*|cf-cg|
// (up to about 0.035 here) when k is small, letting p +/- eps step across the seam even
// though p itself was clear of it.
var diff = F(p) - G(p);
if (Math.Abs(Math.Abs(diff) - k) < Math.Max(0.05f * k, 0.05f))
continue;
tested++;
var analytic = math.smin_grad(new float4(cf, F(p)), new float4(cg, G(p)), k);
// eps=0.01f, not 0.001f: smin is exactly quadratic within a band (no truncation error
// from curvature; see SminGradientTolerance's comment), so the larger eps only dilutes
// float32 rounding noise, at no accuracy cost, while still comfortably inside the band
// exclusion above.
const float eps = 0.01f;
var numeric = new float3(
(Smin(p + new float3(eps, 0.0f, 0.0f)) - Smin(p - new float3(eps, 0.0f, 0.0f))) / (2.0f * eps),
(Smin(p + new float3(0.0f, eps, 0.0f)) - Smin(p - new float3(0.0f, eps, 0.0f))) / (2.0f * eps),
(Smin(p + new float3(0.0f, 0.0f, eps)) - Smin(p - new float3(0.0f, 0.0f, eps))) / (2.0f * eps));
AssertWithinTolerance(numeric.x, analytic.x, SminGradientTolerance);
AssertWithinTolerance(numeric.y, analytic.y, SminGradientTolerance);
AssertWithinTolerance(numeric.z, analytic.z, SminGradientTolerance);
}
}
[Fact]
public void SminGradientIsNotAConstantBlendOrTheSwappedFactor()
{
// Falsifies two named mutants directly: blending the gradient with a constant 0.5 instead of
// h, or with (1 - h) instead of h. Here a is comfortably the minimum (h saturates to 1 in the
// correct formula), so the correct gradient is exactly a's, unmixed. A constant-0.5 blend or a
// swapped (1 - h) blend both pull in b's gradient, which differs from a's here.
var a = new float4(1.0f, 0.0f, 0.0f, -3.0f);
var b = new float4(0.0f, 1.0f, 0.0f, 5.0f);
var result = math.smin_grad(a, b, 0.25f);
Assert.Equal(1.0f, result.x, precision: 5);
Assert.Equal(0.0f, result.y, precision: 5);
Assert.Equal(0.0f, result.z, precision: 5);
Assert.Equal(-3.0f, result.w, precision: 5);
}
// ---- cellular ----
[Fact]
public void F1NeverExceedsF2()
{
var random = new System.Random(0xCE11);
for (var i = 0; i < 256; i++)
{
var p = new float3(random.NextSingle() * 20.0f - 10.0f, random.NextSingle() * 20.0f - 10.0f, random.NextSingle() * 20.0f - 10.0f);
var result = math.cellular(p);
Assert.True(result.edge.w >= -1.0e-4f, $"F2 - F1 = {result.edge.w} at {p.x},{p.y},{p.z}");
// Both returned .w values, checked against the independent 7x7x7 brute-force oracle (F1,
// F6), not a value the test derived some other way.
var (f1, f2) = FindF1F2(p);
Assert.Equal(f1, result.nearest.w, precision: 3);
Assert.Equal(f2 - f1, result.edge.w, precision: 3);
}
}
[Fact]
public void NearestFeatureWinsEvenWhenItIsInTheLastVisitedNeighbourCell()
{
// math.cellular's own search order visits neighbour offsets with dz outermost and dx
// innermost over a 5x5x5 block (F6), and it also prunes cells whose position-only lower
// bound already exceeds the running F2. Offsets with any |component| = 2 are not provably
// unreachable as F1: their minimum possible distance is 1 (e.g. (2, 0, 0), whose lower bound
// per the exactness proof in math.cs is max(0, 2-1) = 1), well inside the query's own cell's
// sqrt(3) ceiling, so a |component| = 2 offset winning F1 is rare rather than impossible
// (Soul, cut 2a-i: measured 3 F1 wins in 2,000,000 samples). It is still rare enough that the
// 3x3x3 corners remain the practically-reachable "last" cases for this test. (dx, dy, dz)
// = (1, 1, 1) is the last of those the fixed iteration order visits, and if pruning or search
// order were broken, this is the offset most likely to be skipped. Rather than hand-picking a
// hash value, search random points for one whose true nearest feature (by independent brute
// force over the same 5x5x5 block) happens to live in exactly that cell, at a comfortable
// margin from the F1 = F2 seam, then check math.cellular agrees.
var random = new System.Random(0xCE14);
var lastOffset = new float3(1.0f, 1.0f, 1.0f);
for (var attempt = 0; attempt < 5000; attempt++)
{
var p = new float3(random.NextSingle() * 8.0f - 4.0f, random.NextSingle() * 8.0f - 4.0f, random.NextSingle() * 8.0f - 4.0f);
var (offset, f1, f2) = FindNearestOffset(p);
if (!offset.Equals(lastOffset) || f1 < 0.05f || f2 - f1 < 0.05f)
continue;
var cell = new float3(MathF.Floor(p.x), MathF.Floor(p.y), MathF.Floor(p.z));
var feature = FeaturePoint(cell + lastOffset);
var expectedGradient = (p - feature) / f1;
var result = math.cellular(p);
Assert.Equal(f1, result.nearest.w, precision: 4);
Assert.Equal(expectedGradient.x, result.nearest.x, precision: 4);
Assert.Equal(expectedGradient.y, result.nearest.y, precision: 4);
Assert.Equal(expectedGradient.z, result.nearest.z, precision: 4);
return;
}
Assert.Fail("did not find a random point whose nearest feature lives in the last-visited neighbour cell");
}
[Fact]
public void NearestGradientIsFiniteZeroAtAFeaturePointInsteadOfNaN()
{
var feature = FeaturePoint(new float3(4.0f, -2.0f, 9.0f));
var result = math.cellular(feature);
Assert.Equal(0.0f, result.nearest.w, precision: 4);
Assert.True(float.IsFinite(result.nearest.x) && float.IsFinite(result.nearest.y) && float.IsFinite(result.nearest.z));
Assert.Equal(0.0f, result.nearest.x, precision: 5);
Assert.Equal(0.0f, result.nearest.y, precision: 5);
Assert.Equal(0.0f, result.nearest.z, precision: 5);
}
[Fact]
public void GradientMatchesCentralDifferencesAwayFromTheF1EqualsF2Seam()
{
var random = new System.Random(0xCE12);
var tested = 0;
var attempts = 0;
while (tested < 120 && attempts < 20000)
{
attempts++;
var p = new float3(random.NextSingle() * 6.0f - 3.0f, random.NextSingle() * 6.0f - 3.0f, random.NextSingle() * 6.0f - 3.0f);
var (f1, f2, f3) = FindF1F2F3(p);
// Exclude a band around the F1 = F2 set: the identity of the nearest and second-nearest
// feature swap there, a genuine kink in F1 and F2 alike, not something central differences
// can approximate. Also exclude a band around F2 = F3: edge (F2 - F1) kinks there too,
// since the second-nearest feature's identity swaps (Soul, cut 2a-i: an unfiltered sample
// hit a 0.78 edge error this way). Also skip the exact feature point itself (F1 = 0),
// covered separately.
if (f2 - f1 < 0.05f || f1 < 0.01f || f3 - f2 < 0.05f)
continue;
tested++;
var result = math.cellular(p);
// eps=0.0005f: unlike smin_grad, cellular's F1/F2 are Euclidean distance fields with real
// curvature, so central differences carry O(eps^2) truncation error and a larger eps is
// worse, not better (see CellularGradientTolerance's comment for the measured sweep).
const float eps = 0.0005f;
// Differentiate cellular's own returned .w values (F1): calling cellular again at p +/-
// eps and reading result.nearest.w / result.edge.w, not a duplicate FindF1F2 search. A
// mutant that corrupts only the returned gradient, leaving the returned value alone,
// would survive a check built on an independent value oracle; it cannot survive this one.
float NearestValueAt(float3 q) => math.cellular(q).nearest.w;
float EdgeValueAt(float3 q) => math.cellular(q).edge.w;
AssertGradientMatchesCentralDifference(NearestValueAt, p, eps, result.nearest, CellularGradientTolerance);
AssertGradientMatchesCentralDifference(EdgeValueAt, p, eps, result.edge, CellularGradientTolerance);
}
Assert.True(tested >= 100, $"only {tested} of {attempts} random points landed away from the F1=F2 band");
}
[Fact]
public void IdIsConstantAcrossAFeaturesWholeRegionAndVariesBetweenCells()
{
var random = new System.Random(0xCE13);
var baseCell = new float3(2.0f, 2.0f, 2.0f);
var feature = FeaturePoint(baseCell);
var idAtFeature = math.cellular(feature).id;
for (var i = 0; i < 32; i++)
{
// Small offsets that stay inside the same Voronoi region (well short of any neighbour).
var jitter = new float3(random.NextSingle() * 0.2f - 0.1f, random.NextSingle() * 0.2f - 0.1f, random.NextSingle() * 0.2f - 0.1f);
var nearby = feature + jitter;
if (!FindNearestCell(nearby).Equals(baseCell))
continue;
Assert.Equal(idAtFeature, math.cellular(nearby).id, precision: 6);
}
var otherCellId = math.cellular(FeaturePoint(baseCell + new float3(5.0f, 0.0f, 0.0f))).id;
Assert.NotEqual(idAtFeature, otherCellId);
}
[Fact]
public void IdIsNotCorrelatedWithAnyJitterComponent()
{
// F4: id must come from hash bits the jitter never uses. A mutant that takes id from
// hash.x (or hash.y) makes id literally equal to that jitter component, which this catches
// as a correlation near 1.0, far past the 0.05 threshold. n = 20000 keeps sampling noise
// (measured under 0.01 at this n) well clear of that threshold.
const int n = 20000;
var random = new System.Random(0xCE15);
var ids = new double[n];
var jitterX = new double[n];
var jitterY = new double[n];
var jitterZ = new double[n];
for (var i = 0; i < n; i++)
{
var cell = new float3(random.Next(-2000, 2000), random.Next(-2000, 2000), random.Next(-2000, 2000));
var hash = math.pcg3d(math.int3(cell));
var jitter = new float3(Unit(hash.x), Unit(hash.y), Unit(hash.z));
ids[i] = math.cellular(cell + jitter).id;
jitterX[i] = jitter.x;
jitterY[i] = jitter.y;
jitterZ[i] = jitter.z;
}
Assert.True(Math.Abs(Correlation(ids, jitterX)) < 0.05, "id correlates with jitter.x");
Assert.True(Math.Abs(Correlation(ids, jitterY)) < 0.05, "id correlates with jitter.y");
Assert.True(Math.Abs(Correlation(ids, jitterZ)) < 0.05, "id correlates with jitter.z");
}
[Fact]
public void JitterAndIdPassAChiSquareUniformityTest()
{
// F4/jitter uniformity: cellular hashes the winning cell's INTEGER coordinate (design.md,
// "Integer hashing"), never that coordinate's float32 bit pattern. Hashing the bit pattern of
// an integer-valued float shares more structure between adjacent cells (same exponent, a
// mantissa one increment apart) than the integers themselves do going into pcg3d's own
// mixing, and that structure survived into the jitter: Soul measured jitter.x's marginal
// chi-square statistic at 203.8 over 1,000,000 cells this way (15 degrees of freedom;
// critical value 37.7 at alpha = 0.001), soundly rejecting uniformity. 16 equal-width bins
// over [0, 1) give those same 15 degrees of freedom for jitter.x/y/z and id alike.
const int n = 200_000;
const int bins = 16;
const double criticalValue = 37.7; // chi-square, 15 df, alpha = 0.001
var jitterXCounts = new int[bins];
var jitterYCounts = new int[bins];
var jitterZCounts = new int[bins];
var idCounts = new int[bins];
var random = new System.Random(0xCE19);
for (var i = 0; i < n; i++)
{
var cell = new float3(random.Next(-100000, 100000), random.Next(-100000, 100000), random.Next(-100000, 100000));
var hash = math.pcg3d(math.int3(cell));
jitterXCounts[Bin(Unit(hash.x), bins)]++;
jitterYCounts[Bin(Unit(hash.y), bins)]++;
jitterZCounts[Bin(Unit(hash.z), bins)]++;
var idHash = math.pcg4d(new int4(math.int3(cell), 0));
idCounts[Bin(Unit(idHash.w), bins)]++;
}
var expected = (double)n / bins;
AssertUniform(jitterXCounts, expected, criticalValue, "jitter.x");
AssertUniform(jitterYCounts, expected, criticalValue, "jitter.y");
AssertUniform(jitterZCounts, expected, criticalValue, "jitter.z");
AssertUniform(idCounts, expected, criticalValue, "id");
}
private static int Bin(float value, int binCount) => Math.Clamp((int)(value * binCount), 0, binCount - 1);
private static void AssertUniform(int[] counts, double expected, double criticalValue, string name)
{
var chiSquare = 0.0;
foreach (var count in counts)
{
var delta = count - expected;
chiSquare += delta * delta / expected;
}
Assert.True(chiSquare < criticalValue, $"{name} chi-square {chiSquare:G6} exceeds critical value {criticalValue} ({counts.Length - 1} df)");
}
[Fact]
public void NoNaNAtIntegerPointsNear2To24()
{
// F3: past |p| of about 2^23, distinct cells' jittered features can round to the same
// float32 value, making F2 = 0 reachable. Soul measured 4096 of 4096 integer points at 2^24
// producing NaN before the guard; this pins the fix over the same sample.
const float scale = 16777216.0f; // 2^24
for (var ix = 0; ix < 16; ix++)
for (var iy = 0; iy < 16; iy++)
for (var iz = 0; iz < 16; iz++)
{
var p = new float3(scale + ix, scale + iy, scale + iz);
var result = math.cellular(p);
Assert.True(
float.IsFinite(result.nearest.x) && float.IsFinite(result.nearest.y) &&
float.IsFinite(result.nearest.z) && float.IsFinite(result.nearest.w),
$"nearest not finite at {p.x},{p.y},{p.z}: {result.nearest.x},{result.nearest.y},{result.nearest.z},{result.nearest.w}");
Assert.True(
float.IsFinite(result.edge.x) && float.IsFinite(result.edge.y) &&
float.IsFinite(result.edge.z) && float.IsFinite(result.edge.w),
$"edge not finite at {p.x},{p.y},{p.z}: {result.edge.x},{result.edge.y},{result.edge.z},{result.edge.w}");
Assert.True(float.IsFinite(result.id), $"id not finite at {p.x},{p.y},{p.z}");
}
}
[Fact]
public void CellularUnitMapsAllOnesBitsBelowOne()
{
// F7: cellular_unit(0xFFFFFFFF) must land in [0, 1), not round up to exactly 1.0f.
var mapped = math.cellular_unit(-1); // -1's bit pattern is 0xFFFFFFFF.
Assert.True(mapped >= 0.0f && mapped < 1.0f, $"cellular_unit(0xFFFFFFFF) = {mapped}");
}
[Fact]
public void ProductionSearchMatchesTheBruteForceOracleExactly()
{
// F6, part 2: the production 5x5x5 search must equal the independent 7x7x7 oracle exactly
// (to the oracle's own float32 precision), over a large seeded sweep in the domain Soul
// measured ([-50, 50]^3, where 3x3x3 missed 9 of 400,000 points).
var random = new System.Random(0xCE16);
for (var i = 0; i < 20000; i++)
{
var p = new float3(random.NextSingle() * 100.0f - 50.0f, random.NextSingle() * 100.0f - 50.0f, random.NextSingle() * 100.0f - 50.0f);
var result = math.cellular(p);
var (f1, f2) = FindF1F2(p);
Assert.Equal(f1, result.nearest.w, precision: 4);
Assert.Equal(f2 - f1, result.edge.w, precision: 4);
}
}
[Fact]
public void SearchReachesEveryAxisAlignedRadiusTwoSlice()
{
// F6, part 2 (Soul finding): a single-slice loop-bound mutant on one axis (e.g. dz < 2
// instead of dz <= 2, silently dropping the +2 layer on that axis alone) still passed the
// random 20,000-point sweeps above. That is not those sweeps being weak; it is the
// exactness proof's own margin (F1, F2 always < sqrt(3), a radius-2 cell's lower bound is
// max(0, 2-1) = 1, not "exactly 2" as an earlier version of this comment claimed) making a
// genuine win from any single |offset component| = 2 cell rare over uniformly random points
// (Soul, cut 2a-i: measured about 1 hit in 28,000 samples, not the "zero in 3,000,000" an
// earlier version of this comment claimed).
// Hunting for one by choosing candidate query POINTS is the wrong end of the search; this
// instead scans candidate integer CELLS and, for each, builds the query point most likely to
// realize that cell's neighbour as the winner: the closest point of the query cell's own unit
// cube to the target neighbour's actual jittered feature. That construction finds a real
// example for every one of the six axis-aligned radius-2 offsets within a couple thousand
// candidate cells, deterministically (pcg3d has no randomness to get lucky or unlucky with).
foreach (var offset in AxisAlignedRadiusTwoOffsets)
{
var p = FindPointRealizingOffset(offset);
var (f1, f2) = FindF1F2(p); // independent 7x7x7 oracle
var result = math.cellular(p);
Assert.Equal(f1, result.nearest.w, precision: 3);
Assert.Equal(f2 - f1, result.edge.w, precision: 3);
}
}
// The six axis-aligned radius-2 offsets (F6, part 2): shared by
// SearchReachesEveryAxisAlignedRadiusTwoSlice and, via FindPointRealizingOffset, by
// HlslSourceCompatibilityTests' mirror-comparison inputs, so both pin the same construction.
internal static readonly float3[] AxisAlignedRadiusTwoOffsets =
{
new(2.0f, 0.0f, 0.0f), new(-2.0f, 0.0f, 0.0f),
new(0.0f, 2.0f, 0.0f), new(0.0f, -2.0f, 0.0f),
new(0.0f, 0.0f, 2.0f), new(0.0f, 0.0f, -2.0f),
};
// Hunting for a query point where a radius-2 offset wins F1 or F2 by choosing candidate query
// POINTS is the wrong end of the search (Soul, cut 2a-i: about 1 hit in 28,000 uniformly random
// points); this instead scans candidate integer CELLS and, for each, builds the query point most
// likely to realize that cell's neighbour as the winner: the closest point of the query cell's
// own unit cube to the target neighbour's actual jittered feature. That construction finds a
// real example for every one of the six axis-aligned radius-2 offsets within a couple thousand
// candidate cells, deterministically (pcg3d has no randomness to get lucky or unlucky with).
internal static float3 FindPointRealizingOffset(float3 offset)
{
for (var qz = -80; qz <= 80; qz++)
for (var qy = -80; qy <= 80; qy++)
for (var qx = -80; qx <= 80; qx++)
{
var q = new float3(qx, qy, qz);
var targetFeature = FeaturePoint(q + offset);
var p = new float3(
Math.Clamp(targetFeature.x, q.x, q.x + 0.999f),
Math.Clamp(targetFeature.y, q.y, q.y + 0.999f),
Math.Clamp(targetFeature.z, q.z, q.z + 0.999f));
var (f1, f2) = FindF1F2(p); // independent 7x7x7 oracle
var d = math.distance(p, targetFeature);
if (MathF.Abs(d - f1) > 1.0e-4f && MathF.Abs(d - f2) > 1.0e-4f)
continue; // this cell didn't realize the target offset as F1 or F2; try another.
return p;
}
throw new InvalidOperationException($"could not construct a query point where offset {offset.x},{offset.y},{offset.z} ever wins F1 or F2");
}
[Fact]
public void PrunedSearchIsBitIdenticalToTheUnprunedFiveCubedSearchBelowTwoToTheTwentyFifth()
{
// F6, part 2: the lower-bound prune inside math.cellular's loop must never change the
// result. Compares against an independent, unpruned 5x5x5 reference that visits every cell
// in the same fixed order, over a large seeded sweep in [-50, 50]^3.
//
// This bit-identity has a domain: it holds for |p| < 2^25 (design.md, "Precision domain").
// Above that, cell + offset starts to round in float32, so the pruned loop's recomputed
// neighbour and the unpruned loop's neighbour at the same nominal (dx, dy, dz) can land on
// different actual cells, and the two searches diverge (Soul, cut 2a-i: 180 of 1,100,000
// points differ at |p| ~= 3.4e7). [-50, 50]^3 is comfortably inside the proven domain.
var random = new System.Random(0xCE18);
for (var i = 0; i < 20000; i++)
{
var p = new float3(random.NextSingle() * 100.0f - 50.0f, random.NextSingle() * 100.0f - 50.0f, random.NextSingle() * 100.0f - 50.0f);
var pruned = math.cellular(p);
var unpruned = UnprunedCellular5x5x5(p);
Assert.Equal(unpruned.nearest.x, pruned.nearest.x);
Assert.Equal(unpruned.nearest.y, pruned.nearest.y);
Assert.Equal(unpruned.nearest.z, pruned.nearest.z);
Assert.Equal(unpruned.nearest.w, pruned.nearest.w);
Assert.Equal(unpruned.edge.x, pruned.edge.x);
Assert.Equal(unpruned.edge.y, pruned.edge.y);
Assert.Equal(unpruned.edge.z, pruned.edge.z);
Assert.Equal(unpruned.edge.w, pruned.edge.w);
Assert.Equal(unpruned.id, pruned.id);
}
}
// Independent of math.cellular's own pruning: visits every cell of the 5x5x5 block, in the same
// fixed order, unconditionally.
private static CultCellular UnprunedCellular5x5x5(float3 p)
{
var cell = new float3(MathF.Floor(p.x), MathF.Floor(p.y), MathF.Floor(p.z));
var f1 = 1.0e30f;
var f2 = 1.0e30f;
var c1 = float3.zero;
var c2 = float3.zero;
var cellId = float3.zero;
for (var dz = -2; dz <= 2; dz++)
for (var dy = -2; dy <= 2; dy++)
for (var dx = -2; dx <= 2; dx++)
{
var neighbor = cell + new float3(dx, dy, dz);
var feature = FeaturePoint(neighbor);
var d = math.distance(p, feature);
if (d < f1) { f2 = f1; c2 = c1; f1 = d; c1 = feature; cellId = neighbor; }
else if (d < f2) { f2 = d; c2 = feature; }
}
var grad1 = f1 > 0.0f ? (p - c1) / f1 : float3.zero;
var grad2 = f2 > 0.0f ? (p - c2) / f2 : float3.zero;
var idHash = math.pcg4d(new int4(math.int3(cellId), 0));
return new CultCellular(new float4(grad1, f1), new float4(grad2 - grad1, f2 - f1), Unit(idHash.w));
}
private static double Correlation(double[] a, double[] b)
{
var n = a.Length;
var meanA = 0.0;
var meanB = 0.0;
for (var i = 0; i < n; i++) { meanA += a[i]; meanB += b[i]; }
meanA /= n;
meanB /= n;
var cov = 0.0;
var varA = 0.0;
var varB = 0.0;
for (var i = 0; i < n; i++)
{
var da = a[i] - meanA;
var db = b[i] - meanB;
cov += da * db;
varA += da * da;
varB += db * db;
}
return cov / Math.Sqrt(varA * varB);
}
// ---- shared fixtures: an independent brute-force Worley search over the same public pcg3d
// hash cellular itself is documented to use (design.md, "Integer hashing"), so these tests do
// not anchor to cellular's own control flow or iteration order. FindF1F2 searches 7x7x7 (radius
// 3), strictly wider than production's own proven-sufficient 5x5x5 (F6), so it is a true
// independent oracle rather than a second copy of the algorithm under test. ----
internal static float3 FeaturePoint(float3 cell)
{
var hash = math.pcg3d(math.int3(cell));
var jitter = new float3(Unit(hash.x), Unit(hash.y), Unit(hash.z));
return cell + jitter;
}
// Mirrors math.cellular_unit's top-24-bit [0, 1) mapping (F7) exactly, so this fixture never
// drifts from production's own jitter/id mapping.
private static float Unit(int bits) => ((uint)bits >> 8) * (1.0f / 16777216.0f);
private static (float f1, float f2) FindF1F2(float3 p)
{
var (f1, f2, _) = FindF1F2F3(p);
return (f1, f2);
}
// Same independent 7x7x7 oracle as FindF1F2, plus the third-nearest distance. Central-difference
// tests need this: excluding only the F1=F2 gap is not enough, because a small F2=F3 gap is also
// a kink (the identity of the second-nearest feature swaps there), and hitting one produces a
// large, curvature-driven finite-difference error unrelated to eps or float noise (Soul, cut
// 2a-i: an unfiltered large sample hit an edge error of 0.78 this way, versus ~2e-4 once F2=F3
// is also excluded).
private static (float f1, float f2, float f3) FindF1F2F3(float3 p)
{
var cell = new float3(MathF.Floor(p.x), MathF.Floor(p.y), MathF.Floor(p.z));
var f1 = float.PositiveInfinity;
var f2 = float.PositiveInfinity;
var f3 = float.PositiveInfinity;
for (var dz = -3; dz <= 3; dz++)
for (var dy = -3; dy <= 3; dy++)
for (var dx = -3; dx <= 3; dx++)
{
var d = math.distance(p, FeaturePoint(cell + new float3(dx, dy, dz)));
if (d < f1) { f3 = f2; f2 = f1; f1 = d; }
else if (d < f2) { f3 = f2; f2 = d; }
else if (d < f3) { f3 = d; }
}
return (f1, f2, f3);
}
private static (float3 offset, float f1, float f2) FindNearestOffset(float3 p)
{
var cell = new float3(MathF.Floor(p.x), MathF.Floor(p.y), MathF.Floor(p.z));
var f1 = float.PositiveInfinity;
var f2 = float.PositiveInfinity;
var bestOffset = float3.zero;
for (var dz = -2; dz <= 2; dz++)
for (var dy = -2; dy <= 2; dy++)
for (var dx = -2; dx <= 2; dx++)
{
var offset = new float3(dx, dy, dz);
var d = math.distance(p, FeaturePoint(cell + offset));
if (d < f1) { f2 = f1; f1 = d; bestOffset = offset; }
else if (d < f2) { f2 = d; }
}
return (bestOffset, f1, f2);
}
private static float3 FindNearestCell(float3 p)
{
var cell = new float3(MathF.Floor(p.x), MathF.Floor(p.y), MathF.Floor(p.z));
var best = float.PositiveInfinity;
var bestCell = cell;
for (var dz = -1; dz <= 1; dz++)
for (var dy = -1; dy <= 1; dy++)
for (var dx = -1; dx <= 1; dx++)
{
var neighbour = cell + new float3(dx, dy, dz);
var d = math.distance(p, FeaturePoint(neighbour));
if (d < best) { best = d; bestCell = neighbour; }
}
return bestCell;
}
private static void AssertGradientMatchesCentralDifference(Func<float3, float> scalarField, float3 p, float eps, float4 valueAndGradient, float tolerance)
{
var numeric = new float3(
(scalarField(p + new float3(eps, 0.0f, 0.0f)) - scalarField(p - new float3(eps, 0.0f, 0.0f))) / (2.0f * eps),
(scalarField(p + new float3(0.0f, eps, 0.0f)) - scalarField(p - new float3(0.0f, eps, 0.0f))) / (2.0f * eps),
(scalarField(p + new float3(0.0f, 0.0f, eps)) - scalarField(p - new float3(0.0f, 0.0f, eps))) / (2.0f * eps));
AssertWithinTolerance(numeric.x, valueAndGradient.x, tolerance);
AssertWithinTolerance(numeric.y, valueAndGradient.y, tolerance);
AssertWithinTolerance(numeric.z, valueAndGradient.z, tolerance);
}
// A fixed decimal-place Assert.Equal snaps two numbers that agree to within a few times 1e-5 to
// different rounded digits whenever their true difference straddles a rounding boundary (seen:
// 0.824987829 vs 0.825006127, agreeing to 2e-5, rounding to 0.82 and 0.83). An absolute-tolerance
// check states the actual claim (central differences agree with the analytic gradient within the
// finite-difference step's own truncation error) without that artifact.
private static void AssertWithinTolerance(float expected, float actual, float tolerance) =>
Assert.True(MathF.Abs(expected - actual) <= tolerance, $"expected {expected}, actual {actual}, tolerance {tolerance}");
private static float3 RandomUnit(System.Random random) =>
new(random.NextSingle() * 2.0f - 1.0f, random.NextSingle() * 2.0f - 1.0f, random.NextSingle() * 2.0f - 1.0f);
private static float RandomScalar(System.Random random) => random.NextSingle() * 20.0f - 10.0f;
}

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

@ -49,6 +49,8 @@ public sealed class HlslMirrorTests
"cultmath_value_noise_bicubic",
"cultmath_value_noise_texture",
"cultmath_value_noise_texture_bicubic",
"cultmath_smin_grad",
"cultmath_cellular",
};
foreach (var symbol in requiredSymbols)

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

@ -25,6 +25,8 @@ public sealed class HlslSourceCompatibilityTests
Assert.NotNull(ShaderFunction(shader, "cultmath_snoise", typeof(float3)));
Assert.NotNull(ShaderFunction(shader, "cultmath_snoise", typeof(float2)));
Assert.NotNull(ShaderFunction(shader, "cultmath_value_noise_bicubic", typeof(float2)));
Assert.NotNull(ShaderFunction(shader, "cultmath_smin_grad", typeof(float4), typeof(float4), typeof(float)));
Assert.NotNull(ShaderFunction(shader, "cultmath_cellular", typeof(float3)));
}
// HLSL functions carry no access modifier, so they compile as private instance methods.
@ -62,6 +64,41 @@ public sealed class HlslSourceCompatibilityTests
for (var k = 0; k < 64; k++) { var values = Enumerable.Range(0, 32).Select(_ => Next()).ToArray(); cases.Add(i => values[i]); }
for (var k = 0; k < 16; k++) { var values = Enumerable.Range(0, 32).Select(_ => Next()).ToArray(); cases.Add(i => values[i & ~1]); }
// Points aimed squarely at cellular's HLSL mirror (F6/F3 Soul findings): the generic random
// and special-value cases above pass ~92 inputs through cultmath_cellular and never happen to
// hit a point where a radius-2 cell wins (about 3.5e-5 of points, Soul measured) or an
// integer point at or past 2^23, so an HLSL mutant of the search radius, the prune, or the
// removed F2=0 guard can survive this test even though it changes cellular's real output.
void AddPointCase(float3 point)
{
var values = new[] { point.x, point.y, point.z };
cases.Add(i => i < values.Length ? values[i] : 0.0f);
}
// The six constructed points from CellularAndSminGradTests.SearchReachesEveryAxisAlignedRadiusTwoSlice:
// each one realizes a different axis-aligned radius-2 offset as F1 or F2, which is exactly
// the case an `if (lowerBound >= f2)` prune written as `>= f1`, or a search radius narrowed on
// one axis (`dz < 2`, `dx` starting at -1), changes.
foreach (var offset in CellularAndSminGradTests.AxisAlignedRadiusTwoOffsets)
AddPointCase(CellularAndSminGradTests.FindPointRealizingOffset(offset));
// Integer points at 2^24, several, both signs: past |p| ~= 2^23 distinct cells' jittered
// features can round onto the same float32 value (design.md, "Precision domain"), making
// F2 = 0 reachable and exercising the grad2 = f2 > 0 ? ... : zero guard this test would
// otherwise never touch.
const float twoTo24 = 16777216.0f;
AddPointCase(new float3(twoTo24, twoTo24, twoTo24));
AddPointCase(new float3(-twoTo24, -twoTo24, -twoTo24));
AddPointCase(new float3(twoTo24 + 3.0f, -twoTo24 + 5.0f, twoTo24 - 7.0f));
AddPointCase(new float3(-twoTo24 + 11.0f, twoTo24 - 13.0f, -twoTo24 + 17.0f));
// One exact F1 = F2 tie point: the midpoint of two adjacent cells' feature points is
// bit-exact-equidistant from both by construction (negating a vector does not change its
// length, and p - a = -(p - b) here by symmetry), without needing to search for one.
var tieMidpoint = (CellularAndSminGradTests.FeaturePoint(new float3(0.0f, 0.0f, 0.0f))
+ CellularAndSminGradTests.FeaturePoint(new float3(1.0f, 0.0f, 0.0f))) * 0.5f;
AddPointCase(tieMidpoint);
var mismatches = new List<string>();
var compared = new HashSet<string>();
foreach (var mirror in shaderType.GetMethods(BindingFlags.Instance | BindingFlags.NonPublic).Where(m => m.Name.StartsWith("cultmath_")))
@ -82,17 +119,93 @@ public sealed class HlslSourceCompatibilityTests
: Activator.CreateInstance(t, Components(t).Select(f => Scalar(f.FieldType, input(slot++))).ToArray())!).ToArray();
var expected = Values(counterpart.Invoke(null, args)!);
var actual = Values(mirror.Invoke(shader, args)!);
if (expected.Length != actual.Length || expected.Zip(actual).Any(p => p.First is float a && p.Second is float b
? !(float.IsNaN(a) && float.IsNaN(b)) && BitConverter.SingleToInt32Bits(a) != BitConverter.SingleToInt32Bits(b)
: !p.First.Equals(p.Second)))
if (!BitwiseEqual(expected, actual))
{
mismatches.Add($"{mirror.Name}({string.Join(", ", args.Select(a => string.Join(" ", Values(a))))}): C# {string.Join(" ", expected)}, mirror {string.Join(" ", actual)}");
mismatches.Add($"{mirror.Name}({string.Join(", ", args.Select(a => string.Join(" ", Values(a).Select(v => v.Value))))}): C# {string.Join(" ", expected.Select(v => v.Value))}, mirror {string.Join(" ", actual.Select(v => v.Value))}");
}
}
}
Assert.True(mismatches.Count == 0, $"{mismatches.Count} mismatches:{Environment.NewLine}{string.Join(Environment.NewLine, mismatches.Take(12))}");
Assert.Equal(28, compared.Count);
Assert.Equal(30, compared.Count);
}
/// <summary>
/// A hand-built pair of structs standing in for a mirror struct return and its C# counterpart
/// (invariant 8's one struct return shape, e.g. <see cref="CultCellular"/>): same field shape,
/// different concrete type, exactly like the real mirror comparison. Proves the field-by-field,
/// bit-for-bit walk actually rejects a struct whose field differs; without this, the walk could
/// have been comparing lengths only, or skipping the nested vector's own components.
/// </summary>
private struct StructMirrorProbeExpected { public float4 first; public float second; }
private struct StructMirrorProbeActual { public float4 first; public float second; }
[Fact]
public void StructReturnComparisonRejectsAMismatchedField()
{
var expected = new StructMirrorProbeExpected { first = new float4(1.0f, 2.0f, 3.0f, 4.0f), second = 5.0f };
var matching = new StructMirrorProbeActual { first = new float4(1.0f, 2.0f, 3.0f, 4.0f), second = 5.0f };
var mismatchedLeaf = new StructMirrorProbeActual { first = new float4(1.0f, 2.0f, 3.0f, 4.0f), second = 5.5f };
var mismatchedNestedComponent = new StructMirrorProbeActual { first = new float4(1.0f, 2.0f, -3.0f, 4.0f), second = 5.0f };
Assert.True(BitwiseEqual(Values(expected), Values(matching)));
Assert.False(BitwiseEqual(Values(expected), Values(mismatchedLeaf)));
Assert.False(BitwiseEqual(Values(expected), Values(mismatchedNestedComponent)));
}
/// <summary>
/// Struct shapes standing in for CultCellular's own field order (<c>nearest; edge; id</c>), used
/// to falsify two named mutants (F2): H2 declares the mirror struct with its first two fields
/// swapped (<c>edge; nearest; id</c>) and swaps the assignments to match, so every leaf VALUE
/// still lands in the same POSITION as the correct struct; only the field NAME at each position
/// differs. H3 renames and retypes every field (id becomes an int carrying the same bits) while
/// keeping the original position order, so only NAME and TYPE differ, never position.
/// </summary>
private struct NearestEdgeIdExpected { public float4 nearest; public float4 edge; public float id; }
private struct EdgeNearestIdSwappedNames { public float4 edge; public float4 nearest; public float id; }
private struct RenamedAndRetypedFields { public float4 a; public float4 b; public int c; }
[Fact]
public void StructReturnComparisonRejectsFieldsSwappedByNameEvenWhenValuesAlign()
{
// H2: the mutant struct's field ORDER is (edge, nearest, id), and its assignments are
// swapped to match, so declaration-position values are identical to the correct struct's.
// Only checking names (not just position) can tell these apart.
var expected = new NearestEdgeIdExpected
{
nearest = new float4(1.0f, 2.0f, 3.0f, 4.0f),
edge = new float4(5.0f, 6.0f, 7.0f, 8.0f),
id = 9.0f,
};
var swapped = new EdgeNearestIdSwappedNames
{
edge = new float4(1.0f, 2.0f, 3.0f, 4.0f), // holds what should be "nearest"'s value.
nearest = new float4(5.0f, 6.0f, 7.0f, 8.0f), // holds what should be "edge"'s value.
id = 9.0f,
};
Assert.False(BitwiseEqual(Values(expected), Values(swapped)));
}
[Fact]
public void StructReturnComparisonRejectsRenamedAndRetypedFields()
{
// H3: field order matches the correct struct exactly, but every field is renamed, and id is
// retyped from float to int carrying the same bit pattern the correct id would round to.
var expected = new NearestEdgeIdExpected
{
nearest = new float4(1.0f, 2.0f, 3.0f, 4.0f),
edge = new float4(5.0f, 6.0f, 7.0f, 8.0f),
id = 9.0f,
};
var renamed = new RenamedAndRetypedFields
{
a = new float4(1.0f, 2.0f, 3.0f, 4.0f),
b = new float4(5.0f, 6.0f, 7.0f, 8.0f),
c = BitConverter.SingleToInt32Bits(9.0f),
};
Assert.False(BitwiseEqual(Values(expected), Values(renamed)));
}
// Integer arguments take the bit pattern of the float input, so specials become 0, 0x80000000, NaN bits, ...
@ -101,10 +214,41 @@ public sealed class HlslSourceCompatibilityTests
private static object Scalar(Type t, float f) => t == typeof(float) ? f
: t == typeof(int) ? BitConverter.SingleToInt32Bits(f) : (object)(uint)BitConverter.SingleToInt32Bits(f);
// Scoped to scalar-typed fields: this only ever builds function ARGUMENTS (float2/3/4, int2/3/4, ...),
// never a struct-of-struct return, so it must not recurse.
private static FieldInfo[] Components(Type t) => t.GetFields().Where(f => !f.IsStatic && Scalars.Contains(f.FieldType)).ToArray();
private static object[] Values(object value) => Scalars.Contains(value.GetType()) ? new[] { value }
: Components(value.GetType()).Select(field => field.GetValue(value)!).ToArray();
/// <summary>
/// One flattened scalar leaf of a (possibly nested) return value: its declaration path (field
/// names joined by '.'), its declared type, and its value. Carrying path and type, not just the
/// value, is what lets the struct-return comparison reject a field that was renamed or retyped
/// but landed in the same position (F2): a positional-only comparison cannot tell "edge" holding
/// nearest's value apart from "nearest" holding it.
/// </summary>
private readonly record struct Leaf(string Path, Type Type, object Value);
/// <summary>
/// Flattens a return value to its scalar leaves, recursing through struct fields of any type
/// (invariant 8's one struct return shape, e.g. <see cref="CultCellular"/>'s float4/float4/float
/// fields) so a struct return is compared field by field, down to bit-for-bit scalars, exactly
/// like a bare vector return already is. Each leaf keeps the field name (and its ancestors') and
/// declared type it was read from.
/// </summary>
private static Leaf[] Values(object value, string path = "") => Scalars.Contains(value.GetType())
? new[] { new Leaf(path, value.GetType(), value) }
: value.GetType().GetFields().Where(f => !f.IsStatic)
.SelectMany(f => Values(f.GetValue(value)!, path.Length == 0 ? f.Name : path + "." + f.Name))
.ToArray();
// Compares leaves by NAME and TYPE as well as value, in declaration order, so a mirror struct
// whose fields are reordered, renamed, or retyped relative to its C# counterpart fails here even
// when the flattened values happen to line up positionally (F2's H2/H3 mutants).
private static bool BitwiseEqual(Leaf[] expected, Leaf[] actual) =>
expected.Length == actual.Length && expected.Zip(actual).All(p =>
p.First.Path == p.Second.Path && p.First.Type == p.Second.Type &&
(p.First.Value is float a && p.Second.Value is float b
? (float.IsNaN(a) && float.IsNaN(b)) || BitConverter.SingleToInt32Bits(a) == BitConverter.SingleToInt32Bits(b)
: p.First.Value.Equals(p.Second.Value)));
/// <summary>The documented HLSL-to-C# transformations, and nothing else.</summary>
internal static string TransformHlslToCSharp(string hlsl)
@ -122,7 +266,13 @@ public sealed class HlslSourceCompatibilityTests
// 4. Floating literals gain the `f` suffix; an unsuffixed C# literal is a double.
source = Regex.Replace(source, @"(?<![\w.])(\d+\.\d*(?:[eE][+-]?\d+)?|\d+[eE][+-]?\d+)(?![\w.])", "$1f");
// 5. The file body is wrapped in a class; C# has no free functions.
// 5. Struct fields gain `public`; HLSL has no field access-modifier concept, and a C# struct's
// fields default to private, which would hide them from the mirror test's reflection-based
// field comparison (invariant 8's one struct return shape, e.g. CultCellular).
source = Regex.Replace(source, @"(?ms)(struct\s+\w+\s*\{)(.*?)(\};)", m =>
m.Groups[1].Value + Regex.Replace(m.Groups[2].Value, @"(?m)^(\s*)(\w+\s+\w+;)", "$1public $2") + m.Groups[3].Value);
// 6. The file body is wrapped in a class; C# has no free functions.
return "using CultMath;\nusing static CultMath.math;\nnamespace CultMathHlsl;\npublic class HlslShader\n{\n" + source + "\n}\n";
}

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

@ -313,6 +313,88 @@ int4 cultmath_pcg4d(int4 value)
int4 cultmath_pcg4d(float4 value) { return cultmath_pcg4d(int4((int)asuint(value.x), (int)asuint(value.y), (int)asuint(value.z), (int)asuint(value.w))); }
// Inigo Quilez, "smooth minimum": quadratic-polynomial smin carried to value-and-gradient form. See
// math.smin_grad's comment (math.cs) for why the gradient blend is exact, not an approximation.
float4 cultmath_smin_grad(float4 a, float4 b, float k)
{
float h = saturate(0.5 + 0.5 * (b.w - a.w) / k);
float value = lerp(b.w, a.w, h) - k * h * (1.0 - h);
float3 gradient = lerp(b.xyz, a.xyz, h);
return float4(gradient, value);
}
// Invariant 8's one struct return shape (design.md, "HLSL Target: Source Transformations"):
// compared field by field and bit for bit against CultMath.CultCellular, same field order.
struct CultCellular
{
float4 nearest;
float4 edge;
float id;
};
// Worley, "A Cellular Texture Basis Function" (SIGGRAPH 1996). See math.cellular's comment
// (math.cs) for the hashing, gradient, search-radius exactness proof, id source and
// degenerate-point conventions this mirrors exactly, including the 5x5x5 lower-bound prune.
CultCellular cultmath_cellular(float3 p)
{
float3 cell = floor(p);
float f1 = 1.0e30;
float f2 = 1.0e30;
float3 c1 = float3(0.0, 0.0, 0.0);
float3 c2 = float3(0.0, 0.0, 0.0);
float3 cellId = float3(0.0, 0.0, 0.0);
for (int dz = -2; dz <= 2; dz++)
{
float lz = max(0.0, abs((float)dz) - 1.0);
for (int dy = -2; dy <= 2; dy++)
{
float ly = max(0.0, abs((float)dy) - 1.0);
for (int dx = -2; dx <= 2; dx++)
{
float lx = max(0.0, abs((float)dx) - 1.0);
float lowerBound = sqrt(lx * lx + ly * ly + lz * lz);
if (lowerBound >= f2)
continue;
float3 neighbor = cell + float3(dx, dy, dz);
int3 hash = cultmath_pcg3d(int3(neighbor));
float3 jitter = float3(
((uint)hash.x >> 8) * (1.0 / 16777216.0),
((uint)hash.y >> 8) * (1.0 / 16777216.0),
((uint)hash.z >> 8) * (1.0 / 16777216.0));
float3 feature = neighbor + jitter;
float d = cultmath_distance(p, feature);
if (d < f1)
{
f2 = f1; c2 = c1;
f1 = d; c1 = feature; cellId = neighbor;
}
else if (d < f2)
{
f2 = d; c2 = feature;
}
}
}
}
// F1 = 0 and F2 = 0 share the same degenerate-gradient guard; see math.cellular's comment
// (math.cs) for why F2 = 0 is reachable past roughly |p| = 2^23.
float3 grad1 = f1 > 0.0 ? (p - c1) / f1 : float3(0.0, 0.0, 0.0);
float3 grad2 = f2 > 0.0 ? (p - c2) / f2 : float3(0.0, 0.0, 0.0);
// id comes from pcg4d over the winning cell's integer coordinate alone (same int-hashing rule
// as the jitter above), never from the pcg3d hash that produced its jitter (F4).
int4 idHash = cultmath_pcg4d(int4(int3(cellId), 0));
CultCellular result;
result.nearest = float4(grad1, f1);
result.edge = float4(grad2 - grad1, f2 - f1);
result.id = ((uint)idHash.w >> 8) * (1.0 / 16777216.0);
return result;
}
float cultmath_value_noise(float2 position)
{
float2 cell = floor(position);