Merge remote-tracking branch 'origin/main' into hands/pub-cut4

docs/cultcache-publication-cut.md

@@ -813,4 +813,16 @@ Recorded, not forked:
   a gap (`ApplyShardSnapshotResponseAsync` `:675-736`). It predates this campaign.
 - **AetheriaEve `Program.cs:3722`.** It reads a committed write's `InvalidOperationException` as a denial. Under
   Q-P5 A, Cut 4's log failures cannot reach it. Observer and `OnUpdate` exceptions still can, so the Cuts 1-2
-  recording stands.
+  recording stands.- **Cut 4 Soul, recorded 2026-09-30 (Self), not fixed in Cut 4:**
+  - **File log store is not crash-safe (hard prerequisite).** `CultNetFileShardMutationLogStore.Append` and
+    `CompactThrough` use `File.WriteAllBytes` (truncate, rewrite, no temp+rename, no fsync). A failure mid-write
+    truncates the log; every later append burns a sequence and a restart cannot construct the database. It also
+    rewrites the whole file per append, which Cut 4 now does under the cache gate (2.2 ms per write at 100
+    entries, 33 ms at 950, a 2.49 s stall). An atomic, append-only store must land before any live daemon
+    enables durable shard logs. Today only AetheriaEve's offline tools do; the live Aetheria daemon passes
+    `enableDurableShardLogs: false`.
+  - **A slow log append blocks cache readers** (`CultCache.Get` takes the gate). Accepted with a note in the
+    store's XML doc.
+  - **Crash between the store write and the log append** leaves the log lagging the store, and nothing detects
+    it; a restarted primary reuses that sequence (A2's divergence reached by a plain kill). The spec called this
+    lag detectable; for a crash it is not.

packages/cultmath/README.md

@@ -33,6 +33,8 @@ shader semantics.
 - `shaders/CultMath.hlsl`, a canonical HLSL mirror include for shader-side
   parity. HLSL already owns `float2`, `float3`, and `float4`; the include
   exposes `cultmath_*` functions for shared semantics HLSL does not provide.
+  It includes `shaders/CultMath.Phacelle.hlsl` (MPL-2.0) from the same
+  directory, so ship both files.
 - `Voronoi.SampleTones`, a C# batch surface that calls the Rust
   `cultmath-core` native kernel when `cultmath_core` is available and falls back
   to the managed parity path otherwise.

packages/cultmath/THIRD-PARTY-NOTICES.md

@@ -19,3 +19,23 @@ permutation and scale as `snoise(float3)`.
 - Copyright 2011 Ashima Arts.
 - Upstream: <https://github.com/ashima/webgl-noise>
 - License: MIT.
+
+## Phacelle Noise
+
+`math.phacelle` and `cultmath_phacelle` follow Rune Skovbo Johansen's Phacelle
+noise (Shadertoy `t3dyWl`; erosion filter `wXcfWn`; blog, March 2026), read
+from the C# port `PhacelleNoise` in `lpmitchell/AdvancedTerrainErosion`
+(`AdvancedTerrainErosion.cs`, converted to burstable C# by Luke Mitchell,
+2026). The cell weight, the phase blend and the normalization follow it. The
+3D-cell generalisation, the caller-supplied stripe wave vector, the exact
+gradient and the exact pruning are new here.
+
+- Copyright (c) 2025 Rune Skovbo Johansen.
+- Upstream: <https://github.com/lpmitchell/AdvancedTerrainErosion>
+- License: MPL-2.0 (<https://mozilla.org/MPL/2.0/>). MPL-2.0 is file-level, so it
+  covers exactly the files that carry the port and nothing else:
+  `src/CultMath/math.Phacelle.cs`, `src/CultMath/CultPhasor.cs`,
+  `shaders/CultMath.Phacelle.hlsl` and the Unity copy of that include. Each
+  carries the MPL-2.0 header and this attribution. The rest of CultMath is MIT,
+  so the `GameCult.Math` NuGet package declares `MIT AND MPL-2.0` and ships this
+  file. The Unity package as a whole is MPL-2.0.

packages/cultmath/docs/design.md

@@ -128,8 +128,8 @@ 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
+as its value. `CultCellular` (`math.cs`, `cellular`) and `CultPhasor`
+(`math.Phacelle.cs`, `phacelle`) are its consumers. A struct return with a mismatched, reordered-by-name, or
 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.
 
@@ -264,6 +264,41 @@ together, laid out as `float4(gradient.xyz, value.w)`, with no value-only twin
   which is the right domain for `pcg3d`/`pcg4d`'s O(1) hash inputs but would
   turn an unclamped octave count into a multi-billion-iteration loop for these
   two mirrors.
+- `phacelle(float3 p, float3 side, float offset, float normalization)` is Rune
+  Skovbo Johansen's Phacelle noise ("phase + cell": a stripe pattern built by
+  blending one cosine/sine wave per jittered cell), generalised to 3D cells. It
+  returns `CultPhasor { cos, sin }`, each `(∇, value)`. `side` is the stripe
+  wave vector, whose length is 2π times the stripe frequency; in 3D the
+  perpendicular to the flow is not unique, so the caller chooses it (on a unit
+  sphere, `cross(p̂, flow)`). `offset` is the phase in cycles. It visits 4×4×4
+  cells with `pcg3d`-hashed jitter in `[-0.5, 0.5)`, weights each wave by
+  `max(0, exp(-2d²) - 0.01111)`, and normalizes the blended phasor as upstream
+  does (length 1 where the raw length is at least `1 - normalization`, scaled
+  by `1 / (1 - normalization)` below). The gradient is exact, weight
+  derivatives and normalization included, rather than upstream's
+  `∇cos ≈ -sin·side`, which ignores the weight gradient. It is discontinuous
+  where the raw length crosses the normalization threshold and where a cell's
+  weight leaves its support (d² = 2.24995, distance 1.49998); the value is
+  continuous at both. At the support edge `∇w` jumps from `-4·0.01111·v` to 0,
+  a per-cell jump of length `0.0444·|v|` (0.0667 at `|v| = 1.5`); measured over
+  random points, the output gradient jumps by a median 1.5% of its length and
+  at most about 27%. `normalization` is meant to lie in `[0, 1]`: 1 or more
+  saturates (every output is stretched to length 1), a negative value scales
+  every output down by `1 / (1 - normalization)`, and the normalized gradient
+  scales as `1 / |I|` for the blended phasor `I`, so it grows where the blend
+  nearly cancels. The total weight `W` is positive for every `pcg3d` input (the
+  smallest found by hill-climbing is 0.056); `W = 0`, and a NaN output, needs
+  adversarial jitter no hash yields.
+  Cells whose per-axis lower bound `Σ max(0, |local - g| - 0.5)²` is at least
+  2.25 lie outside the support and have weight exactly 0, so they are skipped,
+  which changes no output bit. Cells whose bound sits just under 2.25 can still
+  carry a nonzero weight, so the cut cannot be lowered; `PhacelleTests` pins
+  the cut with points next to such cells, against an internal overload with the
+  prune disabled, and the HLSL mirror test pins it on the shader side. The prune
+  keeps about 45 of the 64 cells on average (a saving of about 30%); about 14
+  have a nonzero weight. The port sits in its own MPL-2.0 files
+  (`math.Phacelle.cs`, `CultPhasor.cs`, `CultMath.Phacelle.hlsl`); see
+  `THIRD-PARTY-NOTICES.md`. No 2D variant exists; nothing consumes one.
 
 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

packages/cultmath/scripts/build-unity-package.ps1

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

packages/cultmath/shaders/CultMath.Phacelle.hlsl

@@ -0,0 +1,99 @@
+#ifndef CULTMATH_PHACELLE_HLSL
+#define CULTMATH_PHACELLE_HLSL
+
+// This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0. If a copy of
+// the MPL was not distributed with this file, You can obtain one at https://mozilla.org/MPL/2.0/.
+//
+// Phacelle noise. Copyright (c) 2025 Rune Skovbo Johansen (Shadertoy t3dyWl; MPL-2.0), read from the
+// C# port PhacelleNoise in lpmitchell/AdvancedTerrainErosion, converted to burstable C# by Luke
+// Mitchell (2026). Modified by GameCult: 3D cells, a caller-supplied stripe wave vector, an exact
+// gradient and exact pruning. This file is an MPL-2.0 file; CultMath.hlsl, which includes it after
+// its own definitions of CULTMATH_TAU and cultmath_pcg3d, is MIT (see THIRD-PARTY-NOTICES.md).
+
+// Invariant 8's struct return shape for cultmath_phacelle, compared field by field and bit for bit
+// against CultMath.CultPhasor, same field order.
+struct CultPhasor
+{
+    float4 cos;
+    float4 sin;
+};
+
+// Johansen's Phacelle noise on 3D cells. See math.phacelle's comment (math.Phacelle.cs) for the stripe wave
+// vector, the weight, the normalization, the exact gradient and why the prune is exact. Mirrors
+// the C# body operation for operation; the C#-only internal overload that switches the prune off
+// has no shader twin.
+CultPhasor cultmath_phacelle(float3 p, float3 side, float offset, float normalization)
+{
+    float3 cell = floor(p);
+    float3 local = p - cell;
+    float phaseOffset = offset * CULTMATH_TAU;
+    float sumW = 0.0;
+    float sumCos = 0.0;
+    float sumSin = 0.0;
+    float3 gradW = float3(0.0, 0.0, 0.0);
+    float3 gradCos = float3(0.0, 0.0, 0.0);
+    float3 gradSin = float3(0.0, 0.0, 0.0);
+
+    for (int dz = -1; dz <= 2; dz++)
+    {
+        float gz = max(0.0, abs(local.z - (float)dz) - 0.5);
+        for (int dy = -1; dy <= 2; dy++)
+        {
+            float gy = max(0.0, abs(local.y - (float)dy) - 0.5);
+            for (int dx = -1; dx <= 2; dx++)
+            {
+                float gx = max(0.0, abs(local.x - (float)dx) - 0.5);
+                if (gx * gx + gy * gy + gz * gz >= 2.25)
+                    continue;
+
+                float3 gridStep = float3(dx, dy, dz);
+                int3 hash = cultmath_pcg3d(int3(cell + gridStep));
+                float3 jitter = float3(
+                    ((uint)hash.x >> 8) * (1.0 / 16777216.0) - 0.5,
+                    ((uint)hash.y >> 8) * (1.0 / 16777216.0) - 0.5,
+                    ((uint)hash.z >> 8) * (1.0 / 16777216.0) - 0.5);
+                float3 v = local - gridStep - jitter;
+                float falloff = exp(-2.0 * dot(v, v));
+                float w = max(0.0, falloff - 0.01111);
+                float3 dw = w > 0.0 ? (-4.0 * falloff) * v : float3(0.0, 0.0, 0.0);
+                float phase = dot(v, side) + phaseOffset;
+                float c = cos(phase);
+                float s = sin(phase);
+
+                sumW += w;
+                sumCos += w * c;
+                sumSin += w * s;
+                gradW += dw;
+                gradCos += dw * c - (w * s) * side;
+                gradSin += dw * s + (w * c) * side;
+            }
+        }
+    }
+
+    float invW = 1.0 / sumW;
+    float rawCos = sumCos * invW;
+    float rawSin = sumSin * invW;
+    float3 gradRawCos = (gradCos - rawCos * gradW) * invW;
+    float3 gradRawSin = (gradSin - rawSin * gradW) * invW;
+    float rawLength = sqrt(rawCos * rawCos + rawSin * rawSin);
+    float floorLength = 1.0 - normalization;
+
+    CultPhasor result;
+    if (rawLength > floorLength)
+    {
+        float outCos = rawCos / rawLength;
+        float outSin = rawSin / rawLength;
+        float3 gradLength = (rawCos * gradRawCos + rawSin * gradRawSin) / rawLength;
+        result.cos = float4((gradRawCos - outCos * gradLength) / rawLength, outCos);
+        result.sin = float4((gradRawSin - outSin * gradLength) / rawLength, outSin);
+    }
+    else
+    {
+        result.cos = float4(gradRawCos / floorLength, rawCos / floorLength);
+        result.sin = float4(gradRawSin / floorLength, rawSin / floorLength);
+    }
+
+    return result;
+}
+
+#endif

packages/cultmath/shaders/CultMath.Phacelle.hlsl.meta

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

packages/cultmath/shaders/CultMath.hlsl

@@ -503,6 +503,8 @@ CultCellular cultmath_cellular(float3 p)
     return result;
 }
 
