Affine forms in HLSL: CultMath.Affine.hlsl mirrors math.Affine.cs bit for bit; iv_snoise_ball(float2); the GLSL lowering carries them

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

packages/cultmath/README.md

@@ -34,7 +34,8 @@ shader semantics.
   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) and
-  `shaders/CultMath.Interval.hlsl` from the same directory, so ship all three.
+  `shaders/CultMath.Interval.hlsl` and `shaders/CultMath.Affine.hlsl` from the same
+  directory, so ship all four.
 - interval arithmetic for culling empty space with a proof: `iv_*` functions
   over `float2(lo, hi)` (`iv_add`, `iv_mul`, `iv_exp`, `iv_smoothstep`, ...)
   and `iv_snoise_ball`/`iv_fbm_ball`, which bound `snoise` and `fbm_grad` over a

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

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

packages/cultmath/shaders/CultMath.Affine.hlsl

@@ -0,0 +1,130 @@
+#ifndef CULTMATH_AFFINE_HLSL
+#define CULTMATH_AFFINE_HLSL
+
+// Reduced affine forms with one shared symbol over float3(x0, a, e): the set { x0 + a eps + e delta :
+// delta in [-1, 1] } at the region's shared eps in [-1, 1], e >= 0. Every cultmath_af_* function mirrors
+// math.Affine.cs bit for bit; see that file and docs/design.md ("Affine forms") for each op's enclosure
+// and the float32 rounding allowance it folds into e. Included by CultMath.hlsl after
+// CultMath.Interval.hlsl, whose CULTMATH_SNOISE_LIPSCHITZ it shares, and after its own
+// definition of cultmath_snoise_grad.
+
+// An upper bound on the spectral norm of the Hessian of cultmath_snoise(float3); provenance in
+// math.Affine.cs (SNOISE_HESSIAN).
+static const float CULTMATH_SNOISE_HESSIAN = 56.050385;
+
+float3 cultmath_af_point(float value) { return float3(value, 0.0, 0.0); }
+float3 cultmath_af_symbol(float x0, float a) { return float3(x0, a, 0.0); }
+
+float3 cultmath_af_from_iv(float2 interval)
+{
+    float x0 = interval.x * 0.5 + interval.y * 0.5;
+    float e = interval.y * 0.5 - interval.x * 0.5;
+    return float3(x0, 0.0, e + (abs(x0) + e) * 9.5367431640625e-7 + 1.17549435e-38);
+}
+
+// An endpoint that would be NaN is the infinity on its side.
+float2 cultmath_af_range(float3 x)
+{
+    float r = abs(x.y) + x.z;
+    float w = abs(x.x) * 9.5367431640625e-7 + r * 9.5367431640625e-7;
+    float lo = x.x - r - w;
+    float hi = x.x + r + w;
+    return float2(lo == lo ? lo : asfloat(0xff800000u), hi == hi ? hi : asfloat(0x7f800000u));
+}
+
+float3 cultmath_af_add(float3 x, float3 y)
+{
+    float x0 = x.x + y.x;
+    float a = x.y + y.y;
+    float e = x.z + y.z;
+    return float3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7);
+}
+
+float3 cultmath_af_sub(float3 x, float3 y)
+{
+    float x0 = x.x - y.x;
+    float a = x.y - y.y;
+    float e = x.z + y.z;
+    return float3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7);
+}
+
+float3 cultmath_af_neg(float3 x) { return float3(-x.x, -x.y, x.z); }
+
+float3 cultmath_af_scale(float3 x, float s)
+{
+    float x0 = x.x * s;
+    float a = x.y * s;
+    float e = x.z * abs(s);
+    return float3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7 + 1.17549435e-38);
+}
+
+float3 cultmath_af_add_iv(float3 x, float2 interval) { return cultmath_af_add(x, cultmath_af_from_iv(interval)); }
+
+float3 cultmath_af_mul(float3 x, float3 y)
+{
+    float x0 = x.x * y.x;
+    float p = x.x * y.y;
+    float q = y.x * x.y;
+    float a = p + q;
+    float e = abs(x.x) * y.z + abs(y.x) * x.z + (abs(x.y) + x.z) * (abs(y.y) + y.z);
+    return float3(x0, a, e + (abs(x0) + abs(p) + abs(q) + e) * 9.5367431640625e-7 + 1.17549435e-38);
+}
+
+// The narrower of the Lipschitz form and the centred (Taylor) form over every p = centre + axis eps + w,
+// |w| <= radius, from one cultmath_snoise_grad at the centre; the allowance for the float32 evaluation is
+// in math.Affine.cs.
+float3 cultmath_af_snoise(float3 centre, float3 axis, float radius)
+{
+    float4 g = cultmath_snoise_grad(centre);
+    float3 gradient = float3(g.x, g.y, g.z);
+    float reach = length(axis) + radius;
+    float lipschitz = CULTMATH_SNOISE_LIPSCHITZ * reach;
+    float a = dot(gradient, axis);
+    float remainder = length(gradient) * radius + 0.5 * CULTMATH_SNOISE_HESSIAN * reach * reach;
+    float allowance = 6.103515625e-5 + 2.44140625e-4 * reach
+        + CULTMATH_SNOISE_LIPSCHITZ * 9.5367431640625e-7 * (abs(centre.x) + abs(centre.y) + abs(centre.z) + reach);
+    bool centred = abs(a) + remainder < lipschitz;
+    float slope = centred ? a : 0.0;
+    float e = (centred ? remainder : lipschitz) + allowance;
+    return float3(g.w, slope, e + (abs(g.w) + abs(slope) + e) * 9.5367431640625e-7);
+}
+
+// Octave sum of cultmath_af_snoise, compounding frequency and amplitude as the value-and-gradient fBm does.
+float3 cultmath_af_fbm(float3 centre, float3 axis, float radius, int octaves, float lacunarity, float gain)
+{
+    octaves = clamp(octaves, 0, 16);
+    float amplitude = 1.0;
+    float frequency = 1.0;
+    float3 sum = float3(0.0, 0.0, 0.0);
+    for (int i = 0; i < octaves; i++)
+    {
+        sum = cultmath_af_add(sum, cultmath_af_scale(cultmath_af_snoise(centre * frequency, axis * frequency, radius * abs(frequency)), amplitude));
+        frequency *= lacunarity;
+        amplitude *= gain;
+    }
+
+    return sum;
+}
+
+// The axis of a screen tile's slice over the depth segment [z0, z1]; camera frame of
+// cultmath_iv_frustum_ball.
+float3 cultmath_af_frustum_axis(float2 centreSlope, float z0, float z1)
+{
+    float h = (z1 - z0) * 0.5;
+    return float3(centreSlope.x * h, centreSlope.y * h, h);
+}
+
+// The affine slice of a screen tile over [z0, z1]: centre xyz, radius w. Derivation in math.Affine.cs.
+float4 cultmath_af_frustum_ball(float2 centreSlope, float z0, float z1, float footprintPerDepth, float warpVariation)
+{
+    float zm = (z0 + z1) * 0.5;
+    float h = (z1 - z0) * 0.5;
+    float cx = centreSlope.x * zm;
+    float cy = centreSlope.y * zm;
+    float spread = (abs(centreSlope.x) + abs(centreSlope.y) + 1.0) * h;
+    float radius = z1 * footprintPerDepth + warpVariation;
+    radius += (abs(cx) + abs(cy) + abs(zm) + spread + radius) * 9.5367431640625e-7;
+    return float4(cx, cy, zm, radius);
+}
+
+#endif

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

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

