Promote shader math primitives

README.md

@@ -17,7 +17,13 @@ shader semantics.
 - conversions to and from `System.Numerics`
 - `math` intrinsics: `radians`, `degrees`, `abs`, `floor`, `ceil`, `frac`,
   `min`, `max`, `clamp`, `saturate`, `lerp`, `step`, `smoothstep`, `dot`,
-  `cross`, `length`, `distance`, `normalize`, and `reflect`
+  `cross`, `length`, `distance`, `normalize`, `reflect`, `csum`, exponential
+  `decay`/`damp`, intercept/segment helpers, Catmull-Rom, Bezier curves, and
+  deterministic value-noise primitives
+- `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, including CultMath's
+  safe `normalize` contract.
 - `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.
@@ -38,8 +44,25 @@ using CultMath;
 float3 normal = math.normalize(new float3(0.25f, 1.0f, -0.1f));
 float rim = math.saturate(1.0f - math.dot(normal, new float3(0.0f, 0.0f, -1.0f)));
 float glow = math.smoothstep(0.2f, 1.0f, rim);
+float grain = math.value_noise_bicubic(math.float2(12.0f, 3.5f));
 ```
 
+```hlsl
+#include "CultMath.hlsl"
+
+float3 normal = cultmath_normalize(float3(0.25, 1.0, -0.1));
+float rim = cultmath_saturate(1.0 - dot(normal, float3(0.0, 0.0, -1.0)));
+float glow = cultmath_smoothstep(0.2, 1.0, rim);
+float grain = cultmath_value_noise_bicubic(float2(12.0, 3.5));
+```
+
+## Provenance
+
+The small spline, damping, and pursuit helpers were promoted from Aetheria's
+local `AetheriaMath` surface. The value-noise family follows public shader
+literature from Inigo Quilez's noise and mini-spline articles, expressed in
+CultMath names with tests instead of copied project-local helper drift.
+
 ## Build
 
 ```powershell
@@ -55,3 +78,8 @@ The goal is mechanical sympathy between CPU-side authoring, native kernels, and
 GPU-side shader math. If an intrinsic exists here, it should behave like the
 shader concept unless a runtime's numeric rules make that impossible. Weird
 deviations get tests and documentation, not folklore.
+
+CultMath is still not a shader generator. Projects should include
+`shaders/CultMath.hlsl` directly or add that directory to their shader include
+path. If a shared operation needs CPU/GPU parity, put the primitive here and
+test it here before copying local math helpers into a renderer or daemon.

docs/design.md

@@ -9,16 +9,21 @@ Rust hot kernels, C# host code, and HLSL shader code.
 - Keep Rust as the owner of native hot-path kernels and parity fixtures.
 - Keep C#, shader, and scripting surfaces as runtime wrappers, not competing
   math authorities.
+- Keep `shaders/CultMath.hlsl` as the canonical shader mirror for shared
+  intrinsics. Project-local shader helpers should move here once more than one
+  organ needs them or once CPU/GPU parity matters.
 - Implement public shader semantics directly; do not copy Unity.Mathematics.
 - Prefer small immutable value types.
 - Prefer component-wise overloads over clever generic machinery until the API
   earns that complexity.
-- Keep the first surface boring and exact. Noise, matrices, quaternions, and
-  packing helpers come later.
+- Keep the surface boring and exact. Noise primitives are allowed when they are
+  source-grounded, shader-shaped, and shared by more than one project. Matrices,
+  quaternions, and packing helpers still need an owner before they enter.
 
 ## Non-Goals
 
 - Compiling C# to shaders.
+- Generating HLSL from Rust or Rust from HLSL.
 - Replacing `System.Numerics` for general .NET work.
 - Recreating every Unity.Mathematics type because a spreadsheet somewhere says
   “coverage.” That way lies decorative obesity.

shaders/CultMath.hlsl