+#include "CultMath.Phacelle.hlsl"
+
 float cultmath_value_noise(float2 position)
 {
     float2 cell = floor(position);

packages/cultmath/src/CultMath/CultMath.csproj

@@ -8,6 +8,8 @@
     <RootNamespace>CultMath</RootNamespace>
     <AssemblyName>CultMath</AssemblyName>
     <PackageId>GameCult.Math</PackageId>
+    <!-- MIT overall; math.Phacelle.cs and CultPhasor.cs (and the two Phacelle HLSL files) are MPL-2.0. -->
+    <PackageLicenseExpression>MIT AND MPL-2.0</PackageLicenseExpression>
     <Version>0.2.4</Version>
     <Description>Clean-room C# math shaped for HLSL-style renderer code.</Description>
     <Authors>GameCult</Authors>
@@ -21,6 +23,8 @@
   <ItemGroup>
     <None Include="..\..\README.md" Pack="true" PackagePath="\" />
     <None Include="..\..\shaders\CultMath.hlsl" Pack="true" PackagePath="contentFiles\any\any\shaders\CultMath.hlsl" CopyToOutputDirectory="PreserveNewest" />
+    <None Include="..\..\shaders\CultMath.Phacelle.hlsl" Pack="true" PackagePath="contentFiles\any\any\shaders\CultMath.Phacelle.hlsl" CopyToOutputDirectory="PreserveNewest" />
+    <None Include="..\..\THIRD-PARTY-NOTICES.md" Pack="true" PackagePath="\" />
   </ItemGroup>
 
 </Project>

packages/cultmath/src/CultMath/CultPhasor.cs

@@ -0,0 +1,27 @@
+// This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0. If a copy of
+// the MPL was not distributed with this file, You can obtain one at https://mozilla.org/MPL/2.0/.
+//
+// The result type of math.phacelle (math.Phacelle.cs), Rune Skovbo Johansen's Phacelle noise
+// (Copyright (c) 2025, MPL-2.0; see THIRD-PARTY-NOTICES.md).
+
+namespace CultMath;
+
+/// <summary>
+/// The value-and-gradient result of <see cref="math.phacelle(float3, float3, float, float)"/>
+/// (invariant 8's struct return shape). Both fields use the <c>float4(gradient.xyz, value.w)</c>
+/// layout of a bare value-and-gradient return.
+/// </summary>
+public struct CultPhasor
+{
+    /// <summary>(&#8711;cos, cos): the normalized cosine wave and its gradient.</summary>
+    public float4 cos;
+
+    /// <summary>(&#8711;sin, sin): the normalized sine wave and its gradient.</summary>
+    public float4 sin;
+
+    public CultPhasor(float4 cos, float4 sin)
+    {
+        this.cos = cos;
+        this.sin = sin;
+    }
+}

packages/cultmath/src/CultMath/CultPhasor.cs.meta

@@ -0,0 +1,2 @@
+fileFormatVersion: 2
+guid: 807ae2c8a758cf31896da12f1e5047bd

packages/cultmath/src/CultMath/math.Phacelle.cs

@@ -0,0 +1,127 @@
+// This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0. If a copy of
+// the MPL was not distributed with this file, You can obtain one at https://mozilla.org/MPL/2.0/.
+//
+// Phacelle noise. Copyright (c) 2025 Rune Skovbo Johansen (Shadertoy t3dyWl; MPL-2.0), read from the
+// C# port PhacelleNoise in lpmitchell/AdvancedTerrainErosion, converted to burstable C# by Luke
+// Mitchell (2026). Modified by GameCult: 3D cells, a caller-supplied stripe wave vector, an exact
+// gradient and exact pruning. This file and CultPhasor.cs are the MPL-2.0 files of CultMath; the
+// rest of the package is MIT (see THIRD-PARTY-NOTICES.md).
+
+namespace CultMath;
+
+public static partial class math
+{
+    // Johansen's Phacelle noise (phase + cell), generalised to 3D cells with a caller-supplied
+    // stripe wave vector. Each of the 4x4x4 cells around p carries one jittered feature point
+    // (jitter in [-0.5, 0.5) per axis, from pcg3d over the integer cell coordinate, the same
+    // integer-hashing rule as cellular) and one cosine/sine wave that runs along `side` (the wave
+    // vector: its length is 2*pi times the stripe frequency, and the caller chooses which
+    // perpendicular it is; on a unit sphere that is cross(p_hat, flow)). The waves are blended with
+    // the bell weight max(0, exp(-2 d^2) - 0.01111) (exactly 0 at distance 1.5, the nearest a
+    // feature outside the 4x4x4 block can be in 2D, which is what keeps cell-grid lines out), and
+    // the blended phasor is normalized: its length is stretched to 1 wherever the raw length is at
+    // least 1 - normalization, and scaled by 1 / (1 - normalization) below that. `offset` is the
+    // phase in cycles (1 is a full turn). Value and normalization follow Johansen's
+    // PhacelleNoise as ported in lpmitchell/AdvancedTerrainErosion.
+    //
+    // `normalization` is meant to lie in [0, 1]. At 0 the output is the raw blend (a phasor of
+    // length at most 1); at 1 the floor is 0 and every output is stretched to length 1. Values of
+    // 1 or more saturate at that (the floor is never above the raw length), and a negative value
+    // scales every output down by 1 / (1 - normalization) < 1. The normalized gradient scales as
+    // 1 / |I| where I is the blended phasor (below), so it grows wherever the blend nearly
+    // cancels; a normalization near 1 lets that happen at every point.
+    //
+    // The gradient is exact: it carries the weight derivative (grad w = -4 v exp(-2 |v|^2) where
+    // w > 0) and the normalization, not upstream's convention grad(cos) ~ -sin * side, which
+    // ignores the weight gradient (Asura probe: 16.8 deg off at the median, 97.8 at the 95th
+    // percentile; the exact form is 0.00 and 0.03). With P = sum(w * (cos, sin)) and W = sum(w),
+    // I = P / W has grad I = (grad P - I * grad W) / W; above the threshold the output is I / |I|
+    // with gradient (grad I - out * grad|I|) / |I|, below it grad I / (1 - normalization). The
+    // value is continuous at both seams below; the gradient is not.
+    //  * Where |I| crosses the threshold, the gradient jumps between the two branches.
+    //  * Where a cell's weight leaves its support (d^2 = ln(1 / 0.01111) / 2 = 2.24995, distance
+    //    1.49998), grad w jumps from -4 * 0.01111 * v to 0: a per-cell jump of length
+    //    0.0444 * |v|, 0.0667 at |v| = 1.5. Measured over random points, the resulting jump in
+    //    the output gradient is a median 1.5% of its length and at most about 27%.
+    //
+    // W is the sum of 64 weights and is positive for every input pcg3d produces (the smallest
+    // W found by hill-climbing the hash-driven cell contents is 0.056). It is 0, and the output
+    // NaN, only if all 64 weights vanish, which takes adversarial jitter no hash yields.
+    //
+    // Pruning is exact: per axis, |local - g| - 0.5 (floored at 0) lower-bounds |v| for any jitter,
+    // so a cell whose summed square is at least 2.25 lies outside the support (2.24995 < 2.25) and
+    // has weight exactly 0: it contributes nothing to any sum. Cells whose bound sits just under
+    // 2.25 can still have a nonzero weight, so the cut cannot be lowered. On average the prune
+    // keeps about 45 of the 64 cells (a saving of about 30%); only about 14 have a nonzero weight.
+    // Only the internal overload can switch the prune off; PhacelleTests uses it to pin bit equality.
+    public static CultPhasor phacelle(float3 p, float3 side, float offset, float normalization) =>
+        phacelle(p, side, offset, normalization, true);
+
+    internal static CultPhasor phacelle(float3 p, float3 side, float offset, float normalization, bool prune)
+    {
+        var cell = floor(p);
+        var local = p - cell;
+        var phaseOffset = offset * TAU;
+        var sumW = 0.0f;
+        var sumCos = 0.0f;
+        var sumSin = 0.0f;
+        var gradW = new float3(0.0f, 0.0f, 0.0f);
+        var gradCos = new float3(0.0f, 0.0f, 0.0f);
+        var gradSin = new float3(0.0f, 0.0f, 0.0f);
+
+        for (var dz = -1; dz <= 2; dz++)
+        {
+            var gz = max(0.0f, abs(local.z - (float)dz) - 0.5f);
+            for (var dy = -1; dy <= 2; dy++)
+            {
+                var gy = max(0.0f, abs(local.y - (float)dy) - 0.5f);
+                for (var dx = -1; dx <= 2; dx++)
+                {
+                    var gx = max(0.0f, abs(local.x - (float)dx) - 0.5f);
+                    if (prune && gx * gx + gy * gy + gz * gz >= 2.25f)
+                        continue;
+
+                    var gridStep = new float3(dx, dy, dz);
+                    var hash = pcg3d(int3(cell + gridStep));
+                    var jitter = new float3(cellular_unit(hash.x) - 0.5f, cellular_unit(hash.y) - 0.5f, cellular_unit(hash.z) - 0.5f);
+                    var v = local - gridStep - jitter;
+                    var falloff = exp(-2.0f * dot(v, v));
+                    var w = max(0.0f, falloff - 0.01111f);
+                    var dw = w > 0.0f ? (-4.0f * falloff) * v : new float3(0.0f, 0.0f, 0.0f);
+                    var phase = dot(v, side) + phaseOffset;
+                    var c = cos(phase);
+                    var s = sin(phase);
+
+                    sumW += w;
+                    sumCos += w * c;
+                    sumSin += w * s;
+                    gradW += dw;
+                    gradCos += dw * c - (w * s) * side;
+                    gradSin += dw * s + (w * c) * side;
+                }
+            }
+        }
+
+        var invW = 1.0f / sumW;
+        var rawCos = sumCos * invW;
+        var rawSin = sumSin * invW;
+        var gradRawCos = (gradCos - rawCos * gradW) * invW;
+        var gradRawSin = (gradSin - rawSin * gradW) * invW;
+        var rawLength = sqrt(rawCos * rawCos + rawSin * rawSin);
+        var floorLength = 1.0f - normalization;
+
+        if (rawLength > floorLength)
+        {
+            var outCos = rawCos / rawLength;
+            var outSin = rawSin / rawLength;
+            var gradLength = (rawCos * gradRawCos + rawSin * gradRawSin) / rawLength;
+            return new CultPhasor(
+                new float4((gradRawCos - outCos * gradLength) / rawLength, outCos),
+                new float4((gradRawSin - outSin * gradLength) / rawLength, outSin));
+        }
+
+        return new CultPhasor(
+            new float4(gradRawCos / floorLength, rawCos / floorLength),
+            new float4(gradRawSin / floorLength, rawSin / floorLength));
+    }
+}

packages/cultmath/src/CultMath/math.Phacelle.cs.meta

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

packages/cultmath/tests/CultMath.Tests/HlslMirrorTests.cs

@@ -7,7 +7,7 @@ public sealed class HlslMirrorTests
     [Fact]
     public void HlslMirrorDoesNotPublishGeometryKernels()
     {
-        var source = File.ReadAllText(Path.Combine(GetCultMathRoot(), "shaders", "CultMath.hlsl"));
+        var source = HlslSourceCompatibilityTests.ReadShaderSource(GetCultMathRoot());
         Assert.DoesNotContain("AdvancedErosionFilter.hlsl", source);
         Assert.DoesNotContain("Planetary.hlsl", source);
         Assert.DoesNotContain("spherical_erosion", source);
@@ -17,7 +17,7 @@ public sealed class HlslMirrorTests
     [Fact]
     public void HlslMirrorPublishesCultMathPrimitives()
     {
-        var include = File.ReadAllText(Path.Combine(GetCultMathRoot(), "shaders", "CultMath.hlsl"));
+        var include = HlslSourceCompatibilityTests.ReadShaderSource(GetCultMathRoot());
         var requiredSymbols = new[]
         {
             "cultmath_radians",
@@ -62,6 +62,17 @@ public sealed class HlslMirrorTests
         Assert.DoesNotContain("cultmath_normalize", include);
     }
 
+    [Theory]
+    [InlineData("CultMath.hlsl")]
+    [InlineData("CultMath.Phacelle.hlsl")]
+    public void UnityPackageShaderIsIdenticalToTheCanonicalShader(string name)
+    {
+        var root = GetCultMathRoot();
+        Assert.Equal(
+            File.ReadAllText(Path.Combine(root, "shaders", name)).ReplaceLineEndings("\n"),
+            File.ReadAllText(Path.Combine(root, "unity", "org.gamecult.cultmath", "Shaders", name)).ReplaceLineEndings("\n"));
+    }
+
     private static string GetCultMathRoot()
     {
         var directory = new DirectoryInfo(AppContext.BaseDirectory);

packages/cultmath/tests/CultMath.Tests/HlslSourceCompatibilityTests.cs

@@ -30,6 +30,7 @@ public sealed class HlslSourceCompatibilityTests
         Assert.NotNull(ShaderFunction(shader, "cultmath_snoise_grad", typeof(float3)));
         Assert.NotNull(ShaderFunction(shader, "cultmath_fbm_grad", typeof(float3), typeof(int), typeof(float), typeof(float)));
         Assert.NotNull(ShaderFunction(shader, "cultmath_ridged_grad", typeof(float3), typeof(int), typeof(float), typeof(float)));
+        Assert.NotNull(ShaderFunction(shader, "cultmath_phacelle", typeof(float3), typeof(float3), typeof(float), typeof(float)));
     }
 
     // HLSL functions carry no access modifier, so they compile as private instance methods.
@@ -130,7 +131,63 @@ public sealed class HlslSourceCompatibilityTests
         }
 
         Assert.True(mismatches.Count == 0, $"{mismatches.Count} mismatches:{Environment.NewLine}{string.Join(Environment.NewLine, mismatches.Take(12))}");
-        Assert.Equal(33, compared.Count);
+        Assert.Equal(34, compared.Count);
+    }
+
+    /// <summary>
+    /// The generic comparison above passes about a hundred inputs through <c>cultmath_phacelle</c>, and
+    /// its prune only changes an output for the rare cell whose weight is nonzero but whose lower
+    /// bound sits just under the 2.25 cut (Soul's F6 lesson for cellular): a mirror whose cut is 2.0
+    /// survives a hundred inputs. A dense sweep of realistic arguments closes that gap.
+    /// </summary>
+    [Fact]
+    public void PhacelleMirrorMatchesCSharpBitForBitAcrossADenseSweep()
+    {
+        var (assembly, errors) = CompileShaderMirror();
+        Assert.True(assembly is not null, string.Join(Environment.NewLine, errors));
+        var shaderType = assembly!.GetType("CultMathHlsl.HlslShader")!;
+        var shader = Activator.CreateInstance(shaderType);
+        var mirror = ShaderFunction(shaderType, "cultmath_phacelle", typeof(float3), typeof(float3), typeof(float), typeof(float))!;
+
+        var random = new System.Random(0x5EED2);
+        float Next(float extent) => (random.NextSingle() * 2.0f - 1.0f) * extent;
+        for (var i = 0; i < 20000; i++)
+        {
+            var p = new float3(Next(50.0f), Next(50.0f), Next(50.0f));
+            var side = new float3(Next(9.0f), Next(9.0f), Next(9.0f));
+            var offset = random.NextSingle();
+            var normalization = random.NextSingle();
+            var expected = Values(phacelle(p, side, offset, normalization));
+            var actual = Values(mirror.Invoke(shader, new object[] { p, side, offset, normalization })!);
+            Assert.True(BitwiseEqual(expected, actual), $"mirror differs from C# at p = {p}, side = {side}");
+        }
+    }
+
+    /// <summary>
+    /// The dense sweep almost never meets a cell whose box bound is just under the 2.25 cut and whose
+    /// weight is still nonzero: a mirror cut of 2.245, 2.24 or 2.2 survives it. These points sit next
+    /// to such cells (<see cref="PhacelleTests.PruneCutPoints"/>), so the mirror must match C# bit for
+    /// bit, and must match the unpruned C# sum (the shader has no prune switch, so that is the
+    /// independent statement that the cut is exact).
+    /// </summary>
+    [Fact]
+    public void PhacelleMirrorMatchesCSharpBitForBitNextToCellsJustInsideTheCut()
+    {
+        var (assembly, errors) = CompileShaderMirror();
+        Assert.True(assembly is not null, string.Join(Environment.NewLine, errors));
+        var shaderType = assembly!.GetType("CultMathHlsl.HlslShader")!;
+        var shader = Activator.CreateInstance(shaderType);
+        var mirror = ShaderFunction(shaderType, "cultmath_phacelle", typeof(float3), typeof(float3), typeof(float), typeof(float))!;
+
+        foreach (var p in PhacelleTests.PruneCutPoints)
+        {
+            Assert.True(PhacelleTests.HasCellJustInsideTheCut(p), $"no cell just inside the cut at {p}");
+            var actual = Values(mirror.Invoke(shader, new object[] { p, PhacelleTests.PruneCutSide, PhacelleTests.PruneCutOffset, PhacelleTests.PruneCutNormalization })!);
+            var pruned = Values(phacelle(p, PhacelleTests.PruneCutSide, PhacelleTests.PruneCutOffset, PhacelleTests.PruneCutNormalization));
+            var unpruned = Values(phacelle(p, PhacelleTests.PruneCutSide, PhacelleTests.PruneCutOffset, PhacelleTests.PruneCutNormalization, false));
+            Assert.True(BitwiseEqual(pruned, actual), $"mirror differs from C# at p = {p}");
+            Assert.True(BitwiseEqual(unpruned, actual), $"mirror differs from the unpruned C# sum at p = {p}");
+        }
     }
 
     /// <summary>
@@ -253,6 +310,18 @@ public sealed class HlslSourceCompatibilityTests
                 ? (float.IsNaN(a) && float.IsNaN(b)) || BitConverter.SingleToInt32Bits(a) == BitConverter.SingleToInt32Bits(b)
                 : p.First.Value.Equals(p.Second.Value)));
 
