cultmath: measure the Lipschitz bound of phacelle, with its tool and a pinning test

max singular value of the (cosine, sine) Jacobian <= max(10, 1.25 |side|) / (1 - normalization)
for normalization in [0, 0.9], |side| in [0.25, 25]; method and tightness in design.md.

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

packages/cultmath/docs/design.md

@@ -451,6 +451,26 @@ together, laid out as `float4(gradient.xyz, value.w)`, with no value-only twin
   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.
+  Lipschitz bound of the `(cosine, sine)` field in `p` (one cell is 1 unit): the
+  largest singular value of its 2×3 Jacobian is at most
+  `max(10, 1.25 |side|) / (1 - normalization)` for `normalization` in `[0, 0.9]` and
+  `|side|` in `[0.25, 25]`. A scalar use of `cosine` or `sine` alone is bounded by the same number (its gradient norm is
+  one row of the Jacobian, never above the largest singular value). Method (`tools/phacelle-lipschitz`):
+  random `p` (a 64-unit cube), a uniformly random `side` direction and offset, the
+  exact gradient at each, 4 million samples per (normalization, |side|); then
+  hill-climbing from the 256 best samples over `p`, `side` direction and offset
+  (10 000 steps each); then value differences around the climbed argmax to confirm the
+  gradient is the field's slope (agreement within 3%). Random sampling alone
+  undershoots the climbed maximum by up to 4×, because the maxima sit at near-cancellation
+  points of the blend (`1/|I|`); quote the climbed number. The climbed sup is 0.016 to 0.88 of
+  the envelope over the sweep: tightest (within 1.15×) at `normalization` 0 and `|side| >= 8`,
+  and 1.5 to 60× loose at small `|side|` or large `normalization`, where the
+  envelope's constant 10 is the worst case over all `|side|`. The typical gradient (99.99th percentile of the
+  samples) is 2 to 5× below the sup, so a sphere trace steps safely at `distance / L`
+  and can step farther with the returned gradient. The bound is not a proof (a
+  search found it, not an argument) and it does not hold as `normalization` approaches 1:
+  the 1/|I| factor is unbounded there (sampled 1335 at `normalization` 1, `|side|` 2π, and
+  rising with the sample count). If the caller scales `p` by `f` the bound scales by `f`.
   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

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

@@ -337,4 +337,33 @@ public sealed class PhacelleTests
         Assert.Equal(BitConverter.SingleToInt32Bits(expected.z), BitConverter.SingleToInt32Bits(actual.z));
         Assert.Equal(BitConverter.SingleToInt32Bits(expected.w), BitConverter.SingleToInt32Bits(actual.w));
     }
+
+    // The Lipschitz envelope of the (cosine, sine) field in p, per docs/design.md ("Lipschitz bound"): the largest
+    // singular value of its 2x3 Jacobian never exceeds max(10, 1.25 |side|) / (1 - normalization), for normalization
+    // in [0, 0.9] and |side| in [0.25, 25]. tools/phacelle-lipschitz measures it by sampling and hill-climbing; this
+    // pins that no ordinary sample passes it (the climbed maximum is at most 0.88 of the envelope).
+    [Theory]
+    [InlineData(0.0f)]
+    [InlineData(0.5f)]
+    [InlineData(0.75f)]
+    [InlineData(0.9f)]
+    public void GradientNormStaysUnderTheDocumentedLipschitzEnvelope(float normalization)
+    {
+        var random = new System.Random(0x11B5);
+        foreach (var sideLength in new[] { 0.25f, 1.0f, 1.5f, 4.0f, 12.5f, 25.0f })
+        {
+            var envelope = Math.Max(10.0, 1.25 * sideLength) / (1.0 - normalization);
+            for (var i = 0; i < 20000; i++)
+            {
+                var p = RandomPoint(random, 32.0f);
+                var side = math.normalize(RandomSide(random)) * sideLength;
+                var phasor = math.phacelle(p, side, random.NextSingle(), normalization);
+                var a = (double)math.dot(phasor.cosine.xyz, phasor.cosine.xyz);
+                var d = (double)math.dot(phasor.sine.xyz, phasor.sine.xyz);
+                var b = (double)math.dot(phasor.cosine.xyz, phasor.sine.xyz);
+                var norm = Math.Sqrt((a + d) / 2.0 + Math.Sqrt(Math.Max(0.0, (a + d) * (a + d) / 4.0 - (a * d - b * b))));
+                Assert.True(norm <= envelope, $"|J| {norm} passes the envelope {envelope} at p {p}, side {side}, normalization {normalization}");
+            }
+        }
+    }
 }