@@ -0,0 +1,254 @@
+#ifndef CULTMATH_HLSL
+#define CULTMATH_HLSL
+
+// CultMath HLSL mirror.
+//
+// HLSL already owns the float2/float3/float4 type names and most primitive
+// intrinsics. This include names the CultMath semantic surface explicitly so
+// shader code, C# wrappers, and Rust/native parity fixtures can share the same
+// small vocabulary without local helper drift.
+
+static const float CULTMATH_PI = 3.14159265358979323846;
+static const float CULTMATH_TAU = 6.28318530717958647692;
+static const float CULTMATH_HALF_PI = 1.57079632679489661923;
+static const float CULTMATH_NORMALIZE_EPSILON = 1.0e-20;
+
+float cultmath_radians(float degrees_value) { return degrees_value * (CULTMATH_PI / 180.0); }
+float cultmath_degrees(float radians_value) { return radians_value * (180.0 / CULTMATH_PI); }
+
+float cultmath_frac(float value) { return value - floor(value); }
+float2 cultmath_frac(float2 value) { return value - floor(value); }
+float3 cultmath_frac(float3 value) { return value - floor(value); }
+float4 cultmath_frac(float4 value) { return value - floor(value); }
+
+float cultmath_clamp(float value, float minimum, float maximum) { return min(max(value, minimum), maximum); }
+float2 cultmath_clamp(float2 value, float2 minimum, float2 maximum) { return min(max(value, minimum), maximum); }
+float3 cultmath_clamp(float3 value, float3 minimum, float3 maximum) { return min(max(value, minimum), maximum); }
+float4 cultmath_clamp(float4 value, float4 minimum, float4 maximum) { return min(max(value, minimum), maximum); }
+
+float cultmath_saturate(float value) { return cultmath_clamp(value, 0.0, 1.0); }
+float2 cultmath_saturate(float2 value) { return cultmath_clamp(value, float2(0.0, 0.0), float2(1.0, 1.0)); }
+float3 cultmath_saturate(float3 value) { return cultmath_clamp(value, float3(0.0, 0.0, 0.0), float3(1.0, 1.0, 1.0)); }
+float4 cultmath_saturate(float4 value) { return cultmath_clamp(value, float4(0.0, 0.0, 0.0, 0.0), float4(1.0, 1.0, 1.0, 1.0)); }
+
+float cultmath_lerp(float start, float end, float amount) { return start + (end - start) * amount; }
+float2 cultmath_lerp(float2 start, float2 end, float2 amount) { return start + (end - start) * amount; }
+float3 cultmath_lerp(float3 start, float3 end, float3 amount) { return start + (end - start) * amount; }
+float4 cultmath_lerp(float4 start, float4 end, float4 amount) { return start + (end - start) * amount; }
+
+float cultmath_step(float edge, float value) { return value < edge ? 0.0 : 1.0; }
+float2 cultmath_step(float2 edge, float2 value) { return float2(cultmath_step(edge.x, value.x), cultmath_step(edge.y, value.y)); }
+float3 cultmath_step(float3 edge, float3 value) { return float3(cultmath_step(edge.x, value.x), cultmath_step(edge.y, value.y), cultmath_step(edge.z, value.z)); }
+float4 cultmath_step(float4 edge, float4 value) { return float4(cultmath_step(edge.x, value.x), cultmath_step(edge.y, value.y), cultmath_step(edge.z, value.z), cultmath_step(edge.w, value.w)); }
+
+float cultmath_smoothstep(float minimum, float maximum, float value)
+{
+    float t = cultmath_saturate((value - minimum) / (maximum - minimum));
+    return t * t * (3.0 - 2.0 * t);
+}
+
+float2 cultmath_smoothstep(float2 minimum, float2 maximum, float2 value)
+{
+    float2 t = cultmath_saturate((value - minimum) / (maximum - minimum));
+    return t * t * (float2(3.0, 3.0) - float2(2.0, 2.0) * t);
+}
+
+float3 cultmath_smoothstep(float3 minimum, float3 maximum, float3 value)
+{
+    float3 t = cultmath_saturate((value - minimum) / (maximum - minimum));
+    return t * t * (float3(3.0, 3.0, 3.0) - float3(2.0, 2.0, 2.0) * t);
+}
+
+float4 cultmath_smoothstep(float4 minimum, float4 maximum, float4 value)
+{
+    float4 t = cultmath_saturate((value - minimum) / (maximum - minimum));
+    return t * t * (float4(3.0, 3.0, 3.0, 3.0) - float4(2.0, 2.0, 2.0, 2.0) * t);
+}
+
+float cultmath_smoothstep01(float value) { return cultmath_smoothstep(0.0, 1.0, value); }
+
+float cultmath_smootherstep(float value)
+{
+    value = cultmath_saturate(value);
+    return value * value * value * (value * (value * 6.0 - 15.0) + 10.0);
+}
+
+float cultmath_lengthsq(float2 value) { return dot(value, value); }
+float cultmath_lengthsq(float3 value) { return dot(value, value); }
+float cultmath_lengthsq(float4 value) { return dot(value, value); }
+
+float cultmath_distance(float2 left, float2 right) { return length(left - right); }
+float cultmath_distance(float3 left, float3 right) { return length(left - right); }
+float cultmath_distance(float4 left, float4 right) { return length(left - right); }
+
+float2 cultmath_normalize(float2 value) { return value / max(length(value), CULTMATH_NORMALIZE_EPSILON); }
+float3 cultmath_normalize(float3 value) { return value / max(length(value), CULTMATH_NORMALIZE_EPSILON); }
+float4 cultmath_normalize(float4 value) { return value / max(length(value), CULTMATH_NORMALIZE_EPSILON); }
+
+float2 cultmath_reflect(float2 incident, float2 normal_value)
+{
+    return incident - 2.0 * dot(normal_value, incident) * normal_value;
+}
+
+float3 cultmath_reflect(float3 incident, float3 normal_value)
+{
+    return incident - 2.0 * dot(normal_value, incident) * normal_value;
+}
+
+float cultmath_csum(float2 value) { return value.x + value.y; }
+float cultmath_csum(float3 value) { return value.x + value.y + value.z; }
+float cultmath_csum(float4 value) { return value.x + value.y + value.z + value.w; }
+
+float cultmath_decay(float source, float lambda_value, float dt) { return source * exp(-lambda_value * dt); }
+float2 cultmath_decay(float2 source, float lambda_value, float dt) { return source * exp(-lambda_value * dt); }
+float3 cultmath_decay(float3 source, float lambda_value, float dt) { return source * exp(-lambda_value * dt); }
+float4 cultmath_decay(float4 source, float lambda_value, float dt) { return source * exp(-lambda_value * dt); }
+
+float cultmath_damp(float start, float end, float lambda_value, float dt)
+{
+    return cultmath_lerp(start, end, 1.0 - exp(-lambda_value * dt));
+}
+
+float2 cultmath_damp(float2 start, float2 end, float lambda_value, float dt)
+{
+    return cultmath_lerp(start, end, 1.0 - exp(-lambda_value * dt));
+}
+
+float3 cultmath_damp(float3 start, float3 end, float lambda_value, float dt)
+{
+    return cultmath_lerp(start, end, 1.0 - exp(-lambda_value * dt));
+}
+
+float4 cultmath_damp(float4 start, float4 end, float lambda_value, float dt)
+{
+    return cultmath_lerp(start, end, 1.0 - exp(-lambda_value * dt));
+}
+
+float cultmath_catmullrom(float p0, float p1, float p2, float p3, float t)
+{
+    float t2 = t * t;
+    float t3 = t2 * t;
+    return 0.5 * ((2.0 * p1) + (p2 - p0) * t + (2.0 * p0 - 5.0 * p1 + 4.0 * p2 - p3) * t2 + (-p0 + 3.0 * p1 - 3.0 * p2 + p3) * t3);
+}
+
+float2 cultmath_catmullrom(float2 p0, float2 p1, float2 p2, float2 p3, float t)
+{
+    float t2 = t * t;
+    float t3 = t2 * t;
+    return 0.5 * ((2.0 * p1) + (p2 - p0) * t + (2.0 * p0 - 5.0 * p1 + 4.0 * p2 - p3) * t2 + (-p0 + 3.0 * p1 - 3.0 * p2 + p3) * t3);
+}
+
+float3 cultmath_catmullrom(float3 p0, float3 p1, float3 p2, float3 p3, float t)
+{
+    float t2 = t * t;
+    float t3 = t2 * t;
+    return 0.5 * ((2.0 * p1) + (p2 - p0) * t + (2.0 * p0 - 5.0 * p1 + 4.0 * p2 - p3) * t2 + (-p0 + 3.0 * p1 - 3.0 * p2 + p3) * t3);
+}
+
+float4 cultmath_catmullrom(float4 p0, float4 p1, float4 p2, float4 p3, float t)
+{
+    float t2 = t * t;
+    float t3 = t2 * t;
+    return 0.5 * ((2.0 * p1) + (p2 - p0) * t + (2.0 * p0 - 5.0 * p1 + 4.0 * p2 - p3) * t2 + (-p0 + 3.0 * p1 - 3.0 * p2 + p3) * t3);
+}
+
+float cultmath_quadratic_bezier(float p0, float p1, float p2, float t)
+{
+    t = cultmath_saturate(t);
+    float inv = 1.0 - t;
+    return inv * inv * p0 + 2.0 * inv * t * p1 + t * t * p2;
+}
+
+float3 cultmath_quadratic_bezier(float3 p0, float3 p1, float3 p2, float t)
+{
+    t = cultmath_saturate(t);
+    float inv = 1.0 - t;
+    return inv * inv * p0 + 2.0 * inv * t * p1 + t * t * p2;
+}
+
+float cultmath_cubic_bezier(float p0, float p1, float p2, float p3, float t)
+{
+    t = cultmath_saturate(t);
+    float inv = 1.0 - t;
+    return inv * inv * inv * p0 + 3.0 * inv * inv * t * p1 + 3.0 * inv * t * t * p2 + t * t * t * p3;
+}
+
+float3 cultmath_cubic_bezier(float3 p0, float3 p1, float3 p2, float3 p3, float t)
+{
+    t = cultmath_saturate(t);
+    float inv = 1.0 - t;
+    return inv * inv * inv * p0 + 3.0 * inv * inv * t * p1 + 3.0 * inv * t * t * p2 + t * t * t * p3;
+}
+
+float cultmath_hash(float value) { return cultmath_frac(sin(value) * 43758.5453); }
+float cultmath_hash(float2 value) { return cultmath_hash(dot(value, float2(127.1, 311.7))); }
+float cultmath_hash(float3 value) { return cultmath_hash(dot(value, float3(127.1, 311.7, 74.7))); }
+
+float cultmath_value_noise(float2 position)
+{
+    float2 cell = floor(position);
+    float2 local = cultmath_frac(position);
+    float2 u = local * local * (float2(3.0, 3.0) - float2(2.0, 2.0) * local);
+
+    float a = cultmath_hash(cell);
+    float b = cultmath_hash(cell + float2(1.0, 0.0));
+    float c = cultmath_hash(cell + float2(0.0, 1.0));
+    float d = cultmath_hash(cell + float2(1.0, 1.0));
+
+    return cultmath_lerp(cultmath_lerp(a, b, u.x), cultmath_lerp(c, d, u.x), u.y);
+}
+
+float cultmath_value_noise_bicubic(float2 position)
+{
+    float2 cell = floor(position);
+    float2 local = cultmath_frac(position);
+
+    float y0 = cultmath_catmullrom(cultmath_hash(cell + float2(-1.0, -1.0)), cultmath_hash(cell + float2(0.0, -1.0)), cultmath_hash(cell + float2(1.0, -1.0)), cultmath_hash(cell + float2(2.0, -1.0)), local.x);
+    float y1 = cultmath_catmullrom(cultmath_hash(cell + float2(-1.0, 0.0)), cultmath_hash(cell + float2(0.0, 0.0)), cultmath_hash(cell + float2(1.0, 0.0)), cultmath_hash(cell + float2(2.0, 0.0)), local.x);
+    float y2 = cultmath_catmullrom(cultmath_hash(cell + float2(-1.0, 1.0)), cultmath_hash(cell + float2(0.0, 1.0)), cultmath_hash(cell + float2(1.0, 1.0)), cultmath_hash(cell + float2(2.0, 1.0)), local.x);
+    float y3 = cultmath_catmullrom(cultmath_hash(cell + float2(-1.0, 2.0)), cultmath_hash(cell + float2(0.0, 2.0)), cultmath_hash(cell + float2(1.0, 2.0)), cultmath_hash(cell + float2(2.0, 2.0)), local.x);
+
+    return cultmath_catmullrom(y0, y1, y2, y3, local.y);
+}
+
+float cultmath_value_noise_texture(Texture2D<float> noise_texture, SamplerState noise_sampler, float2 uv, float scale)
+{
+    return noise_texture.SampleLevel(noise_sampler, uv * scale, 0.0).r;
+}
+
+float cultmath_value_noise_texture_bicubic(Texture2D<float> noise_texture, SamplerState noise_sampler, float2 uv, float scale, float2 texel_size)
+{
+    float2 position = uv * scale;
+    float2 cell = floor(position);
+    float2 local = cultmath_frac(position);
+    float2 base_uv = cell / scale;
+
+    float row0 = cultmath_catmullrom(
+        noise_texture.SampleLevel(noise_sampler, base_uv + texel_size * float2(-1.0, -1.0), 0.0).r,
+        noise_texture.SampleLevel(noise_sampler, base_uv + texel_size * float2(0.0, -1.0), 0.0).r,
+        noise_texture.SampleLevel(noise_sampler, base_uv + texel_size * float2(1.0, -1.0), 0.0).r,
+        noise_texture.SampleLevel(noise_sampler, base_uv + texel_size * float2(2.0, -1.0), 0.0).r,
+        local.x);
+    float row1 = cultmath_catmullrom(
+        noise_texture.SampleLevel(noise_sampler, base_uv + texel_size * float2(-1.0, 0.0), 0.0).r,
+        noise_texture.SampleLevel(noise_sampler, base_uv + texel_size * float2(0.0, 0.0), 0.0).r,
+        noise_texture.SampleLevel(noise_sampler, base_uv + texel_size * float2(1.0, 0.0), 0.0).r,
+        noise_texture.SampleLevel(noise_sampler, base_uv + texel_size * float2(2.0, 0.0), 0.0).r,
+        local.x);
+    float row2 = cultmath_catmullrom(
+        noise_texture.SampleLevel(noise_sampler, base_uv + texel_size * float2(-1.0, 1.0), 0.0).r,
+        noise_texture.SampleLevel(noise_sampler, base_uv + texel_size * float2(0.0, 1.0), 0.0).r,
+        noise_texture.SampleLevel(noise_sampler, base_uv + texel_size * float2(1.0, 1.0), 0.0).r,
+        noise_texture.SampleLevel(noise_sampler, base_uv + texel_size * float2(2.0, 1.0), 0.0).r,
+        local.x);
+    float row3 = cultmath_catmullrom(
+        noise_texture.SampleLevel(noise_sampler, base_uv + texel_size * float2(-1.0, 2.0), 0.0).r,
+        noise_texture.SampleLevel(noise_sampler, base_uv + texel_size * float2(0.0, 2.0), 0.0).r,
+        noise_texture.SampleLevel(noise_sampler, base_uv + texel_size * float2(1.0, 2.0), 0.0).r,
+        noise_texture.SampleLevel(noise_sampler, base_uv + texel_size * float2(2.0, 2.0), 0.0).r,
+        local.x);
+
+    return cultmath_catmullrom(row0, row1, row2, row3, local.y);
+}
+
+#endif