+    /// <summary>
+    /// The text of <c>CultMath.hlsl</c> with each <c>#include "name"</c> replaced in place by the
+    /// included file from the same directory, as the shader compiler would resolve it. Only the
+    /// package's own includes exist, so a missing file is an error rather than a skip.
+    /// </summary>
+    internal static string ReadShaderSource(string cultMathRoot)
+    {
+        var shaders = Path.Combine(cultMathRoot, "shaders");
+        return Regex.Replace(File.ReadAllText(Path.Combine(shaders, "CultMath.hlsl")), @"(?m)^#include ""([^""]+)""[ \t]*\r?$",
+            match => File.ReadAllText(Path.Combine(shaders, match.Groups[1].Value)));
+    }
+
     /// <summary>The documented HLSL-to-C# transformations, and nothing else.</summary>
     internal static string TransformHlslToCSharp(string hlsl)
     {
@@ -281,7 +350,7 @@ public sealed class HlslSourceCompatibilityTests
 
     private static (Assembly? Assembly, IReadOnlyList<string> Errors) CompileShaderMirror()
     {
-        var hlsl = File.ReadAllText(Path.Combine(FindCultMathRoot(), "shaders", "CultMath.hlsl"));
+        var hlsl = ReadShaderSource(FindCultMathRoot());
         var tree = CSharpSyntaxTree.ParseText(TransformHlslToCSharp(hlsl));
         var references = ((string)AppContext.GetData("TRUSTED_PLATFORM_ASSEMBLIES")!)
             .Split(Path.PathSeparator)

packages/cultmath/tests/CultMath.Tests/PhacelleTests.cs

@@ -0,0 +1,340 @@
+using System;
+using CultMath;
+using Xunit;
+
+namespace CultMath.Tests;
+
+public sealed class PhacelleTests
+{
+    // Every tolerance below is at least 5x the largest error the same test measures at its seed.
+    private const float GradientEps = 1.0e-3f;
+    private const float GradientTolerance = 5.0e-3f;
+    private const double ValueTolerance = 1.6e-5;
+    private const double ThresholdBand = 0.02;
+
+    private static float3 RandomPoint(System.Random random, float extent = 3.0f) => new(
+        (random.NextSingle() * 2.0f - 1.0f) * extent, (random.NextSingle() * 2.0f - 1.0f) * extent, (random.NextSingle() * 2.0f - 1.0f) * extent);
+
+    // A wave vector of length 2*pi times a stripe frequency in [0.5, 1.5], in a random direction.
+    private static float3 RandomSide(System.Random random)
+    {
+        float3 direction;
+        do
+        {
+            direction = new float3(random.NextSingle() - 0.5f, random.NextSingle() - 0.5f, random.NextSingle() - 0.5f);
+        }
+        while (math.length(direction) < 0.1f);
+
+        return math.normalize(direction) * (math.TAU * (0.5f + random.NextSingle()));
+    }
+
+    // ---- An independent double-precision evaluation: all 64 cells, no prune, no gradient. ----
+    //
+    // It shares only the hash (pcg3d over the integer cell coordinate, top 24 bits) and the constants
+    // the function is specified by: weight max(0, exp(-2 d^2) - 0.01111), jitter in [-0.5, 0.5).
+
+    private readonly record struct Reference(double Cos, double Sin, double RawLength, double SupportEdgeGap);
+
+    private static Reference Evaluate(float3 p, float3 side, float offset, float normalization)
+    {
+        double fx = Math.Floor(p.x), fy = Math.Floor(p.y), fz = Math.Floor(p.z);
+        double lx = p.x - fx, ly = p.y - fy, lz = p.z - fz;
+        double sumW = 0.0, sumCos = 0.0, sumSin = 0.0, edgeGap = double.MaxValue;
+        var supportRadius = Math.Sqrt(Math.Log(1.0 / 0.01111) / 2.0);
+        for (var dz = -1; dz <= 2; dz++)
+        for (var dy = -1; dy <= 2; dy++)
+        for (var dx = -1; dx <= 2; dx++)
+        {
+            var hash = math.pcg3d(new int3((int)fx + dx, (int)fy + dy, (int)fz + dz));
+            double vx = lx - dx - (math.cellular_unit(hash.x) - 0.5);
+            double vy = ly - dy - (math.cellular_unit(hash.y) - 0.5);
+            double vz = lz - dz - (math.cellular_unit(hash.z) - 0.5);
+            var d2 = vx * vx + vy * vy + vz * vz;
+            edgeGap = Math.Min(edgeGap, Math.Abs(Math.Sqrt(d2) - supportRadius));
+            var w = Math.Max(0.0, Math.Exp(-2.0 * d2) - 0.01111);
+            var phase = vx * side.x + vy * side.y + vz * side.z + (double)offset * math.TAU;
+            sumW += w;
+            sumCos += w * Math.Cos(phase);
+            sumSin += w * Math.Sin(phase);
+        }
+
+        double rawCos = sumCos / sumW, rawSin = sumSin / sumW;
+        var raw = Math.Sqrt(rawCos * rawCos + rawSin * rawSin);
+        var length = Math.Max(1.0 - normalization, raw);
+        return new Reference(rawCos / length, rawSin / length, raw, edgeGap);
+    }
+
+    private static float3 CentralDifference(Func<float3, float> f, float3 p, float eps) => new(
+        (f(p + new float3(eps, 0.0f, 0.0f)) - f(p - new float3(eps, 0.0f, 0.0f))) / (2.0f * eps),
+        (f(p + new float3(0.0f, eps, 0.0f)) - f(p - new float3(0.0f, eps, 0.0f))) / (2.0f * eps),
+        (f(p + new float3(0.0f, 0.0f, eps)) - f(p - new float3(0.0f, 0.0f, eps))) / (2.0f * eps));
+
+    private static float MaxComponentError(float3 numeric, float4 analytic) => MathF.Max(
+        MathF.Abs(numeric.x - analytic.x), MathF.Max(MathF.Abs(numeric.y - analytic.y), MathF.Abs(numeric.z - analytic.z)));
+
+    // ---- Value: the function is Johansen's, on 3D cells ----
+
+    [Theory]
+    [InlineData(0.0f)]
+    [InlineData(0.5f)]
+    [InlineData(0.9f)]
+    [InlineData(1.0f)]
+    public void ValueMatchesTheIndependentReferenceAndNormalizesAboveTheThreshold(float normalization)
+    {
+        var random = new System.Random(0x7ACE1 + (int)(normalization * 10.0f));
+        var maxError = 0.0;
+        var above = 0;
+        var below = 0;
+        for (var i = 0; i < 3000; i++)
+        {
+            var p = RandomPoint(random, 40.0f);
+            var side = RandomSide(random);
+            var offset = random.NextSingle();
+            var actual = math.phacelle(p, side, offset, normalization);
+            var expected = Evaluate(p, side, offset, normalization);
+
+            maxError = Math.Max(maxError, Math.Max(Math.Abs(actual.cos.w - expected.Cos), Math.Abs(actual.sin.w - expected.Sin)));
+            if (expected.RawLength >= 1.0 - normalization + ThresholdBand)
+            {
+                // The output magnitude is 1 wherever the raw magnitude reaches 1 - normalization.
+                above++;
+                var magnitude = Math.Sqrt((double)actual.cos.w * actual.cos.w + (double)actual.sin.w * actual.sin.w);
+                Assert.True(Math.Abs(magnitude - 1.0) < ValueTolerance, $"magnitude {magnitude} above the threshold at {p}");
+            }
+            else if (expected.RawLength <= 1.0 - normalization - ThresholdBand)
+            {
+                below++;
+            }
+        }
+
+        Assert.True(maxError < ValueTolerance, $"max value error against the reference {maxError}");
+        if (normalization > 0.0f)
+            Assert.True(above >= 500, $"only {above} points above the threshold");
+        if (normalization is > 0.0f and < 1.0f)
+            Assert.True(below >= 50, $"only {below} points below the threshold");
+    }
+
+    [Fact]
+    public void ValueIsPeriodicInOffsetAndAQuarterCycleRotatesCosIntoSin()
+    {
+        // offset is in cycles: a quarter turn maps cos to -sin and sin to cos, and a whole turn is the
+        // identity. The weights and the normalization do not depend on the phase, so this is exact
+        // up to float rounding.
+        var random = new System.Random(0x7ACE2);
+        for (var i = 0; i < 200; i++)
+        {
+            var p = RandomPoint(random);
+            var side = RandomSide(random);
+            var offset = random.NextSingle();
+            var basePhasor = math.phacelle(p, side, offset, 0.5f);
+            var quarter = math.phacelle(p, side, offset + 0.25f, 0.5f);
+            var whole = math.phacelle(p, side, offset + 1.0f, 0.5f);
+
+            Assert.True(MathF.Abs(quarter.cos.w + basePhasor.sin.w) < 2.0e-5f);
+            Assert.True(MathF.Abs(quarter.sin.w - basePhasor.cos.w) < 2.0e-5f);
+            Assert.True(MathF.Abs(whole.cos.w - basePhasor.cos.w) < 2.0e-5f);
+            Assert.True(MathF.Abs(whole.sin.w - basePhasor.sin.w) < 2.0e-5f);
+        }
+    }
+
+    // ---- Invariant 8: the gradient is exact ----
+
+    [Fact]
+    public void GradientOfCosAndSinMatchesCentralDifferencesOnBothSidesOfTheNormalizationThreshold()
+    {
+        // Above the threshold the output is I / |I|; below it, I / (1 - normalization). The gradient
+        // is discontinuous where |I| crosses the threshold, and at the edge of each cell's weight support
+        // (distance 1.49998, where max(0, ...) starts: a jump of 4 d 0.01111 in grad w), so a band
+        // around each is excluded (by the
+        // reference's raw length, since the output alone cannot say how close it is).
+        const float normalization = 0.5f;
+        var random = new System.Random(0x7ACE3);
+        var above = 0;
+        var below = 0;
+        var maxError = 0.0f;
+        for (var attempts = 0; attempts < 20000 && (above < 150 || below < 150); attempts++)
+        {
+            var p = RandomPoint(random);
+            var side = RandomSide(random);
+            var offset = random.NextSingle();
+            var reference = Evaluate(p, side, offset, normalization);
+            var distance = reference.RawLength - (1.0 - normalization);
+            if (Math.Abs(distance) < ThresholdBand || reference.SupportEdgeGap < 3.0 * GradientEps)
+                continue;
+            if (distance > 0.0 ? above >= 150 : below >= 150)
+                continue;
+            if (distance > 0.0) above++; else below++;
+
+            var analytic = math.phacelle(p, side, offset, normalization);
+            var numericCos = CentralDifference(q => math.phacelle(q, side, offset, normalization).cos.w, p, GradientEps);
+            var numericSin = CentralDifference(q => math.phacelle(q, side, offset, normalization).sin.w, p, GradientEps);
+            maxError = MathF.Max(maxError, MathF.Max(MaxComponentError(numericCos, analytic.cos), MaxComponentError(numericSin, analytic.sin)));
+        }
+
+        Assert.True(above >= 150 && below >= 150, $"{above} points above and {below} below the threshold");
+        Assert.True(maxError < GradientTolerance, $"max gradient error against central differences {maxError}");
+    }
+
+    // ---- Stripes run along the flow ----
+
+    [Fact]
+    public void OnTheUnitSphereStripesRunAlongTheFlow()
+    {
+        // side = cross(p_hat, flow) * 2*pi is tangent to the sphere and perpendicular to the flow, so the
+        // cosine should change across the stripes (along side) and hardly at all along them (along
+        // flow). Asura's probe measured 0.212 for this ratio, against 0.202 for the planar upstream.
+        const float radiusInCells = 10.0f;
+        var random = new System.Random(0x7ACE4);
+        var reference = math.normalize(new float3(0.3f, 0.5f, 0.8f));
+        var alongFlow = 0.0;
+        var acrossFlow = 0.0;
+        var used = 0;
+        while (used < 2000)
+        {
+            var hat = math.normalize(new float3(random.NextSingle() - 0.5f, random.NextSingle() - 0.5f, random.NextSingle() - 0.5f));
+            var tangent = reference - hat * math.dot(hat, reference);
+            if (math.length(tangent) < 0.2f)
+                continue;
+
+            var flow = math.normalize(tangent);
+            var sideDirection = math.cross(hat, flow);
+            var phasor = math.phacelle(hat * radiusInCells, sideDirection * math.TAU, 0.25f, 0.5f);
+            var gradient = new float3(phasor.cos.x, phasor.cos.y, phasor.cos.z);
+            alongFlow += Math.Abs(math.dot(gradient, flow));
+            acrossFlow += Math.Abs(math.dot(gradient, sideDirection));
+            used++;
+        }
+
+        var ratio = alongFlow / acrossFlow;
+        Assert.True(ratio <= 0.25, $"mean|d cos/d flow| / mean|d cos/d side| = {ratio}");
+    }
+
+    // ---- Continuity ----
+
+    [Fact]
+    public void ValueStepDifferenceShrinksTenfoldPerTenfoldStepIncludingAcrossCellPlanes()
+    {
+        // Pairs straddle an integer plane on each axis (where the 4x4x4 window moves by one cell and a
+        // weight that did not vanish at the window edge would jump) and also start at random points.
+        // Lipschitz continuity means the largest difference over a step scales with the step.
+        var random = new System.Random(0x7ACE5);
+        var steps = new[] { 1.0e-3f, 1.0e-4f, 1.0e-5f };
+        var maxDifference = new float[steps.Length];
+        for (var i = 0; i < 3000; i++)
+        {
+            var p = RandomPoint(random);
+            var side = RandomSide(random);
+            var offset = random.NextSingle();
+            var axis = i % 3;
+            if (i % 2 == 0)
+                p = SetAxis(p, axis, MathF.Round(GetAxis(p, axis)));
+            for (var k = 0; k < steps.Length; k++)
+            {
+                var low = math.phacelle(SetAxis(p, axis, GetAxis(p, axis) - steps[k] * 0.5f), side, offset, 0.5f);
+                var high = math.phacelle(SetAxis(p, axis, GetAxis(p, axis) + steps[k] * 0.5f), side, offset, 0.5f);
+                maxDifference[k] = MathF.Max(maxDifference[k], MathF.Max(MathF.Abs(low.cos.w - high.cos.w), MathF.Abs(low.sin.w - high.sin.w)));
+            }
+        }
+
+        Assert.True(maxDifference[1] * 8.0f <= maxDifference[0], $"1e-3: {maxDifference[0]}, 1e-4: {maxDifference[1]}");
+        Assert.True(maxDifference[2] * 8.0f <= maxDifference[1], $"1e-4: {maxDifference[1]}, 1e-5: {maxDifference[2]}");
+    }
+
+    private static float GetAxis(float3 v, int axis) => axis == 0 ? v.x : axis == 1 ? v.y : v.z;
+
+    private static float3 SetAxis(float3 v, int axis, float value) =>
+        axis == 0 ? new float3(value, v.y, v.z) : axis == 1 ? new float3(v.x, value, v.z) : new float3(v.x, v.y, value);
+
+    // ---- Exact pruning ----
+
+    [Fact]
+    public void PrunedOutputIsBitEqualToTheUnprunedSum()
+    {
+        var random = new System.Random(0x7ACE6);
+        for (var i = 0; i < 20000; i++)
+        {
+            var p = RandomPoint(random, 50.0f);
+            var side = RandomSide(random);
+            var offset = random.NextSingle();
+            var normalization = random.NextSingle();
+            var pruned = math.phacelle(p, side, offset, normalization);
+            var unpruned = math.phacelle(p, side, offset, normalization, false);
+
+            AssertBitEqual(unpruned.cos, pruned.cos);
+            AssertBitEqual(unpruned.sin, pruned.sin);
+        }
+    }
+
+    // Points next to a cell whose jitter sits at the far corner of its cube: the cell's box bound is
+    // just under the 2.25 cut (2.2451 to 2.25) while its real distance is still under the support
+    // radius (d^2 < 2.24995), so its weight is small but nonzero and a prune that cuts even slightly
+    // low (2.245, 2.24, 2.2) drops a contribution and changes output bits. Random points meet such
+    // a cell about once in 10^4 evaluations, which is why the random sweep above let those cuts
+    // through. Each was found once by an offline search (cells K in [-12, 12]^3 whose jitter on one
+    // axis lies within 0.0008 of the cube face; p is the feature moved 1.4996 out along that axis,
+    // so |p - feature| = 1.4996 and the bound is about 1.4988). The witness check below re-derives
+    // that property from the hash, so a change to pcg3d fails loudly instead of degenerating.
+    // Shared with the HLSL mirror test.
+    internal static readonly float3[] PruneCutPoints =
+    {
+        new(1.77850008f, -2.16020727f, -0.9991166f),
+        new(1.09346282f, -2.37882328f, -3.00082445f),
+        new(-6.39144945f, -3.99953747f, -0.610023975f),
+        new(-0.99940908f, 4.16302729f, 4.07277536f),
+        new(3.99912024f, 7.21890545f, 1.09053993f),
+        new(5.43959713f, 1.99937367f, -6.14929914f),
+    };
+
+    internal static readonly float3 PruneCutSide = new(2.1f, -3.7f, 1.3f);
+    internal const float PruneCutOffset = 0.3f;
+    internal const float PruneCutNormalization = 0.5f;
+
+    // True when some cell of the 4x4x4 window has a box bound in [2.2451, 2.25) and a positive weight.
+    internal static bool HasCellJustInsideTheCut(float3 p)
+    {
+        double fx = Math.Floor(p.x), fy = Math.Floor(p.y), fz = Math.Floor(p.z);
+        double lx = p.x - fx, ly = p.y - fy, lz = p.z - fz;
+        for (var dz = -1; dz <= 2; dz++)
+        for (var dy = -1; dy <= 2; dy++)
+        for (var dx = -1; dx <= 2; dx++)
+        {
+            var gx = Math.Max(0.0, Math.Abs(lx - dx) - 0.5);
+            var gy = Math.Max(0.0, Math.Abs(ly - dy) - 0.5);
+            var gz = Math.Max(0.0, Math.Abs(lz - dz) - 0.5);
+            var bound = gx * gx + gy * gy + gz * gz;
+            if (bound < 2.2451 || bound >= 2.25)
+                continue;
+
+            var hash = math.pcg3d(new int3((int)fx + dx, (int)fy + dy, (int)fz + dz));
+            double vx = lx - dx - (math.cellular_unit(hash.x) - 0.5);
+            double vy = ly - dy - (math.cellular_unit(hash.y) - 0.5);
+            double vz = lz - dz - (math.cellular_unit(hash.z) - 0.5);
+            if (Math.Exp(-2.0 * (vx * vx + vy * vy + vz * vz)) > 0.01111)
+                return true;
+        }
+
+        return false;
+    }
+
+    [Fact]
+    public void PrunedOutputIsBitEqualToTheUnprunedSumNextToCellsJustInsideTheCut()
+    {
+        foreach (var p in PruneCutPoints)
+        {
+            Assert.True(HasCellJustInsideTheCut(p), $"no cell just inside the cut at {p}: the fixture no longer exercises the prune");
+            var pruned = math.phacelle(p, PruneCutSide, PruneCutOffset, PruneCutNormalization);
+            var unpruned = math.phacelle(p, PruneCutSide, PruneCutOffset, PruneCutNormalization, false);
+
+            AssertBitEqual(unpruned.cos, pruned.cos);
+            AssertBitEqual(unpruned.sin, pruned.sin);
+        }
+    }
+
+    private static void AssertBitEqual(float4 expected, float4 actual)
+    {
+        Assert.Equal(BitConverter.SingleToInt32Bits(expected.x), BitConverter.SingleToInt32Bits(actual.x));
+        Assert.Equal(BitConverter.SingleToInt32Bits(expected.y), BitConverter.SingleToInt32Bits(actual.y));
+        Assert.Equal(BitConverter.SingleToInt32Bits(expected.z), BitConverter.SingleToInt32Bits(actual.z));
+        Assert.Equal(BitConverter.SingleToInt32Bits(expected.w), BitConverter.SingleToInt32Bits(actual.w));
+    }
+}

packages/cultmath/unity/org.gamecult.cultmath/LICENSE.md

@@ -4,7 +4,7 @@ CultMath is licensed under the Mozilla Public License 2.0 (`MPL-2.0`).
 
 The complete license text is available from the Mozilla Foundation at
 <https://www.mozilla.org/MPL/2.0/>. Source and modification history for the
-precompiled assembly in this package are available in the GameCult/CultMath
-repository at the corresponding package tag.
+precompiled assembly in this package are available in the GameCult/CultLib
+repository (`packages/cultmath`) at the corresponding package tag.
 
 Third-party code retains its original license. See `THIRD-PARTY-NOTICES.md`.

packages/cultmath/unity/org.gamecult.cultmath/README.md

@@ -2,7 +2,7 @@
 
 This is CultMath's Unity 2021.3-compatible Git-UPM surface. It contains the
 precompiled `CultMath.dll`, portable symbols, the numeric
-`Shaders/CultMath.hlsl` mirror, and the `CultMath.UnityBridge` assembly.
+`Shaders/CultMath.hlsl` mirror (with its `CultMath.Phacelle.hlsl` include), and the `CultMath.UnityBridge` assembly.
 
 `CultMath.UnityBridge` is the only C# source Unity compiles here. It holds the
 UnityEngine conversions the engine-free core cannot declare: `ToCultMath()` on

packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Phacelle.hlsl

@@ -0,0 +1,99 @@
+#ifndef CULTMATH_PHACELLE_HLSL
+#define CULTMATH_PHACELLE_HLSL
+
+// This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0. If a copy of
+// the MPL was not distributed with this file, You can obtain one at https://mozilla.org/MPL/2.0/.
+//
+// Phacelle noise. Copyright (c) 2025 Rune Skovbo Johansen (Shadertoy t3dyWl; MPL-2.0), read from the
+// C# port PhacelleNoise in lpmitchell/AdvancedTerrainErosion, converted to burstable C# by Luke
+// Mitchell (2026). Modified by GameCult: 3D cells, a caller-supplied stripe wave vector, an exact
+// gradient and exact pruning. This file is an MPL-2.0 file; CultMath.hlsl, which includes it after
+// its own definitions of CULTMATH_TAU and cultmath_pcg3d, is MIT (see THIRD-PARTY-NOTICES.md).
+
+// Invariant 8's struct return shape for cultmath_phacelle, compared field by field and bit for bit
+// against CultMath.CultPhasor, same field order.
+struct CultPhasor
+{
+    float4 cos;
+    float4 sin;
+};
+
+// Johansen's Phacelle noise on 3D cells. See math.phacelle's comment (math.Phacelle.cs) for the stripe wave
+// vector, the weight, the normalization, the exact gradient and why the prune is exact. Mirrors
+// the C# body operation for operation; the C#-only internal overload that switches the prune off
+// has no shader twin.
+CultPhasor cultmath_phacelle(float3 p, float3 side, float offset, float normalization)
+{
+    float3 cell = floor(p);
+    float3 local = p - cell;
+    float phaseOffset = offset * CULTMATH_TAU;
+    float sumW = 0.0;
+    float sumCos = 0.0;
+    float sumSin = 0.0;
+    float3 gradW = float3(0.0, 0.0, 0.0);
+    float3 gradCos = float3(0.0, 0.0, 0.0);
+    float3 gradSin = float3(0.0, 0.0, 0.0);
+
+    for (int dz = -1; dz <= 2; dz++)
+    {
+        float gz = max(0.0, abs(local.z - (float)dz) - 0.5);
+        for (int dy = -1; dy <= 2; dy++)
+        {
+            float gy = max(0.0, abs(local.y - (float)dy) - 0.5);
+            for (int dx = -1; dx <= 2; dx++)
+            {
+                float gx = max(0.0, abs(local.x - (float)dx) - 0.5);
+                if (gx * gx + gy * gy + gz * gz >= 2.25)
+                    continue;
+
+                float3 gridStep = float3(dx, dy, dz);
+                int3 hash = cultmath_pcg3d(int3(cell + gridStep));
+                float3 jitter = float3(
+                    ((uint)hash.x >> 8) * (1.0 / 16777216.0) - 0.5,
+                    ((uint)hash.y >> 8) * (1.0 / 16777216.0) - 0.5,
+                    ((uint)hash.z >> 8) * (1.0 / 16777216.0) - 0.5);
+                float3 v = local - gridStep - jitter;
+                float falloff = exp(-2.0 * dot(v, v));
+                float w = max(0.0, falloff - 0.01111);
+                float3 dw = w > 0.0 ? (-4.0 * falloff) * v : float3(0.0, 0.0, 0.0);
+                float phase = dot(v, side) + phaseOffset;
+                float c = cos(phase);
+                float s = sin(phase);
+
+                sumW += w;
+                sumCos += w * c;
+                sumSin += w * s;
+                gradW += dw;
+                gradCos += dw * c - (w * s) * side;
+                gradSin += dw * s + (w * c) * side;
+            }
+        }
+    }
+
+    float invW = 1.0 / sumW;
+    float rawCos = sumCos * invW;
+    float rawSin = sumSin * invW;
+    float3 gradRawCos = (gradCos - rawCos * gradW) * invW;
+    float3 gradRawSin = (gradSin - rawSin * gradW) * invW;
+    float rawLength = sqrt(rawCos * rawCos + rawSin * rawSin);
+    float floorLength = 1.0 - normalization;
+
+    CultPhasor result;
+    if (rawLength > floorLength)
+    {
+        float outCos = rawCos / rawLength;
+        float outSin = rawSin / rawLength;
+        float3 gradLength = (rawCos * gradRawCos + rawSin * gradRawSin) / rawLength;
+        result.cos = float4((gradRawCos - outCos * gradLength) / rawLength, outCos);
+        result.sin = float4((gradRawSin - outSin * gradLength) / rawLength, outSin);
+    }
+    else
+    {
+        result.cos = float4(gradRawCos / floorLength, rawCos / floorLength);
+        result.sin = float4(gradRawSin / floorLength, rawSin / floorLength);
+    }
+
+    return result;
+}
+
+#endif

packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.Phacelle.hlsl.meta

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

packages/cultmath/unity/org.gamecult.cultmath/Shaders/CultMath.hlsl

@@ -503,6 +503,8 @@ CultCellular cultmath_cellular(float3 p)
     return result;
 }
 
+#include "CultMath.Phacelle.hlsl"
+
 float cultmath_value_noise(float2 position)
 {
     float2 cell = floor(position);

packages/cultmath/unity/org.gamecult.cultmath/THIRD-PARTY-NOTICES.md

@@ -19,3 +19,23 @@ permutation and scale as `snoise(float3)`.
 - Copyright 2011 Ashima Arts.
 - Upstream: <https://github.com/ashima/webgl-noise>
 - License: MIT.
+
+## Phacelle Noise
+
+`math.phacelle` and `cultmath_phacelle` follow Rune Skovbo Johansen's Phacelle
+noise (Shadertoy `t3dyWl`; erosion filter `wXcfWn`; blog, March 2026), read
+from the C# port `PhacelleNoise` in `lpmitchell/AdvancedTerrainErosion`
+(`AdvancedTerrainErosion.cs`, converted to burstable C# by Luke Mitchell,
+2026). The cell weight, the phase blend and the normalization follow it. The
+3D-cell generalisation, the caller-supplied stripe wave vector, the exact
+gradient and the exact pruning are new here.
+
+- Copyright (c) 2025 Rune Skovbo Johansen.
+- Upstream: <https://github.com/lpmitchell/AdvancedTerrainErosion>
+- License: MPL-2.0 (<https://mozilla.org/MPL/2.0/>). MPL-2.0 is file-level, so it
+  covers exactly the files that carry the port and nothing else:
+  `src/CultMath/math.Phacelle.cs`, `src/CultMath/CultPhasor.cs`,
+  `shaders/CultMath.Phacelle.hlsl` and the Unity copy of that include. Each
+  carries the MPL-2.0 header and this attribution. The rest of CultMath is MIT,
+  so the `GameCult.Math` NuGet package declares `MIT AND MPL-2.0` and ships this
+  file. The Unity package as a whole is MPL-2.0.

src/GameCult.Caching.MessagePack/CultMathJson.cs

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

src/GameCult.Caching.MessagePack/CultMathResolver.cs

@@ -180,7 +180,7 @@ public sealed class CultMathResolver : IFormatterResolver
                     return value;
                 }),
 
-            // Compositions: matrices are rows, rect is [min, max], CultCellular is [nearest, edge, id].
+            // Compositions: matrices are rows, rect is [min, max], CultCellular is [nearest, edge, id], CultPhasor is [cos, sin].
             Shape<float2x2>((ref MessagePackWriter w, float2x2 v) => { w.WriteArrayHeader(2); PutFloat2(ref w, v[0]); PutFloat2(ref w, v[1]); },
                 (ref MessagePackReader r) =>
                 {
@@ -221,6 +221,21 @@ public sealed class CultMathResolver : IFormatterResolver
                     }
                     return value;
                 }),