packages/cultmath/shaders/CultMath.Interval.hlsl

@@ -9,6 +9,9 @@
 // Lipschitz constant of cultmath_snoise(float3); provenance in math.Interval.cs (SNOISE_LIPSCHITZ).
 static const float CULTMATH_SNOISE_LIPSCHITZ = 7.9640074;
 
+// Lipschitz constant of cultmath_snoise(float2); provenance in math.Interval.cs (SNOISE2_LIPSCHITZ).
+static const float CULTMATH_SNOISE2_LIPSCHITZ = 8.117002;
+
 float2 cultmath_iv_point(float value) { return float2(value, value); }
 float2 cultmath_iv_add(float2 a, float2 b) { return float2(a.x + b.x, a.y + b.y); }
 float2 cultmath_iv_sub(float2 a, float2 b) { return float2(a.x - b.y, a.y - b.x); }
@@ -87,6 +90,14 @@ float2 cultmath_iv_snoise_ball(float3 centre, float radius)
     return float2(max(n - e, -1.0), min(n + e, 1.0));
 }
 
+// The 2D overload: [n - L2 r, n + L2 r] intersected with [-1, 1], n = cultmath_snoise(centre). One snoise, value only.
+float2 cultmath_iv_snoise_ball(float2 centre, float radius)
+{
+    float n = cultmath_snoise(centre);
+    float e = CULTMATH_SNOISE2_LIPSCHITZ * radius;
+    return float2(max(n - e, -1.0), min(n + e, 1.0));
+}
+
 // Octave sum of cultmath_iv_snoise_ball, compounding frequency and amplitude as the value-and-gradient fBm does.
 float2 cultmath_iv_fbm_ball(float3 centre, float radius, int octaves, float lacunarity, float gain)
 {

packages/cultmath/shaders/CultMath.glsl

@@ -1,4 +1,4 @@
-// Generated from shaders/CultMath.hlsl and CultMath.Interval.hlsl by GlslLowering (packages/cultmath/tests/CultMath.Tests/GlslLowering.cs). Do not edit; regenerate with CULTMATH_WRITE_GLSL=1 dotnet test --filter GlslMirrorTests. MIT. cultmath_phacelle is in CultMath.Phacelle.glsl, to be concatenated after this file by consumers that call it.
+// Generated from shaders/CultMath.hlsl and CultMath.Interval.hlsl and CultMath.Affine.hlsl by GlslLowering (packages/cultmath/tests/CultMath.Tests/GlslLowering.cs). Do not edit; regenerate with CULTMATH_WRITE_GLSL=1 dotnet test --filter GlslMirrorTests. MIT. cultmath_phacelle is in CultMath.Phacelle.glsl, to be concatenated after this file by consumers that call it.
 #ifndef CULTMATH_GLSL
 #define CULTMATH_GLSL
 // CultMath HLSL mirror.
@@ -518,6 +518,9 @@ CultCellular cultmath_cellular(vec3 p)
 // Lipschitz constant of cultmath_snoise(vec3); provenance in math.Interval.cs (SNOISE_LIPSCHITZ).
 const float CULTMATH_SNOISE_LIPSCHITZ = 7.9640074;
 
+// Lipschitz constant of cultmath_snoise(vec2); provenance in math.Interval.cs (SNOISE2_LIPSCHITZ).
+const float CULTMATH_SNOISE2_LIPSCHITZ = 8.117002;
+
 vec2 cultmath_iv_point(float value) { return vec2(value, value); }
 vec2 cultmath_iv_add(vec2 a, vec2 b) { return vec2(a.x + b.x, a.y + b.y); }
 vec2 cultmath_iv_sub(vec2 a, vec2 b) { return vec2(a.x - b.y, a.y - b.x); }
@@ -596,6 +599,14 @@ vec2 cultmath_iv_snoise_ball(vec3 centre, float radius)
     return vec2(max(n - e, -1.0), min(n + e, 1.0));
 }
 