src/CultMath/CultMath.csproj

@@ -17,6 +17,7 @@
 
   <ItemGroup>
     <None Include="..\..\README.md" Pack="true" PackagePath="\" />
+    <None Include="..\..\shaders\CultMath.hlsl" Pack="true" PackagePath="contentFiles\any\any\shaders\CultMath.hlsl" CopyToOutputDirectory="PreserveNewest" />
   </ItemGroup>
 
 </Project>

src/CultMath/math.cs

@@ -6,6 +6,16 @@ public static class math
     public const float TAU = MathF.PI * 2.0f;
     public const float HALF_PI = MathF.PI * 0.5f;
 
+    public static float2 float2(float x, float y) => new(x, y);
+    public static float2 float2(float value) => new(value, value);
+    public static float3 float3(float x, float y, float z) => new(x, y, z);
+    public static float3 float3(float2 xy, float z) => new(xy, z);
+    public static float3 float3(float value) => new(value, value, value);
+    public static float4 float4(float x, float y, float z, float w) => new(x, y, z, w);
+    public static float4 float4(float2 xy, float z, float w) => new(xy, z, w);
+    public static float4 float4(float3 xyz, float w) => new(xyz, w);
+    public static float4 float4(float value) => new(value, value, value, value);
+
     public static float radians(float degrees) => degrees * (PI / 180.0f);
     public static float degrees(float radians) => radians * (180.0f / PI);
     public static float sin(float value) => MathF.Sin(value);