+            Shape<CultPhasor>((ref MessagePackWriter w, CultPhasor v) => { w.WriteArrayHeader(2); PutFloat4(ref w, v.cos); PutFloat4(ref w, v.sin); },
+                (ref MessagePackReader r) =>
+                {
+                    var value = default(CultPhasor);
+                    for (int i = 0, n = r.ReadArrayHeader(); i < n; i++)
+                    {
+                        switch (i)
+                        {
+                            case 0: value.cos = GetFloat4(ref r); break;
+                            case 1: value.sin = GetFloat4(ref r); break;
+                            default: r.Skip(); break;
+                        }
+                    }
+                    return value;
+                }),
         };
 
         var map = new Dictionary<Type, object>();

src/GameCult.Caching.MessagePack/README.md

@@ -93,6 +93,7 @@ registered ahead of these wins.
 | `float2x2`, `float3x3` | array of row vectors (row-major) | array of row objects |
 | `rect` | `[min, max]` of `float2` | `{"min":{..},"max":{..}}` |
 | `CultCellular` | `[nearest, edge, id]` | `{"nearest":{..},"edge":{..},"id":..}` |
+| `CultPhasor` | `[cos, sin]` of `float4` | `{"cos":{..},"sin":{..}}` |
 
 Floats are MessagePack float32 (`0xCA`), doubles float64 (`0xCB`), integers the
 shortest MessagePack integer.

