CultMath shaders: HLSL twins of the four lattice steps, both copies, regenerated GLSL and parity fixture (cut lattice-weighted-shaders)

Co-Authored-By: Claude Sonnet 5.5 <noreply@anthropic.com>

packages/cultmath/README.md

@@ -34,10 +34,12 @@ shader semantics.
   `soft_gabriel_gate` gates a pair of sites by the empty-disc margin. Each is a loop over the shader steps
   `cultmath_power_score` / `cultmath_multiplicative_score`, `cultmath_site_cells_empty`, `cultmath_site_cells_insert`,
   `cultmath_site_cells_border`, `cultmath_power_border_step` (the exact power cell border, `power_border`) and `cultmath_gabriel_margin`; a shader loops over its own buffer.
-- Weighted lattice cells (C# only until the shader cut): `weighted_lattice_cells` returns `float4(F1, F2, smooth, id)` over a jittered,
+- Weighted lattice cells: `weighted_lattice_cells` returns `float4(F1, F2, smooth, id)` over a jittered,
   drifting lattice whose cell weights come from the exact power border of the body span at each cell centre, with a stretched
   metric and an exact ring search over the stated range (`LatticeMaxBoost`, `LatticeMaxStretch`, `LatticeMaxSmoothing`). It composes the
-  public steps `lattice_site`, `lattice_weight`, `lattice_distance` and `lattice_smooth_step`.
+  public steps `lattice_site`, `lattice_weight`, `lattice_distance` and `lattice_smooth_step`. The shaders carry the four steps as
+  `cultmath_lattice_site`, `cultmath_lattice_weight`, `cultmath_lattice_distance` and `cultmath_lattice_smooth_step`; the ring and box loops
+  are the shader caller's, written as the comment above them in `CultMath.hlsl` gives them, and the argument checks become its preconditions.
 - `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.

packages/cultmath/shaders/CultMath.glsl

@@ -688,6 +688,52 @@ float cultmath_gabriel_margin(vec3 a, vec3 b, vec3 other)
     return (cultmath_distance(other, (a + b) * 0.5) - r) / r;
 }
 
+// Weighted lattice cells: the four per-cell steps of math.Lattice.cs, statement for statement. A shader function cannot take a
+// span (design.md, "Sites"), so weighted_lattice_cells is the consumer's loop over its own body buffer, in two passes:
+//   vec2 q = p / cellSize; ivec2 c0 = (ivec2)floor(q); float f1 = 1.0e30, f2 = 1.0e30; int reach = 0;
+//   for (int k = 0; ; k++) {   // ring pass: Chebyshev rings 0, 1, 2, ... around c0, exact exit
+//       if (k > 1 && 0.5 * (k - 1) >= max(f2, f1 + 2.0 * smoothing)) break;   // k never reaches 38 in the accepted range
+//       for each cell on ring k: d = cultmath_lattice_distance(q, cultmath_lattice_site(cell, time, drift).xy, w, normal, stretch),
+//           w = cultmath_lattice_weight(site.z, border, boost, falloff), border = the .w of cultmath_power_border_step
+//           folded over the bodies (owner = the nearest body of the cell centre (cell + 0.5) * cellSize, z = 0);
+//           if (d < f1) { f2 = f1; f1 = d; } else if (d < f2) f2 = d;   store d for the box pass
+//       reach = k; }
+//   float m = f1 + smoothing;   // box pass: every stored cell, rows dz = -reach..reach, then dx = -reach..reach
+//   for each (dz, dx): m = cultmath_lattice_smooth_step(m, d, f1, smoothing);   // strict < for the winning id (earlier cell wins ties)
+// The preconditions the C# function throws for are the caller's in a shader: finite arguments, cellSize > 0, boost in [0, 4],
+// falloff > 0, smoothing in (0, 1], stretch in [1, 2], p within 2^23 cells. The normal is normalised as C# does: divide by its
+// largest component, then by its length (a zero normal means (1, 0) and needs stretch == 1). The ring bound is 37.
+vec3 cultmath_lattice_site(vec2 cell, float time, float drift)
+{
+    ivec3 h = cultmath_pcg3d(ivec3(int(cell.x), int(cell.y), 0));
+    vec2 phase = vec2(float(uint(h.x) >> 8) * (1.0 / 16777216.0), float(uint(h.y) >> 8) * (1.0 / 16777216.0));
+    vec2 site = cell + 0.5 + 0.5 * vec2(sin(drift * time + CULTMATH_TAU * phase.x), sin(drift * time + CULTMATH_TAU * phase.y));
+    return vec3(site, 0.5 + 0.5 * (float(uint(h.z) >> 8) * (1.0 / 16777216.0)));
+}
+
+float cultmath_lattice_weight(float base_weight, float border, float boost, float falloff)
+{
+    float t = border / falloff;
+    return base_weight * (1.0 + boost * exp(-(t * t)));
+}
+
+float cultmath_lattice_distance(vec2 q, vec2 site, float weight, vec2 normal, float stretch)
+{
+    vec2 v = site - q;
+    float a = dot(v, normal);
+    float c = v.x * normal.y - v.y * normal.x;
+    return weight * length(vec2(stretch * a, c));
+}
+
+float cultmath_lattice_smooth_step(float m, float d, float f1, float smoothing)
+{
+    if (d >= f1 + 2.0 * smoothing)
+        return m;
+    float h = cultmath_smoothstep(0.0, 1.0, 0.5 + 0.5 * (m - d) / smoothing);
+    float next = m + (d - m) * h - h * (1.0 - h) * smoothing / (1.0 + 3.0 * smoothing);
+    return min(next, f1 + smoothing);
+}
+
 vec2 cultmath_segment_closest(vec2 point_value, vec2 start, vec2 end)
 {
     vec2 segment = end - start;

packages/cultmath/shaders/CultMath.hlsl

@@ -688,6 +688,52 @@ float cultmath_gabriel_margin(float3 a, float3 b, float3 other)
     return (cultmath_distance(other, (a + b) * 0.5) - r) / r;
 }
 
+// Weighted lattice cells: the four per-cell steps of math.Lattice.cs, statement for statement. A shader function cannot take a
+// span (design.md, "Sites"), so weighted_lattice_cells is the consumer's loop over its own body buffer, in two passes:
+//   float2 q = p / cellSize; int2 c0 = (int2)floor(q); float f1 = 1.0e30, f2 = 1.0e30; int reach = 0;
+//   for (int k = 0; ; k++) {   // ring pass: Chebyshev rings 0, 1, 2, ... around c0, exact exit
+//       if (k > 1 && 0.5 * (k - 1) >= max(f2, f1 + 2.0 * smoothing)) break;   // k never reaches 38 in the accepted range
+//       for each cell on ring k: d = cultmath_lattice_distance(q, cultmath_lattice_site(cell, time, drift).xy, w, normal, stretch),
+//           w = cultmath_lattice_weight(site.z, border, boost, falloff), border = the .w of cultmath_power_border_step
+//           folded over the bodies (owner = the nearest body of the cell centre (cell + 0.5) * cellSize, z = 0);
+//           if (d < f1) { f2 = f1; f1 = d; } else if (d < f2) f2 = d;   store d for the box pass
+//       reach = k; }
+//   float m = f1 + smoothing;   // box pass: every stored cell, rows dz = -reach..reach, then dx = -reach..reach
+//   for each (dz, dx): m = cultmath_lattice_smooth_step(m, d, f1, smoothing);   // strict < for the winning id (earlier cell wins ties)
+// The preconditions the C# function throws for are the caller's in a shader: finite arguments, cellSize > 0, boost in [0, 4],
+// falloff > 0, smoothing in (0, 1], stretch in [1, 2], p within 2^23 cells. The normal is normalised as C# does: divide by its
+// largest component, then by its length (a zero normal means (1, 0) and needs stretch == 1). The ring bound is 37.
+float3 cultmath_lattice_site(float2 cell, float time, float drift)
+{
+    int3 h = cultmath_pcg3d(int3((int)cell.x, (int)cell.y, 0));
+    float2 phase = float2((float)((uint)h.x >> 8) * (1.0 / 16777216.0), (float)((uint)h.y >> 8) * (1.0 / 16777216.0));
+    float2 site = cell + 0.5 + 0.5 * float2(sin(drift * time + CULTMATH_TAU * phase.x), sin(drift * time + CULTMATH_TAU * phase.y));
+    return float3(site, 0.5 + 0.5 * ((float)((uint)h.z >> 8) * (1.0 / 16777216.0)));
+}
+
+float cultmath_lattice_weight(float base_weight, float border, float boost, float falloff)
+{
+    float t = border / falloff;
+    return base_weight * (1.0 + boost * exp(-(t * t)));
+}
+
+float cultmath_lattice_distance(float2 q, float2 site, float weight, float2 normal, float stretch)
+{
+    float2 v = site - q;
+    float a = dot(v, normal);
+    float c = v.x * normal.y - v.y * normal.x;
+    return weight * length(float2(stretch * a, c));
+}
+
+float cultmath_lattice_smooth_step(float m, float d, float f1, float smoothing)
+{
+    if (d >= f1 + 2.0 * smoothing)
+        return m;
+    float h = cultmath_smoothstep(0.0, 1.0, 0.5 + 0.5 * (m - d) / smoothing);
+    float next = m + (d - m) * h - h * (1.0 - h) * smoothing / (1.0 + 3.0 * smoothing);
+    return min(next, f1 + smoothing);
+}
+
 float2 cultmath_segment_closest(float2 point_value, float2 start, float2 end)
 {
     float2 segment = end - start;

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

@@ -852,6 +852,34 @@ static int4 clamp(int4 x, int4 a, int4 b) { if (a.x > b.x || a.y > b.y || a.z >
             var a = Point3(r);
             return new object[] { a, r.Next(8) == 0 ? a : Point3(r), Point3(r) };
         }),
+        // The lattice steps. The cell is an integer-valued vec2 (negative and large included), time and drift span
+        // zero and the accepted range; the weight's border covers 0, near the falloff, and the 1.0e30 of no body.
+        new Family("vec3 cultmath_lattice_site(vec2, float, float)", "ulp-bounded", r =>
+            new object[] { new float2(r.Next(-4000, 4000), r.Next(-4000, 4000)), r.Next(6) == 0 ? 0.0f : Uniform(r, -200.0f, 200.0f), r.Next(6) == 0 ? 0.0f : Uniform(r, -3.0f, 3.0f) }),
+        new Family("float cultmath_lattice_weight(float, float, float, float)", "ulp-bounded", r =>
+        {
+            var falloff = Uniform(r, 0.05f, 40.0f);
+            var border = r.Next(5) switch { 0 => 0.0f, 1 => 1.0e30f, 2 => falloff * Uniform(r, 0.9f, 1.1f), _ => Uniform(r, 0.0f, 3.0f * falloff) };
+            var boost = r.Next(5) switch { 0 => 0.0f, 1 => 4.0f, _ => Uniform(r, 0.0f, 4.0f) };
+            return new object[] { Uniform(r, 0.5f, 1.0f), border, boost, falloff };
+        }),
+        new Family("float cultmath_lattice_distance(vec2, vec2, float, vec2, float)", "ulp-bounded", r =>
+        {
+            var angle = Uniform(r, 0.0f, 6.2831855f);
+            var normal = r.Next(6) == 0 ? new float2(1.0f, 0.0f) : new float2(MathF.Cos(angle), MathF.Sin(angle));
+            var stretch = r.Next(5) switch { 0 => 1.0f, 1 => 2.0f, _ => Uniform(r, 1.0f, 2.0f) };
+            return new object[] { Point2(r), Point2(r), Uniform(r, 0.5f, 5.0f), normal, stretch };
+        }),
+        // The smooth step: d exactly at f1 + 2 smoothing one in four (the equal branch), near it one in four, the
+        // smoothing near 0 and at 1, m anywhere up to f1 + smoothing.
+        new Family("float cultmath_lattice_smooth_step(float, float, float, float)", "ulp-bounded", r =>
+        {
+            var smoothing = r.Next(5) switch { 0 => 1.0f, 1 => Uniform(r, 1.0e-3f, 1.0e-2f), _ => Uniform(r, 0.01f, 1.0f) };
+            var f1 = Uniform(r, 0.0f, 10.0f);
+            var edge = f1 + 2.0f * smoothing;
+            var d = r.Next(4) switch { 0 => edge, 1 => edge * Uniform(r, 0.99f, 1.01f), _ => Uniform(r, f1, edge + 1.0f) };
+            return new object[] { Uniform(r, f1, f1 + smoothing), d, f1, smoothing };
+        }),
     }).ToArray();
 
     // A point of magnitude m, m log-uniform in [lo, hi]: one component at +-m, the other two uniform in [-m, m],

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

@@ -35,6 +35,10 @@ public sealed class HlslSourceCompatibilityTests
         Assert.NotNull(ShaderFunction(shader, "cultmath_site_cells_border", siteCells));
         Assert.NotNull(ShaderFunction(shader, "cultmath_gabriel_margin", typeof(float3), typeof(float3), typeof(float3)));
         Assert.NotNull(ShaderFunction(shader, "cultmath_power_border_step", typeof(float3), typeof(float4), typeof(float4), typeof(float4)));
+        Assert.NotNull(ShaderFunction(shader, "cultmath_lattice_site", typeof(float2), typeof(float), typeof(float)));
+        Assert.NotNull(ShaderFunction(shader, "cultmath_lattice_weight", typeof(float), typeof(float), typeof(float), typeof(float)));
+        Assert.NotNull(ShaderFunction(shader, "cultmath_lattice_distance", typeof(float2), typeof(float2), typeof(float), typeof(float2), typeof(float)));
+        Assert.NotNull(ShaderFunction(shader, "cultmath_lattice_smooth_step", typeof(float), typeof(float), typeof(float), typeof(float)));
         Assert.NotNull(ShaderFunction(shader, "cultmath_segment_closest", typeof(float2), typeof(float2), typeof(float2)));
         Assert.NotNull(ShaderFunction(shader, "cultmath_distance_to_segment", typeof(float2), typeof(float2), typeof(float2)));
         Assert.NotNull(ShaderFunction(shader, "cultmath_segment_closest", typeof(float3), typeof(float3), typeof(float3)));
@@ -200,7 +204,7 @@ public sealed class HlslSourceCompatibilityTests
         }
 
         Assert.True(mismatches.Count == 0, $"{mismatches.Count} mismatches:{Environment.NewLine}{string.Join(Environment.NewLine, mismatches.Take(12))}");
-        Assert.Equal(136, compared.Count);
+        Assert.Equal(140, compared.Count);
     }
 
     /// <summary>

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

@@ -688,6 +688,52 @@ float cultmath_gabriel_margin(float3 a, float3 b, float3 other)
     return (cultmath_distance(other, (a + b) * 0.5) - r) / r;
 }
 
+// Weighted lattice cells: the four per-cell steps of math.Lattice.cs, statement for statement. A shader function cannot take a
+// span (design.md, "Sites"), so weighted_lattice_cells is the consumer's loop over its own body buffer, in two passes:
+//   float2 q = p / cellSize; int2 c0 = (int2)floor(q); float f1 = 1.0e30, f2 = 1.0e30; int reach = 0;
+//   for (int k = 0; ; k++) {   // ring pass: Chebyshev rings 0, 1, 2, ... around c0, exact exit
+//       if (k > 1 && 0.5 * (k - 1) >= max(f2, f1 + 2.0 * smoothing)) break;   // k never reaches 38 in the accepted range
+//       for each cell on ring k: d = cultmath_lattice_distance(q, cultmath_lattice_site(cell, time, drift).xy, w, normal, stretch),
+//           w = cultmath_lattice_weight(site.z, border, boost, falloff), border = the .w of cultmath_power_border_step
+//           folded over the bodies (owner = the nearest body of the cell centre (cell + 0.5) * cellSize, z = 0);
+//           if (d < f1) { f2 = f1; f1 = d; } else if (d < f2) f2 = d;   store d for the box pass
+//       reach = k; }
+//   float m = f1 + smoothing;   // box pass: every stored cell, rows dz = -reach..reach, then dx = -reach..reach
+//   for each (dz, dx): m = cultmath_lattice_smooth_step(m, d, f1, smoothing);   // strict < for the winning id (earlier cell wins ties)
+// The preconditions the C# function throws for are the caller's in a shader: finite arguments, cellSize > 0, boost in [0, 4],
+// falloff > 0, smoothing in (0, 1], stretch in [1, 2], p within 2^23 cells. The normal is normalised as C# does: divide by its
+// largest component, then by its length (a zero normal means (1, 0) and needs stretch == 1). The ring bound is 37.
+float3 cultmath_lattice_site(float2 cell, float time, float drift)
+{
+    int3 h = cultmath_pcg3d(int3((int)cell.x, (int)cell.y, 0));
+    float2 phase = float2((float)((uint)h.x >> 8) * (1.0 / 16777216.0), (float)((uint)h.y >> 8) * (1.0 / 16777216.0));
+    float2 site = cell + 0.5 + 0.5 * float2(sin(drift * time + CULTMATH_TAU * phase.x), sin(drift * time + CULTMATH_TAU * phase.y));
+    return float3(site, 0.5 + 0.5 * ((float)((uint)h.z >> 8) * (1.0 / 16777216.0)));
+}
+
+float cultmath_lattice_weight(float base_weight, float border, float boost, float falloff)
+{
+    float t = border / falloff;
+    return base_weight * (1.0 + boost * exp(-(t * t)));
+}
+
+float cultmath_lattice_distance(float2 q, float2 site, float weight, float2 normal, float stretch)
+{
+    float2 v = site - q;
+    float a = dot(v, normal);
+    float c = v.x * normal.y - v.y * normal.x;
+    return weight * length(float2(stretch * a, c));
+}
+
+float cultmath_lattice_smooth_step(float m, float d, float f1, float smoothing)
+{
+    if (d >= f1 + 2.0 * smoothing)
+        return m;
+    float h = cultmath_smoothstep(0.0, 1.0, 0.5 + 0.5 * (m - d) / smoothing);
+    float next = m + (d - m) * h - h * (1.0 - h) * smoothing / (1.0 + 3.0 * smoothing);
+    return min(next, f1 + smoothing);
+}
+
 float2 cultmath_segment_closest(float2 point_value, float2 start, float2 end)
 {
     float2 segment = end - start;