+// The 2D overload: [n - L2 r, n + L2 r] intersected with [-1, 1], n = cultmath_snoise(centre). One snoise, value only.
+vec2 cultmath_iv_snoise_ball(vec2 centre, float radius)
+{
+    float n = cultmath_snoise(centre);
+    float e = CULTMATH_SNOISE2_LIPSCHITZ * radius;
+    return vec2(max(n - e, -1.0), min(n + e, 1.0));
+}
+
 // Octave sum of cultmath_iv_snoise_ball, compounding frequency and amplitude as the value-and-gradient fBm does.
 vec2 cultmath_iv_fbm_ball(vec3 centre, float radius, int octaves, float lacunarity, float gain)
 {
@@ -631,6 +642,137 @@ vec4 cultmath_iv_frustum_ball(vec2 centreSlope, float z0, float z1, float footpr
 
 
 
+
+
+// Reduced affine forms with one shared symbol over vec3(x0, a, e): the set { x0 + a eps + e delta :
+// delta in [-1, 1] } at the region's shared eps in [-1, 1], e >= 0. Every cultmath_af_* function mirrors
+// math.Affine.cs bit for bit; see that file and docs/design.md ("Affine forms") for each op's enclosure
+// and the float32 rounding allowance it folds into e. Included by CultMath.hlsl after
+// CultMath.Interval.hlsl, whose CULTMATH_SNOISE_LIPSCHITZ it shares, and after its own
+// definition of cultmath_snoise_grad.
+
+// An upper bound on the spectral norm of the Hessian of cultmath_snoise(vec3); provenance in
+// math.Affine.cs (SNOISE_HESSIAN).
+const float CULTMATH_SNOISE_HESSIAN = 56.050385;
+
+vec3 cultmath_af_point(float value) { return vec3(value, 0.0, 0.0); }
+vec3 cultmath_af_symbol(float x0, float a) { return vec3(x0, a, 0.0); }
+
+vec3 cultmath_af_from_iv(vec2 interval)
+{
+    float x0 = interval.x * 0.5 + interval.y * 0.5;
+    float e = interval.y * 0.5 - interval.x * 0.5;
+    return vec3(x0, 0.0, e + (abs(x0) + e) * 9.5367431640625e-7 + 1.17549435e-38);
+}
+
+// An endpoint that would be NaN is the infinity on its side.
+vec2 cultmath_af_range(vec3 x)
+{
+    float r = abs(x.y) + x.z;
+    float w = abs(x.x) * 9.5367431640625e-7 + r * 9.5367431640625e-7;
+    float lo = x.x - r - w;
+    float hi = x.x + r + w;
+    return vec2(lo == lo ? lo : uintBitsToFloat(0xff800000u), hi == hi ? hi : uintBitsToFloat(0x7f800000u));
+}
+
+vec3 cultmath_af_add(vec3 x, vec3 y)
+{
+    float x0 = x.x + y.x;
+    float a = x.y + y.y;
+    float e = x.z + y.z;
+    return vec3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7);
+}
+
+vec3 cultmath_af_sub(vec3 x, vec3 y)
+{
+    float x0 = x.x - y.x;
+    float a = x.y - y.y;
+    float e = x.z + y.z;
+    return vec3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7);
+}
+
+vec3 cultmath_af_neg(vec3 x) { return vec3(-x.x, -x.y, x.z); }
+
+vec3 cultmath_af_scale(vec3 x, float s)
+{
+    float x0 = x.x * s;
+    float a = x.y * s;
+    float e = x.z * abs(s);
+    return vec3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7 + 1.17549435e-38);
+}
+
+vec3 cultmath_af_add_iv(vec3 x, vec2 interval) { return cultmath_af_add(x, cultmath_af_from_iv(interval)); }
+
+vec3 cultmath_af_mul(vec3 x, vec3 y)
+{
+    float x0 = x.x * y.x;
+    float p = x.x * y.y;
+    float q = y.x * x.y;
+    float a = p + q;
+    float e = abs(x.x) * y.z + abs(y.x) * x.z + (abs(x.y) + x.z) * (abs(y.y) + y.z);
+    return vec3(x0, a, e + (abs(x0) + abs(p) + abs(q) + e) * 9.5367431640625e-7 + 1.17549435e-38);
+}
+
+// The narrower of the Lipschitz form and the centred (Taylor) form over every p = centre + axis eps + w,
+// |w| <= radius, from one cultmath_snoise_grad at the centre; the allowance for the float32 evaluation is
+// in math.Affine.cs.
+vec3 cultmath_af_snoise(vec3 centre, vec3 axis, float radius)
+{
+    vec4 g = cultmath_snoise_grad(centre);
+    vec3 gradient = vec3(g.x, g.y, g.z);
+    float reach = length(axis) + radius;
+    float lipschitz = CULTMATH_SNOISE_LIPSCHITZ * reach;
+    float a = dot(gradient, axis);
+    float remainder = length(gradient) * radius + 0.5 * CULTMATH_SNOISE_HESSIAN * reach * reach;
+    float allowance = 6.103515625e-5 + 2.44140625e-4 * reach
+        + CULTMATH_SNOISE_LIPSCHITZ * 9.5367431640625e-7 * (abs(centre.x) + abs(centre.y) + abs(centre.z) + reach);
+    bool centred = abs(a) + remainder < lipschitz;
+    float slope = centred ? a : 0.0;
+    float e = (centred ? remainder : lipschitz) + allowance;
+    return vec3(g.w, slope, e + (abs(g.w) + abs(slope) + e) * 9.5367431640625e-7);
+}
+
+// Octave sum of cultmath_af_snoise, compounding frequency and amplitude as the value-and-gradient fBm does.
+vec3 cultmath_af_fbm(vec3 centre, vec3 axis, float radius, int octaves, float lacunarity, float gain)
+{
+    octaves = clamp(octaves, 0, 16);
+    float amplitude = 1.0;
+    float frequency = 1.0;
+    vec3 sum = vec3(0.0, 0.0, 0.0);
+    for (int i = 0; i < octaves; i++)
+    {
+        sum = cultmath_af_add(sum, cultmath_af_scale(cultmath_af_snoise(centre * frequency, axis * frequency, radius * abs(frequency)), amplitude));
+        frequency *= lacunarity;
+        amplitude *= gain;
+    }
+
+    return sum;
+}
+
+// The axis of a screen tile's slice over the depth segment [z0, z1]; camera frame of
+// cultmath_iv_frustum_ball.
+vec3 cultmath_af_frustum_axis(vec2 centreSlope, float z0, float z1)
+{
+    float h = (z1 - z0) * 0.5;
+    return vec3(centreSlope.x * h, centreSlope.y * h, h);
+}
+
+// The affine slice of a screen tile over [z0, z1]: centre xyz, radius w. Derivation in math.Affine.cs.
+vec4 cultmath_af_frustum_ball(vec2 centreSlope, float z0, float z1, float footprintPerDepth, float warpVariation)
+{
+    float zm = (z0 + z1) * 0.5;
+    float h = (z1 - z0) * 0.5;
+    float cx = centreSlope.x * zm;
+    float cy = centreSlope.y * zm;
+    float spread = (abs(centreSlope.x) + abs(centreSlope.y) + 1.0) * h;
+    float radius = z1 * footprintPerDepth + warpVariation;
+    radius += (abs(cx) + abs(cy) + abs(zm) + spread + radius) * 9.5367431640625e-7;
+    return vec4(cx, cy, zm, radius);
+}
+
+
+
+
 float cultmath_value_noise(vec2 position)
 {
     vec2 cell = floor(position);

packages/cultmath/shaders/CultMath.hlsl

@@ -509,6 +509,7 @@ CultCellular cultmath_cellular(float3 p)
 
 #include "CultMath.Phacelle.hlsl"
 #include "CultMath.Interval.hlsl"
+#include "CultMath.Affine.hlsl"
 
 float cultmath_value_noise(float2 position)
 {

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

@@ -66,6 +66,7 @@ public sealed class HlslMirrorTests
     [InlineData("CultMath.hlsl")]
     [InlineData("CultMath.Phacelle.hlsl")]
     [InlineData("CultMath.Interval.hlsl")]
+    [InlineData("CultMath.Affine.hlsl")]
     public void UnityPackageShaderIsIdenticalToTheCanonicalShader(string name)
     {
         var root = GetCultMathRoot();

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

@@ -33,6 +33,10 @@ public sealed class HlslSourceCompatibilityTests
         Assert.NotNull(ShaderFunction(shader, "cultmath_phacelle", typeof(float3), typeof(float3), typeof(float), typeof(float)));
         Assert.NotNull(ShaderFunction(shader, "cultmath_iv_snoise_ball", typeof(float3), typeof(float)));
         Assert.NotNull(ShaderFunction(shader, "cultmath_iv_fbm_ball", typeof(float3), typeof(float), typeof(int), typeof(float), typeof(float)));
+        Assert.NotNull(ShaderFunction(shader, "cultmath_iv_snoise_ball", typeof(float2), typeof(float)));
+        Assert.NotNull(ShaderFunction(shader, "cultmath_af_snoise", typeof(float3), typeof(float3), typeof(float)));
+        Assert.NotNull(ShaderFunction(shader, "cultmath_af_fbm", typeof(float3), typeof(float3), typeof(float), typeof(int), typeof(float), typeof(float)));
+        Assert.NotNull(ShaderFunction(shader, "cultmath_af_frustum_ball", typeof(float2), typeof(float), typeof(float), typeof(float), typeof(float)));
     }
 
     // HLSL functions carry no access modifier, so they compile as private instance methods.
@@ -133,7 +137,7 @@ public sealed class HlslSourceCompatibilityTests
         }
 
         Assert.True(mismatches.Count == 0, $"{mismatches.Count} mismatches:{Environment.NewLine}{string.Join(Environment.NewLine, mismatches.Take(12))}");
-        Assert.Equal(53, compared.Count);
+        Assert.Equal(67, compared.Count);
     }
 
     /// <summary>

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

