diff --git a/packages/cultmath/docs/design.md b/packages/cultmath/docs/design.md index f5f5d5a5..30cef06c 100644 --- a/packages/cultmath/docs/design.md +++ b/packages/cultmath/docs/design.md @@ -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", + ), 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 diff --git a/packages/cultmath/shaders/CultMath.hlsl b/packages/cultmath/shaders/CultMath.hlsl index 6179d45a..5c00bda0 100644 --- a/packages/cultmath/shaders/CultMath.hlsl +++ b/packages/cultmath/shaders/CultMath.hlsl @@ -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); diff --git a/packages/cultmath/src/CultMath/CultCellular.cs b/packages/cultmath/src/CultMath/CultCellular.cs new file mode 100644 index 00000000..f856236c --- /dev/null +++ b/packages/cultmath/src/CultMath/CultCellular.cs @@ -0,0 +1,26 @@ +namespace CultMath; + +/// +/// The value-and-gradient result of (invariant 8's one struct +/// return shape). Every field follows the same float4(gradient.xyz, value.w) layout as a +/// bare value-and-gradient return. +/// +public struct CultCellular +{ + /// (∇F1, F1): the gradient and distance to the nearest feature point. + public float4 nearest; + + /// (∇F2 - ∇F1, F2 - F1): for lineae along cell borders. + public float4 edge; + + /// The nearest cell's pcg4d hash, mapped to [0, 1); independent of the pcg3d hash that + /// produced that cell's jitter. + public float id; + + public CultCellular(float4 nearest, float4 edge, float id) + { + this.nearest = nearest; + this.edge = edge; + this.id = id; + } +} diff --git a/packages/cultmath/src/CultMath/CultCellular.cs.meta b/packages/cultmath/src/CultMath/CultCellular.cs.meta new file mode 100644 index 00000000..8c9ed1ae --- /dev/null +++ b/packages/cultmath/src/CultMath/CultCellular.cs.meta @@ -0,0 +1,2 @@ +fileFormatVersion: 2 +guid: 76be9498b63e410891715b3ec965caf8 diff --git a/packages/cultmath/src/CultMath/math.cs b/packages/cultmath/src/CultMath/math.cs index b5d5339c..1b03f4ed 100644 --- a/packages/cultmath/src/CultMath/math.cs +++ b/packages/cultmath/src/CultMath/math.cs @@ -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); diff --git a/packages/cultmath/tests/CultMath.Tests/CellularAndSminGradTests.cs b/packages/cultmath/tests/CultMath.Tests/CellularAndSminGradTests.cs new file mode 100644 index 00000000..465f6961 --- /dev/null +++ b/packages/cultmath/tests/CultMath.Tests/CellularAndSminGradTests.cs @@ -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 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; +} diff --git a/packages/cultmath/tests/CultMath.Tests/HlslMirrorTests.cs b/packages/cultmath/tests/CultMath.Tests/HlslMirrorTests.cs index d73c6d7b..2b695fcb 100644 --- a/packages/cultmath/tests/CultMath.Tests/HlslMirrorTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/HlslMirrorTests.cs @@ -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) diff --git a/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs b/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs index 9426e4c8..fb269ee2 100644 --- a/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs +++ b/packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs @@ -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(); var compared = new HashSet(); 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); + } + + /// + /// 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. ): 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. + /// + 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))); + } + + /// + /// Struct shapes standing in for CultCellular's own field order (nearest; edge; id), used + /// to falsify two named mutants (F2): H2 declares the mirror struct with its first two fields + /// swapped (edge; nearest; id) 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. + /// + 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(); + /// + /// 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. + /// + private readonly record struct Leaf(string Path, Type Type, object Value); + + /// + /// Flattens a return value to its scalar leaves, recursing through struct fields of any type + /// (invariant 8's one struct return shape, e.g. '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. + /// + 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))); /// The documented HLSL-to-C# transformations, and nothing else. 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, @"(? + 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"; } diff --git a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl index 6179d45a..5c00bda0 100644 --- a/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl +++ b/packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl @@ -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);