@@ -16,6 +26,10 @@ public static class math
     public static float atan(float value) => MathF.Atan(value);
     public static float atan2(float y, float x) => MathF.Atan2(y, x);
     public static float sqrt(float value) => MathF.Sqrt(value);
+    public static float exp(float value) => MathF.Exp(value);
+    public static float2 exp(float2 value) => new(exp(value.x), exp(value.y));
+    public static float3 exp(float3 value) => new(exp(value.x), exp(value.y), exp(value.z));
+    public static float4 exp(float4 value) => new(exp(value.x), exp(value.y), exp(value.z), exp(value.w));
 
     public static float abs(float value) => MathF.Abs(value);
     public static float2 abs(float2 value) => new(abs(value.x), abs(value.y));
@@ -119,4 +133,192 @@ public static class math
 
     public static float2 reflect(float2 incident, float2 normal) => incident - 2.0f * dot(normal, incident) * normal;
     public static float3 reflect(float3 incident, float3 normal) => incident - 2.0f * dot(normal, incident) * normal;
+
+    public static float csum(float2 value) => value.x + value.y;
+    public static float csum(float3 value) => value.x + value.y + value.z;
+    public static float csum(float4 value) => value.x + value.y + value.z + value.w;
+
+    public static float decay(float source, float lambda, float dt) => source * exp(-lambda * dt);
+    public static float2 decay(float2 source, float lambda, float dt) => source * exp(-lambda * dt);
+    public static float3 decay(float3 source, float lambda, float dt) => source * exp(-lambda * dt);
+    public static float4 decay(float4 source, float lambda, float dt) => source * exp(-lambda * dt);
+    public static float2 decay(float2 source, float2 lambda, float dt) => source * exp(-lambda * dt);
+    public static float3 decay(float3 source, float3 lambda, float dt) => source * exp(-lambda * dt);
+    public static float4 decay(float4 source, float4 lambda, float dt) => source * exp(-lambda * dt);
+
+    public static float damp(float start, float end, float lambda, float dt) => lerp(start, end, 1.0f - exp(-lambda * dt));
+    public static float2 damp(float2 start, float2 end, float lambda, float dt) => lerp(start, end, 1.0f - exp(-lambda * dt));
+    public static float3 damp(float3 start, float3 end, float lambda, float dt) => lerp(start, end, 1.0f - exp(-lambda * dt));
+    public static float4 damp(float4 start, float4 end, float lambda, float dt) => lerp(start, end, 1.0f - exp(-lambda * dt));
+
+    public static float first_order_intercept_time(float shotSpeed, float3 targetRelativePosition, float3 targetRelativeVelocity)
+    {
+        var velocitySquared = lengthsq(targetRelativeVelocity);
+        if (velocitySquared < 0.001f)
+        {
+            return 0.0f;
+        }
+
+        var a = velocitySquared - shotSpeed * shotSpeed;
+        if (abs(a) < 0.001f)
+        {
+            var time = -lengthsq(targetRelativePosition) / (2.0f * dot(targetRelativeVelocity, targetRelativePosition));
+            return max(time, 0.0f);
+        }
+
+        var b = 2.0f * dot(targetRelativeVelocity, targetRelativePosition);
+        var c = lengthsq(targetRelativePosition);
+        var determinant = b * b - 4.0f * a * c;
+
+        if (determinant > 0.0f)
+        {
+            var root = sqrt(determinant);
+            var t1 = (-b + root) / (2.0f * a);
+            var t2 = (-b - root) / (2.0f * a);
+            if (t1 > 0.0f)
+            {
+                return t2 > 0.0f ? min(t1, t2) : t1;
+            }
+
+            return max(t2, 0.0f);
+        }
+
+        if (determinant < 0.0f)
+        {
+            return 0.0f;
+        }
+
+        return max(-b / (2.0f * a), 0.0f);
+    }
+
+    public static float3 first_order_intercept(
+        float3 shooterPosition,
+        float3 shooterVelocity,
+        float shotSpeed,
+        float3 targetPosition,
+        float3 targetVelocity)
+    {
+        var targetRelativePosition = targetPosition - shooterPosition;
+        var targetRelativeVelocity = targetVelocity - shooterVelocity;
+        var time = first_order_intercept_time(shotSpeed, targetRelativePosition, targetRelativeVelocity);
+        return targetPosition + time * targetRelativeVelocity;
+    }
+
+    public static float distance_to_segment(float2 point, float2 start, float2 end, out float2 closest)
+    {
+        var segment = end - start;
+        var segmentLengthSq = lengthsq(segment);
+        if (segmentLengthSq <= 0.0f)
+        {
+            closest = start;
+            return distance(point, start);
+        }
+
+        var t = saturate(dot(point - start, segment) / segmentLengthSq);
+        closest = start + segment * t;
+        return distance(point, closest);
+    }
+
+    public static float catmullrom(float p0, float p1, float p2, float p3, float t)
+    {
+        var t2 = t * t;
+        var t3 = t2 * t;
+        return 0.5f * ((2.0f * p1) + (p2 - p0) * t + (2.0f * p0 - 5.0f * p1 + 4.0f * p2 - p3) * t2 + (-p0 + 3.0f * p1 - 3.0f * p2 + p3) * t3);
+    }
+
+    public static float2 catmullrom(float2 p0, float2 p1, float2 p2, float2 p3, float t) =>
+        0.5f * ((2.0f * p1) + (p2 - p0) * t + (2.0f * p0 - 5.0f * p1 + 4.0f * p2 - p3) * (t * t) + (-p0 + 3.0f * p1 - 3.0f * p2 + p3) * (t * t * t));
+
+    public static float3 catmullrom(float3 p0, float3 p1, float3 p2, float3 p3, float t) =>
+        0.5f * ((2.0f * p1) + (p2 - p0) * t + (2.0f * p0 - 5.0f * p1 + 4.0f * p2 - p3) * (t * t) + (-p0 + 3.0f * p1 - 3.0f * p2 + p3) * (t * t * t));
+
+    public static float4 catmullrom(float4 p0, float4 p1, float4 p2, float4 p3, float t) =>
+        0.5f * ((2.0f * p1) + (p2 - p0) * t + (2.0f * p0 - 5.0f * p1 + 4.0f * p2 - p3) * (t * t) + (-p0 + 3.0f * p1 - 3.0f * p2 + p3) * (t * t * t));
+
+    public static float quadratic_bezier(float p0, float p1, float p2, float t)
+    {
+        t = saturate(t);
+        var inv = 1.0f - t;
+        return inv * inv * p0 + 2.0f * inv * t * p1 + t * t * p2;
+    }
+
+    public static float3 quadratic_bezier(float3 p0, float3 p1, float3 p2, float t)
+    {
+        t = saturate(t);
+        var inv = 1.0f - t;
+        return inv * inv * p0 + 2.0f * inv * t * p1 + t * t * p2;
+    }
+
+    public static float cubic_bezier(float p0, float p1, float p2, float p3, float t)
+    {
+        t = saturate(t);
+        var inv = 1.0f - t;
+        return inv * inv * inv * p0 + 3.0f * inv * inv * t * p1 + 3.0f * inv * t * t * p2 + t * t * t * p3;
+    }
+
+    public static float3 cubic_bezier(float3 p0, float3 p1, float3 p2, float3 p3, float t)
+    {
+        t = saturate(t);
+        var inv = 1.0f - t;
+        return inv * inv * inv * p0 + 3.0f * inv * inv * t * p1 + 3.0f * inv * t * t * p2 + t * t * t * p3;
+    }
+
+    public static float smoothstep01(float value) => smoothstep(0.0f, 1.0f, value);
+
+    public static float smootherstep(float value)
+    {
+        value = saturate(value);
+        return value * value * value * (value * (value * 6.0f - 15.0f) + 10.0f);
+    }
+
+    public static float hash(float value) => frac(sin(value) * 43758.5453f);
+    public static float hash(float2 value) => hash(dot(value, new float2(127.1f, 311.7f)));
+    public static float hash(float3 value) => hash(dot(value, new float3(127.1f, 311.7f, 74.7f)));
+
+    public static float value_noise(float2 position)
+    {
+        var cell = floor(position);
+        var local = frac(position);
+        var u = local * local * (3.0f - 2.0f * local);
+
+        var a = hash(cell);
+        var b = hash(cell + new float2(1.0f, 0.0f));
+        var c = hash(cell + new float2(0.0f, 1.0f));
+        var d = hash(cell + new float2(1.0f, 1.0f));
+
+        return lerp(lerp(a, b, u.x), lerp(c, d, u.x), u.y);
+    }
+
+    public static float value_noise_bicubic(float2 position)
+    {
+        var cell = floor(position);
+        var local = frac(position);
+
+        var y0 = catmullrom(
+            hash(cell + new float2(-1.0f, -1.0f)),
+            hash(cell + new float2(0.0f, -1.0f)),
+            hash(cell + new float2(1.0f, -1.0f)),
+            hash(cell + new float2(2.0f, -1.0f)),
+            local.x);
+        var y1 = catmullrom(
+            hash(cell + new float2(-1.0f, 0.0f)),
+            hash(cell + new float2(0.0f, 0.0f)),
+            hash(cell + new float2(1.0f, 0.0f)),
+            hash(cell + new float2(2.0f, 0.0f)),
+            local.x);
+        var y2 = catmullrom(
+            hash(cell + new float2(-1.0f, 1.0f)),
+            hash(cell + new float2(0.0f, 1.0f)),
+            hash(cell + new float2(1.0f, 1.0f)),
+            hash(cell + new float2(2.0f, 1.0f)),
+            local.x);
+        var y3 = catmullrom(
+            hash(cell + new float2(-1.0f, 2.0f)),
+            hash(cell + new float2(0.0f, 2.0f)),
+            hash(cell + new float2(1.0f, 2.0f)),
+            hash(cell + new float2(2.0f, 2.0f)),
+            local.x);
+
+        return catmullrom(y0, y1, y2, y3, local.y);
+    }
 }