@@ -0,0 +1,130 @@
+#ifndef CULTMATH_AFFINE_HLSL
+#define CULTMATH_AFFINE_HLSL
+
+// Reduced affine forms with one shared symbol over float3(x0, a, e): the set { x0 + a eps + e delta :
+// delta in [-1, 1] } at the region's shared eps in [-1, 1], e >= 0. Every cultmath_af_* function mirrors
+// math.Affine.cs bit for bit; see that file and docs/design.md ("Affine forms") for each op's enclosure
+// and the float32 rounding allowance it folds into e. Included by CultMath.hlsl after
+// CultMath.Interval.hlsl, whose CULTMATH_SNOISE_LIPSCHITZ it shares, and after its own
+// definition of cultmath_snoise_grad.
+
+// An upper bound on the spectral norm of the Hessian of cultmath_snoise(float3); provenance in
+// math.Affine.cs (SNOISE_HESSIAN).
+static const float CULTMATH_SNOISE_HESSIAN = 56.050385;
+
+float3 cultmath_af_point(float value) { return float3(value, 0.0, 0.0); }
+float3 cultmath_af_symbol(float x0, float a) { return float3(x0, a, 0.0); }
+
+float3 cultmath_af_from_iv(float2 interval)
+{
+    float x0 = interval.x * 0.5 + interval.y * 0.5;
+    float e = interval.y * 0.5 - interval.x * 0.5;
+    return float3(x0, 0.0, e + (abs(x0) + e) * 9.5367431640625e-7 + 1.17549435e-38);
+}
+
+// An endpoint that would be NaN is the infinity on its side.
+float2 cultmath_af_range(float3 x)
+{
+    float r = abs(x.y) + x.z;
+    float w = abs(x.x) * 9.5367431640625e-7 + r * 9.5367431640625e-7;
+    float lo = x.x - r - w;
+    float hi = x.x + r + w;
+    return float2(lo == lo ? lo : asfloat(0xff800000u), hi == hi ? hi : asfloat(0x7f800000u));
+}
+
+float3 cultmath_af_add(float3 x, float3 y)
+{
+    float x0 = x.x + y.x;
+    float a = x.y + y.y;
+    float e = x.z + y.z;
+    return float3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7);
+}
+
+float3 cultmath_af_sub(float3 x, float3 y)
+{
+    float x0 = x.x - y.x;
+    float a = x.y - y.y;
+    float e = x.z + y.z;
+    return float3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7);
+}
+
+float3 cultmath_af_neg(float3 x) { return float3(-x.x, -x.y, x.z); }
+
+float3 cultmath_af_scale(float3 x, float s)
+{
+    float x0 = x.x * s;
+    float a = x.y * s;
+    float e = x.z * abs(s);
+    return float3(x0, a, e + (abs(x0) + abs(a) + e) * 9.5367431640625e-7 + 1.17549435e-38);
+}
+
+float3 cultmath_af_add_iv(float3 x, float2 interval) { return cultmath_af_add(x, cultmath_af_from_iv(interval)); }
+
+float3 cultmath_af_mul(float3 x, float3 y)
+{
+    float x0 = x.x * y.x;
+    float p = x.x * y.y;
+    float q = y.x * x.y;
+    float a = p + q;
+    float e = abs(x.x) * y.z + abs(y.x) * x.z + (abs(x.y) + x.z) * (abs(y.y) + y.z);
+    return float3(x0, a, e + (abs(x0) + abs(p) + abs(q) + e) * 9.5367431640625e-7 + 1.17549435e-38);
+}
+
+// The narrower of the Lipschitz form and the centred (Taylor) form over every p = centre + axis eps + w,
+// |w| <= radius, from one cultmath_snoise_grad at the centre; the allowance for the float32 evaluation is
+// in math.Affine.cs.
+float3 cultmath_af_snoise(float3 centre, float3 axis, float radius)
+{
+    float4 g = cultmath_snoise_grad(centre);
+    float3 gradient = float3(g.x, g.y, g.z);
+    float reach = length(axis) + radius;
+    float lipschitz = CULTMATH_SNOISE_LIPSCHITZ * reach;
+    float a = dot(gradient, axis);
+    float remainder = length(gradient) * radius + 0.5 * CULTMATH_SNOISE_HESSIAN * reach * reach;
+    float allowance = 6.103515625e-5 + 2.44140625e-4 * reach
+        + CULTMATH_SNOISE_LIPSCHITZ * 9.5367431640625e-7 * (abs(centre.x) + abs(centre.y) + abs(centre.z) + reach);
+    bool centred = abs(a) + remainder < lipschitz;
+    float slope = centred ? a : 0.0;
+    float e = (centred ? remainder : lipschitz) + allowance;
+    return float3(g.w, slope, e + (abs(g.w) + abs(slope) + e) * 9.5367431640625e-7);
+}
+
+// Octave sum of cultmath_af_snoise, compounding frequency and amplitude as the value-and-gradient fBm does.
+float3 cultmath_af_fbm(float3 centre, float3 axis, float radius, int octaves, float lacunarity, float gain)
+{
+    octaves = clamp(octaves, 0, 16);
+    float amplitude = 1.0;
+    float frequency = 1.0;
+    float3 sum = float3(0.0, 0.0, 0.0);
+    for (int i = 0; i < octaves; i++)
+    {
+        sum = cultmath_af_add(sum, cultmath_af_scale(cultmath_af_snoise(centre * frequency, axis * frequency, radius * abs(frequency)), amplitude));
+        frequency *= lacunarity;
+        amplitude *= gain;
+    }
+
+    return sum;
+}
+
+// The axis of a screen tile's slice over the depth segment [z0, z1]; camera frame of
+// cultmath_iv_frustum_ball.
+float3 cultmath_af_frustum_axis(float2 centreSlope, float z0, float z1)
+{
+    float h = (z1 - z0) * 0.5;
+    return float3(centreSlope.x * h, centreSlope.y * h, h);
+}
+
+// The affine slice of a screen tile over [z0, z1]: centre xyz, radius w. Derivation in math.Affine.cs.
+float4 cultmath_af_frustum_ball(float2 centreSlope, float z0, float z1, float footprintPerDepth, float warpVariation)
+{
+    float zm = (z0 + z1) * 0.5;
+    float h = (z1 - z0) * 0.5;
+    float cx = centreSlope.x * zm;
+    float cy = centreSlope.y * zm;
+    float spread = (abs(centreSlope.x) + abs(centreSlope.y) + 1.0) * h;
+    float radius = z1 * footprintPerDepth + warpVariation;
+    radius += (abs(cx) + abs(cy) + abs(zm) + spread + radius) * 9.5367431640625e-7;
+    return float4(cx, cy, zm, radius);
+}
+
+#endif

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

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

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