packages/cultmath/tools/phacelle-lipschitz/Program.cs

@@ -0,0 +1,85 @@
+// Measures the Lipschitz constant of math.phacelle's (cosine, sine) field, the number docs/design.md quotes and
+// PhacelleTests.GradientNormStaysUnderTheDocumentedLipschitzEnvelope pins.
+//   dotnet run -c Release --project packages/cultmath/tools/phacelle-lipschitz -- <samples per cell> <normalizations> <side lengths>
+//   e.g. ... -- 4000000 0,0.5,0.75,0.9 0.25,0.5,1,1.5,2,4,8,12.5,25
+// Columns: sampledMax (random p), climbedMax (hill-climbed from the 256 best samples; the number to quote), secantMax
+// (value differences around the climbed argmax, the check that the analytic gradient is the field's slope),
+// scalarMax (the larger of |grad cosine|, |grad sine| over the samples), p99.99 (the typical tail), and climbed*(1-n)/|side|.
+using System;
+using System.Linq;
+using System.Threading.Tasks;
+using CultMath;
+
+// Lipschitz bound of math.phacelle's (cosine, sine) field: sup over p of the operator norm of its 2x3 Jacobian,
+// found by (1) random sampling of the exact analytic gradient, (2) hill-climbing from the top samples over
+// (p, off, side direction), (3) a secant check on the VALUES (|f(p+d) - f(p)| / |d|) so the gradient itself is not trusted.
+static double Op(float3 p, float3 side, float off, float n, out double comp)
+{
+    var r = math.phacelle(p, side, off, n);
+    double a = r.cosine.x*r.cosine.x + r.cosine.y*r.cosine.y + r.cosine.z*r.cosine.z;
+    double d = r.sine.x*r.sine.x + r.sine.y*r.sine.y + r.sine.z*r.sine.z;
+    double b = r.cosine.x*r.sine.x + r.cosine.y*r.sine.y + r.cosine.z*r.sine.z;
+    double tr = a + d, det = a*d - b*b;
+    comp = Math.Sqrt(Math.Max(a, d));
+    var op = Math.Sqrt(tr/2 + Math.Sqrt(Math.Max(0, tr*tr/4 - det)));
+    return double.IsFinite(op) ? op : 0;
+}
+static float3 Dir(double u, double th, double len) { var s = Math.Sqrt(Math.Max(0, 1-u*u)); return new float3((float)(s*Math.Cos(th)*len), (float)(s*Math.Sin(th)*len), (float)(u*len)); }
+
+long N = long.Parse(args[0]);
+var norms = args[1].Split(',').Select(float.Parse).ToArray();
+var sides = args[2].Split(',').Select(float.Parse).ToArray();
+Console.WriteLine("n\t|side|\tN\tsampledMax\tclimbedMax\tsecantMax\tscalarMax(cos|sin)\tp99.99\tclimbed*(1-n)/|side|");
+foreach (var n in norms) foreach (var sl in sides)
+{
+    int T = Environment.ProcessorCount;
+    var tops = new (double v, float3 p, double u, double th, float off)[T][];
+    var allv = new double[T][]; var scal = new double[T];
+    Parallel.For(0, T, t =>
+    {
+        var rng = new System.Random(unchecked(99991 * t + (int)(n * 1000) * 131 + (int)(sl * 100)));
+        var top = new (double v, float3 p, double u, double th, float off)[16];
+        var list = new double[N / T];
+        double sc = 0;
+        for (long i = 0; i < N / T; i++)
+        {
+            var p = new float3((float)(rng.NextDouble()*64-32), (float)(rng.NextDouble()*64-32), (float)(rng.NextDouble()*64-32));
+            double u = rng.NextDouble()*2-1, th = rng.NextDouble()*Math.PI*2; float off = (float)rng.NextDouble();
+            var v = Op(p, Dir(u, th, sl), off, n, out var c);
+            list[i] = v; if (c > sc) sc = c;
+            int m = 0; for (int k = 1; k < 16; k++) if (top[k].v < top[m].v) m = k;
+            if (v > top[m].v) top[m] = (v, p, u, th, off);
+        }
+        tops[t] = top; allv[t] = list; scal[t] = sc;
+    });
+    var all = allv.SelectMany(x => x).ToArray(); Array.Sort(all);
+    double sampled = all[^1], p9999 = all[(int)(all.Length * 0.9999)];
+    // hill-climb the best 16*T by random local search
+    var seeds = tops.SelectMany(x => x).OrderByDescending(x => x.v).Take(256).ToArray();
+    var climbed = new (double v, float3 p, double u, double th, float off)[seeds.Length];
+    Parallel.For(0, seeds.Length, s =>
+    {
+        var rng = new System.Random(s + 17);
+        var cur = seeds[s]; double step = 0.05;
+        for (int it = 0; it < 10000; it++)
+        {
+            var q = new float3(cur.p.x + (float)((rng.NextDouble()-.5)*step), cur.p.y + (float)((rng.NextDouble()-.5)*step), cur.p.z + (float)((rng.NextDouble()-.5)*step));
+            var u = Math.Clamp(cur.u + (rng.NextDouble()-.5)*step, -1, 1); var th = cur.th + (rng.NextDouble()-.5)*step; var off = cur.off + (float)((rng.NextDouble()-.5)*step);
+            var v = Op(q, Dir(u, th, sl), off, n, out _);
+            if (v > cur.v) cur = (v, q, u, th, off); else if (it % 1000 == 999) step *= 0.5;
+        }
+        climbed[s] = cur;
+    });
+    var best = climbed.OrderByDescending(x => x.v).First();
+    // secant check around the climbed argmax: values only, tiny displacements in random directions
+    double secant = 0; var sr = new System.Random(5);
+    var side0 = Dir(best.u, best.th, sl);
+    for (int i = 0; i < 20000; i++)
+    {
+        double e = Math.Pow(10, -2 - 2*sr.NextDouble()); var d = Dir(sr.NextDouble()*2-1, sr.NextDouble()*6.2832, e);
+        var a = math.phacelle(best.p, side0, best.off, n); var b = math.phacelle(best.p + d, side0, best.off, n);
+        double dc = b.cosine.w - a.cosine.w, ds = b.sine.w - a.sine.w; var dl = Math.Sqrt((double)d.x*d.x + (double)d.y*d.y + (double)d.z*d.z);
+        secant = Math.Max(secant, Math.Sqrt(dc*dc + ds*ds) / dl);
+    }
+    Console.WriteLine($"{n}\t{sl}\t{N}\t{sampled:F3}\t{best.v:F3}\t{secant:F3}\t{scal.Max():F3}\t{p9999:F3}\t{best.v*(1-n)/sl:F3}");
+}

packages/cultmath/tools/phacelle-lipschitz/phacelle-lipschitz.csproj

@@ -0,0 +1,4 @@
+<Project Sdk="Microsoft.NET.Sdk">
+  <PropertyGroup><OutputType>Exe</OutputType><TargetFramework>net10.0</TargetFramework><Nullable>enable</Nullable><AllowUnsafeBlocks>true</AllowUnsafeBlocks><InvariantGlobalization>true</InvariantGlobalization><ServerGarbageCollector>false</ServerGarbageCollector></PropertyGroup>
+  <ItemGroup><ProjectReference Include="../../src/CultMath/CultMath.csproj" /></ItemGroup>
+</Project>