tests/CultMath.Tests/HlslMirrorTests.cs

@@ -0,0 +1,63 @@
+using Xunit;
+
+namespace CultMath.Tests;
+
+public sealed class HlslMirrorTests
+{
+    [Fact]
+    public void HlslMirrorPublishesCultMathPrimitives()
+    {
+        var include = File.ReadAllText(Path.Combine(GetRepositoryRoot(), "shaders", "CultMath.hlsl"));
+        var requiredSymbols = new[]
+        {
+            "cultmath_radians",
+            "cultmath_degrees",
+            "cultmath_frac",
+            "cultmath_clamp",
+            "cultmath_saturate",
+            "cultmath_lerp",
+            "cultmath_step",
+            "cultmath_smoothstep",
+            "cultmath_smootherstep",
+            "cultmath_lengthsq",
+            "cultmath_distance",
+            "cultmath_normalize",
+            "cultmath_reflect",
+            "cultmath_csum",
+            "cultmath_decay",
+            "cultmath_damp",
+            "cultmath_catmullrom",
+            "cultmath_quadratic_bezier",
+            "cultmath_cubic_bezier",
+            "cultmath_hash",
+            "cultmath_value_noise",
+            "cultmath_value_noise_bicubic",
+            "cultmath_value_noise_texture",
+            "cultmath_value_noise_texture_bicubic",
+        };
+
+        foreach (var symbol in requiredSymbols)
+        {
+            Assert.Contains(symbol, include);
+        }
+
+        Assert.Contains("CULTMATH_NORMALIZE_EPSILON", include);
+        Assert.Contains("max(length(value), CULTMATH_NORMALIZE_EPSILON)", include);
+    }
+
+    private static string GetRepositoryRoot()
+    {
+        var directory = new DirectoryInfo(AppContext.BaseDirectory);
+        while (directory is not null)
+        {
+            if (File.Exists(Path.Combine(directory.FullName, "CultMath.slnx")))
+            {
+                return directory.FullName;
+            }
+
+            directory = directory.Parent;
+        }
+
+        throw new InvalidOperationException("Could not locate CultMath repository root.");
+    }
+}