@@ -9,6 +9,9 @@
 // Lipschitz constant of cultmath_snoise(float3); provenance in math.Interval.cs (SNOISE_LIPSCHITZ).
 static const float CULTMATH_SNOISE_LIPSCHITZ = 7.9640074;
 
+// Lipschitz constant of cultmath_snoise(float2); provenance in math.Interval.cs (SNOISE2_LIPSCHITZ).
+static const float CULTMATH_SNOISE2_LIPSCHITZ = 8.117002;
+
 float2 cultmath_iv_point(float value) { return float2(value, value); }
 float2 cultmath_iv_add(float2 a, float2 b) { return float2(a.x + b.x, a.y + b.y); }
 float2 cultmath_iv_sub(float2 a, float2 b) { return float2(a.x - b.y, a.y - b.x); }
@@ -87,6 +90,14 @@ float2 cultmath_iv_snoise_ball(float3 centre, float radius)
     return float2(max(n - e, -1.0), min(n + e, 1.0));
 }
 
+// The 2D overload: [n - L2 r, n + L2 r] intersected with [-1, 1], n = cultmath_snoise(centre). One snoise, value only.
+float2 cultmath_iv_snoise_ball(float2 centre, float radius)
+{
+    float n = cultmath_snoise(centre);
+    float e = CULTMATH_SNOISE2_LIPSCHITZ * radius;
+    return float2(max(n - e, -1.0), min(n + e, 1.0));
+}
+
 // Octave sum of cultmath_iv_snoise_ball, compounding frequency and amplitude as the value-and-gradient fBm does.
 float2 cultmath_iv_fbm_ball(float3 centre, float radius, int octaves, float lacunarity, float gain)
 {

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

@@ -509,6 +509,7 @@ CultCellular cultmath_cellular(float3 p)
 
 #include "CultMath.Phacelle.hlsl"
 #include "CultMath.Interval.hlsl"
+#include "CultMath.Affine.hlsl"
 
 float cultmath_value_noise(float2 position)
 {