src/GameCult.Geometry/CultGeometryDocuments.cs

@@ -153,8 +153,8 @@ namespace GameCult.Geometry
         public string RequestId { get; set; } = string.Empty;
 
         [Key(1)]
-        [CultReference(typeof(CultGeometryDomainDocument))]
-        public string DomainKey { get; set; } = string.Empty;
+        [MessagePackFormatter(typeof(CultRecordRefFormatter<CultGeometryDomainDocument>))]
+        public CultRecordRef<CultGeometryDomainDocument> DomainKey { get; set; }
 
         [Key(2)]
         [CultIndex]
@@ -200,7 +200,7 @@ namespace GameCult.Geometry
         {
             if (request == null) throw new ArgumentNullException(nameof(request));
             return new CultRecordKey("geometry:request:" + StableHash(
-                request.DomainKey,
+                request.DomainKey.Key.Value,
                 request.WorkerGroup,
                 StableVector(request.CameraPosition),
                 StableVector(request.FrustumMin),
@@ -228,8 +228,8 @@ namespace GameCult.Geometry
         public string CutId { get; set; } = string.Empty;
 
         [Key(1)]
-        [CultReference(typeof(CultGeometryBuildRequest))]
-        public string RequestKey { get; set; } = string.Empty;
+        [MessagePackFormatter(typeof(CultRecordRefFormatter<CultGeometryBuildRequest>))]
+        public CultRecordRef<CultGeometryBuildRequest> RequestKey { get; set; }
 
         [Key(2)]
         public string[] SelectedNodes { get; set; } = Array.Empty<string>();
@@ -247,7 +247,7 @@ namespace GameCult.Geometry
         {
             if (cut == null) throw new ArgumentNullException(nameof(cut));
             return new CultRecordKey("geometry:cut:" + StableHash(
-                cut.RequestKey,
+                cut.RequestKey.Key.Value,
                 cut.CutId,
                 StableArray(cut.SelectedNodes),
                 StableArray(cut.DeferredChildRequests),
@@ -313,8 +313,8 @@ namespace GameCult.Geometry
         public string ChunkId { get; set; } = string.Empty;
 
         [Key(1)]
-        [CultReference(typeof(CultGeometrySelectedCutManifest))]
-        public string CutKey { get; set; } = string.Empty;
+        [MessagePackFormatter(typeof(CultRecordRefFormatter<CultGeometrySelectedCutManifest>))]
+        public CultRecordRef<CultGeometrySelectedCutManifest> CutKey { get; set; }
 
         [Key(2)]
         [CultIndex]
@@ -357,7 +357,7 @@ namespace GameCult.Geometry
         {
             if (chunk == null) throw new ArgumentNullException(nameof(chunk));
             return new CultRecordKey("geometry:chunk:" + StableHash(
-                chunk.CutKey,
+                chunk.CutKey.Key.Value,
                 chunk.ChunkId,
                 chunk.SelectedCutId,
                 StableArray(chunk.SourceDomainKeys),

src/GameCult.Geometry/GameCult.Geometry.csproj

@@ -10,6 +10,7 @@
   <ItemGroup>
     <ProjectReference Include="..\..\packages\cultmath\src\CultMath\CultMath.csproj" />
     <ProjectReference Include="..\GameCult.Caching\GameCult.Caching.csproj" />
+    <ProjectReference Include="..\GameCult.Caching.MessagePack\GameCult.Caching.MessagePack.csproj" />
     <PackageReference Include="MessagePack" Version="3.1.7" />
     <PackageReference Include="System.Text.Json" Version="9.0.10" />
   </ItemGroup>

tests/GameCult.Caching.Tests/CultMathJsonTests.cs

@@ -43,6 +43,8 @@ namespace GameCult.Caching.Tests
             Json(new rect(5f, 7f, -1f, -2f)).Should().Be("""{"min":{"x":-1,"y":-2},"max":{"x":5,"y":7}}""");
             Json(new CultCellular(new float4(1f, 2f, 3f, 4f), new float4(5f, 6f, 7f, 8f), 0.25f)).Should().Be(
                 """{"nearest":{"x":1,"y":2,"z":3,"w":4},"edge":{"x":5,"y":6,"z":7,"w":8},"id":0.25}""");
+            Json(new CultPhasor(new float4(1f, 2f, 3f, 4f), new float4(5f, 6f, 7f, 8f))).Should().Be(
+                """{"cos":{"x":1,"y":2,"z":3,"w":4},"sin":{"x":5,"y":6,"z":7,"w":8}}""");
         }
 
         [Test]
@@ -66,6 +68,7 @@ namespace GameCult.Caching.Tests
             RoundTrip(new float3x3(1f, 2f, 3f, 4f, 5f, 6f, 7f, 8f, 9f));
             RoundTrip(new rect(-1f, -2f, 5f, 7f));
             RoundTrip(new CultCellular(new float4(1, 2, 3, 4), new float4(5, 6, 7, 8), 0.25f));
+            RoundTrip(new CultPhasor(new float4(1, 2, 3, 4), new float4(5, 6, 7, 8)));
         }
 
         [Test]

tests/GameCult.Caching.Tests/CultMathSerializationTests.cs

@@ -101,6 +101,7 @@ namespace GameCult.Caching.Tests
             RoundTrip(new float3x3(1f, 2f, 3f, 4f, 5f, 6f, 7f, 8f, 9f));
             RoundTrip(new rect(-1f, -2f, 5f, 7f));
             RoundTrip(new CultCellular(new float4(1, 2, 3, 4), new float4(5, 6, 7, 8), 0.25f));
+            RoundTrip(new CultPhasor(new float4(1, 2, 3, 4), new float4(5, 6, 7, 8)));
         }
 
         [Test]
@@ -132,6 +133,10 @@ namespace GameCult.Caching.Tests
                 "94CA3F800000CA40000000CA40400000CA40800000" +
                 "94CA40A00000CA40C00000CA40E00000CA41000000" +
                 "CA3E800000");
+            Hex(new CultPhasor(new float4(1f, 2f, 3f, 4f), new float4(5f, 6f, 7f, 8f))).Should().Be(
+                "92" +
+                "94CA3F800000CA40000000CA40400000CA40800000" +
+                "94CA40A00000CA40C00000CA40E00000CA41000000");
         }
 
         [Test]

tests/GameCult.Geometry.Tests/GeometryDocumentTests.cs

@@ -10,6 +10,7 @@ using R3;
 using System;
 using System.Collections.Generic;
 using System.IO;
+using System.Linq;
 using System.Text.Json;
 using System.Threading.Tasks;
 
@@ -104,7 +105,7 @@ namespace GameCult.Geometry.Tests
             var domainKey = CultGeometryDomainDocument.CreateRecordKey(domain);
             var request = new CultGeometryBuildRequest
             {
-                DomainKey = domainKey.Value,
+                DomainKey = new CultRecordRef<CultGeometryDomainDocument>(domainKey),
                 WorkerGroup = "ragnarok-column-workers",
                 CameraPosition = [36f, -42f, 30f],
                 FrustumMin = [-75f, -75f, -30f],
@@ -129,6 +130,26 @@ namespace GameCult.Geometry.Tests
             first.Value.Should().StartWith("geometry:request:");
         }
 
+        // The geometry documents reference each other by record key. The reference members are
+        // CultRecordRef<T> so the cache can walk them, and they stay a bare msgpack string on the wire.
+        [Test]
+        public void GeometryReferences_AreWalkableByTheCache_AndStayStringsOnTheWire()
+        {
+            var cache = new CultCache();
+            var descriptor = cache.Registry.GetRequired<CultGeometryBuildRequest>();
+            descriptor.DeclaredMembers.Should().Contain(m => m.MemberName == "DomainKey" && m.IsReference && !m.IsMany);
+
+            var chunk = SampleChunk();
+            var references = cache.Registry.GetRequired<CultGeometryChunkArtifact>().ReferencesOf(chunk, "CutKey").ToArray();
+            references.Should().ContainSingle().Which.Target.Value.Should().Be("geometry:cut:test");
+
+            var payload = CultDocumentMessagePackSerialization.SerializeUntyped(chunk, typeof(CultGeometryChunkArtifact));
+            var reader = new MessagePackReader(payload);
+            reader.ReadArrayHeader();
+            reader.ReadString().Should().Be(chunk.ChunkId);
+            reader.ReadString().Should().Be("geometry:cut:test");
+        }
+
         [Test]
         public void ChunkArtifact_RoundTrips_ThroughCultCacheMessagePack()
         {
@@ -382,7 +403,7 @@ namespace GameCult.Geometry.Tests
             return new CultGeometryChunkArtifact
             {
                 ChunkId = "chunk/ragnarok-column/column-00",
-                CutKey = "geometry:cut:test",
+                CutKey = new CultRecordRef<CultGeometrySelectedCutManifest>(new CultRecordKey("geometry:cut:test")),
                 SelectedCutId = "cut-test",
                 BoundsMin = [-1f, -1f, 0f],
                 BoundsMax = [1f, 1f, 1f],