tests/CultMath.Tests/MathTests.cs

@@ -23,6 +23,9 @@ public sealed class MathTests
     [Fact]
     public void VectorConstructorsAcceptLowerDimensionVectors()
     {
+        Assert.Equal(new float2(1.0f, 2.0f), math.float2(1.0f, 2.0f));
+        Assert.Equal(new float3(1.0f, 1.0f, 1.0f), math.float3(1.0f));
+        Assert.Equal(new float4(1.0f, 2.0f, 3.0f, 4.0f), math.float4(new float3(1.0f, 2.0f, 3.0f), 4.0f));
         Assert.Equal(new float3(1.0f, 2.0f, 3.0f), new float3(new float2(1.0f, 2.0f), 3.0f));
         Assert.Equal(new float4(1.0f, 2.0f, 3.0f, 4.0f), new float4(new float2(1.0f, 2.0f), 3.0f, 4.0f));
         Assert.Equal(new float4(1.0f, 2.0f, 3.0f, 4.0f), new float4(new float3(1.0f, 2.0f, 3.0f), 4.0f));
@@ -42,6 +45,7 @@ public sealed class MathTests
         Assert.Equal(0.0f, math.smoothstep(0.0f, 1.0f, -1.0f));
         Assert.Equal(0.5f, math.smoothstep(0.0f, 1.0f, 0.5f), precision: 5);
         Assert.Equal(1.0f, math.smoothstep(0.0f, 1.0f, 2.0f));
+        Assert.Equal(0.5f, math.smootherstep(0.5f), precision: 5);
     }
 
     [Fact]
@@ -53,6 +57,50 @@ public sealed class MathTests
         Assert.Equal(math.HALF_PI, math.atan2(1.0f, 0.0f), precision: 5);
     }
 
+    [Fact]
+    public void AetheriaPrimitivesRemainSmallAndDeterministic()
+    {
+        Assert.Equal(6.0f, math.csum(new float3(1.0f, 2.0f, 3.0f)));
+        Assert.Equal(MathF.Exp(-2.0f), math.decay(1.0f, 2.0f, 1.0f), precision: 5);
+        Assert.Equal(1.0f - MathF.Exp(-2.0f), math.damp(0.0f, 1.0f, 2.0f, 1.0f), precision: 5);
+
+        Assert.Equal(15.0f, math.catmullrom(0.0f, 10.0f, 20.0f, 30.0f, 0.5f), precision: 5);
+        Assert.Equal(0.25f, math.quadratic_bezier(0.0f, 0.0f, 1.0f, 0.5f), precision: 5);
+        Assert.Equal(0.5f, math.cubic_bezier(0.0f, 0.0f, 1.0f, 1.0f, 0.5f), precision: 5);
+    }
+
+    [Fact]
+    public void SegmentDistanceReportsClosestPoint()
+    {
+        var distance = math.distance_to_segment(new float2(2.0f, 3.0f), new float2(0.0f, 0.0f), new float2(4.0f, 0.0f), out var closest);
+
+        Assert.Equal(3.0f, distance, precision: 5);
+        Assert.Equal(new float2(2.0f, 0.0f), closest);
+    }
+
+    [Fact]
+    public void FirstOrderInterceptSolvesSimpleClosingTarget()
+    {
+        var time = math.first_order_intercept_time(
+            2.0f,
+            new float3(10.0f, 0.0f, 0.0f),
+            new float3(-1.0f, 0.0f, 0.0f));
+
+        Assert.Equal(10.0f / 3.0f, time, precision: 5);
+    }
+
+    [Fact]
+    public void ValueNoiseIsDeterministicAndSmoothAtCellCorners()
+    {
+        var position = new float2(12.25f, -4.5f);
+
+        Assert.Equal(math.value_noise(position), math.value_noise(position));
+        Assert.Equal(math.hash(new float2(3.0f, 5.0f)), math.value_noise(new float2(3.0f, 5.0f)), precision: 5);
+
+        var bicubic = math.value_noise_bicubic(position);
+        Assert.True(bicubic > -0.25f && bicubic < 1.25f);
+    }
+
     [Fact]
     public void VoronoiBatchProducesToneMappedColors()
     {