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

using System;

// Cut 2a-i fix batch 1: cellular_unit's [0,1) mapping bug (F7) is only reachable from cellular's
// private call sites; this lets CultMath.Tests probe it directly instead of hunting for a hash
// input that happens to land on the bit pattern under test.
[assembly: System.Runtime.CompilerServices.InternalsVisibleTo("CultMath.Tests")]

namespace CultMath;

public static partial class math
{
    public const float PI = MathF.PI;
    public const float TAU = MathF.PI * 2.0f;
    public const float HALF_PI = MathF.PI * 0.5f;

    /// <summary>
    /// Extrinsic Euler order for <see cref="CultMath.float3x3.Euler(CultMath.float3, RotationOrder)"/>:
    /// the first letter is the axis applied first. Nested here so <c>using static CultMath.math;</c>
    /// brings it into scope with the rest of the shader vocabulary.
    /// </summary>
    public enum RotationOrder : byte
    {
        XYZ,
        XZY,
        YXZ,
        YZX,
        ZXY,
        ZYX,
        Default = ZXY,
    }

    public static float2 float2(float x, float y) => new(x, y);
    public static float2 float2(float value) => new(value, value);
    public static bool2 bool2(bool x, bool y) => new(x, y);
    public static bool2 bool2(bool value) => new(value, value);
    public static int2 int2(int x, int y) => new(x, y);
    public static int2 int2(int value) => new(value, value);
    // Mixed constructor functions (float4(float2, float2), ...) are generated in Swizzles.g.cs.
    public static float3 float3(float x, float y, float z) => new(x, y, 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(float value) => new(value, value, value, value);
    public static bool3 bool3(bool x, bool y, bool z) => new(x, y, z);
    public static bool3 bool3(bool value) => new(value, value, value);
    public static bool4 bool4(bool x, bool y, bool z, bool w) => new(x, y, z, w);
    public static bool4 bool4(bool value) => new(value, value, value, value);
    public static int3 int3(int x, int y, int z) => new(x, y, z);
    public static int3 int3(int value) => new(value, value, value);
    public static int4 int4(int x, int y, int z, int w) => new(x, y, z, w);
    public static int4 int4(int value) => new(value, value, value, value);

    // Matrix constructor functions take rows, as HLSL does. Under `using static CultMath.math;`
    // these method groups hide the type names in member access, so static members need a
    // qualified type: `CultMath.float3x3.Euler(...)`, `CultMath.float2x2.Rotate(...)`.
    public static float2x2 float2x2(float m00, float m01, float m10, float m11) => new(m00, m01, m10, m11);
    public static float2x2 float2x2(float2 row0, float2 row1) => new(row0, row1);
    public static float3x3 float3x3(float m00, float m01, float m02, float m10, float m11, float m12, float m20, float m21, float m22) =>
        new(m00, m01, m02, m10, m11, m12, m20, m21, m22);
    public static float3x3 float3x3(float3 row0, float3 row1, float3 row2) => new(row0, row1, row2);
    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);
    public static float cos(float value) => MathF.Cos(value);
    public static float tan(float value) => MathF.Tan(value);
    public static float asin(float value) => MathF.Asin(value);
    public static float acos(float value) => MathF.Acos(value);
    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 pow(float value, float power) => MathF.Pow(value, power);
    public static bool isinf(float value) => float.IsInfinity(value);
    public static float exp(float value) => MathF.Exp(value);
    public static double exp(double value) => Math.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 double2 exp(double2 value) => new(exp(value.x), exp(value.y));
    public static double3 exp(double3 value) => new(exp(value.x), exp(value.y), exp(value.z));

    public static float abs(float value) => MathF.Abs(value);
    public static float2 abs(float2 value) => new(abs(value.x), abs(value.y));
    public static float3 abs(float3 value) => new(abs(value.x), abs(value.y), abs(value.z));
    public static float4 abs(float4 value) => new(abs(value.x), abs(value.y), abs(value.z), abs(value.w));
    public static int abs(int value) => Math.Abs(value);
    public static int2 abs(int2 value) => new(abs(value.x), abs(value.y));
    public static int3 abs(int3 value) => new(abs(value.x), abs(value.y), abs(value.z));
    public static int4 abs(int4 value) => new(abs(value.x), abs(value.y), abs(value.z), abs(value.w));

    // HLSL sign returns int. NaN compares false both ways, so it yields 0 instead of throwing.
    public static int sign(float value) => value > 0.0f ? 1 : value < 0.0f ? -1 : 0;
    public static int2 sign(float2 value) => new(sign(value.x), sign(value.y));
    public static int3 sign(float3 value) => new(sign(value.x), sign(value.y), sign(value.z));
    public static int4 sign(float4 value) => new(sign(value.x), sign(value.y), sign(value.z), sign(value.w));
    public static int sign(int value) => value > 0 ? 1 : value < 0 ? -1 : 0;
    public static int2 sign(int2 value) => new(sign(value.x), sign(value.y));
    public static int3 sign(int3 value) => new(sign(value.x), sign(value.y), sign(value.z));
    public static int4 sign(int4 value) => new(sign(value.x), sign(value.y), sign(value.z), sign(value.w));

    public static float floor(float value) => MathF.Floor(value);
    public static float2 floor(float2 value) => new(floor(value.x), floor(value.y));
    public static float3 floor(float3 value) => new(floor(value.x), floor(value.y), floor(value.z));
    public static float4 floor(float4 value) => new(floor(value.x), floor(value.y), floor(value.z), floor(value.w));

    public static float ceil(float value) => MathF.Ceiling(value);
    public static float2 ceil(float2 value) => new(ceil(value.x), ceil(value.y));
    public static float3 ceil(float3 value) => new(ceil(value.x), ceil(value.y), ceil(value.z));
    public static float4 ceil(float4 value) => new(ceil(value.x), ceil(value.y), ceil(value.z), ceil(value.w));

    public static float frac(float value) => value - floor(value);
    public static float2 frac(float2 value) => value - floor(value);
    public static float3 frac(float3 value) => value - floor(value);
    public static float4 frac(float4 value) => value - floor(value);
    // dxc has no double Frc: it demotes to float, as it does for exp and lerp. CultMath's double
    // overloads keep double precision instead, so frac(double) is x - floor(x) in double.
    public static double frac(double value) => value - Math.Floor(value);

    // dxc lowers float min/max/clamp to DXIL FMin/FMax (saturate to Saturate, defined as FMin(1, FMax(0, x))).
    // DXIL.rst: FMin(a, b) is a < b ? a : b, FMax(a, b) is a >= b ? a : b, and a NaN operand returns the other.
    public static float min(float left, float right) => left < right || float.IsNaN(right) ? left : right;
    public static int min(int left, int right) => Math.Min(left, right);
    public static float2 min(float2 left, float2 right) => new(min(left.x, right.x), min(left.y, right.y));
    public static float3 min(float3 left, float3 right) => new(min(left.x, right.x), min(left.y, right.y), min(left.z, right.z));
    public static float4 min(float4 left, float4 right) => new(min(left.x, right.x), min(left.y, right.y), min(left.z, right.z), min(left.w, right.w));
    public static int2 min(int2 left, int2 right) => new(min(left.x, right.x), min(left.y, right.y));
    public static int3 min(int3 left, int3 right) => new(min(left.x, right.x), min(left.y, right.y), min(left.z, right.z));
    public static int4 min(int4 left, int4 right) => new(min(left.x, right.x), min(left.y, right.y), min(left.z, right.z), min(left.w, right.w));

    public static float max(float left, float right) => left >= right || float.IsNaN(right) ? left : right;
    public static int max(int left, int right) => Math.Max(left, right);
    public static float2 max(float2 left, float2 right) => new(max(left.x, right.x), max(left.y, right.y));
    public static float3 max(float3 left, float3 right) => new(max(left.x, right.x), max(left.y, right.y), max(left.z, right.z));
    public static float4 max(float4 left, float4 right) => new(max(left.x, right.x), max(left.y, right.y), max(left.z, right.z), max(left.w, right.w));
    public static int2 max(int2 left, int2 right) => new(max(left.x, right.x), max(left.y, right.y));
    public static int3 max(int3 left, int3 right) => new(max(left.x, right.x), max(left.y, right.y), max(left.z, right.z));
    public static int4 max(int4 left, int4 right) => new(max(left.x, right.x), max(left.y, right.y), max(left.z, right.z), max(left.w, right.w));

    public static float clamp(float value, float minimum, float maximum) => min(max(value, minimum), maximum);
    public static float2 clamp(float2 value, float2 minimum, float2 maximum) => min(max(value, minimum), maximum);
    public static float3 clamp(float3 value, float3 minimum, float3 maximum) => min(max(value, minimum), maximum);
    public static float4 clamp(float4 value, float4 minimum, float4 maximum) => min(max(value, minimum), maximum);
    // dxc lowers int clamp to IMin(IMax(x, a), b): inverted bounds return b rather than throwing like Math.Clamp.
    public static int clamp(int value, int minimum, int maximum) => min(max(value, minimum), maximum);
    public static int2 clamp(int2 value, int2 minimum, int2 maximum) => min(max(value, minimum), maximum);
    public static int3 clamp(int3 value, int3 minimum, int3 maximum) => min(max(value, minimum), maximum);
    public static int4 clamp(int4 value, int4 minimum, int4 maximum) => min(max(value, minimum), maximum);

    // DXIL Saturate is FMin(1, FMax(0, x)); that operand order, unlike clamp(x, 0, 1), maps -0 to +0.
    public static float saturate(float value) => min(1.0f, max(0.0f, value));
    public static float2 saturate(float2 value) => min(1.0f, max(0.0f, value));
    public static float3 saturate(float3 value) => min(1.0f, max(0.0f, value));
    public static float4 saturate(float4 value) => min(1.0f, max(0.0f, value));

    public static float lerp(float start, float end, float amount) => start + (end - start) * amount;
    public static double lerp(double start, double end, double amount) => start + (end - start) * amount;
    public static float2 lerp(float2 start, float2 end, float2 amount) => start + (end - start) * amount;
    public static float3 lerp(float3 start, float3 end, float3 amount) => start + (end - start) * amount;
    public static float4 lerp(float4 start, float4 end, float4 amount) => start + (end - start) * amount;
    public static double2 lerp(double2 start, double2 end, double2 amount) => start + (end - start) * amount;
    public static double3 lerp(double3 start, double3 end, double3 amount) => start + (end - start) * amount;

    public static float unlerp(float start, float end, float value) => (value - start) / (end - start);
    public static double unlerp(double start, double end, double value) => (value - start) / (end - start);
    public static float2 unlerp(float2 start, float2 end, float2 value) => (value - start) / (end - start);
    public static float3 unlerp(float3 start, float3 end, float3 value) => (value - start) / (end - start);
    public static float4 unlerp(float4 start, float4 end, float4 value) => (value - start) / (end - start);
    public static double2 unlerp(double2 start, double2 end, double2 value) => (value - start) / (end - start);
    public static double3 unlerp(double3 start, double3 end, double3 value) => (value - start) / (end - start);

    // dxc lowers step(y, x) to x < y ? 0 : 1 (fcmp olt, then select), so a NaN on either side yields 1.
    public static float step(float edge, float value) => value < edge ? 0.0f : 1.0f;
    public static float2 step(float2 edge, float2 value) => new(step(edge.x, value.x), step(edge.y, value.y));
    public static float3 step(float3 edge, float3 value) => new(step(edge.x, value.x), step(edge.y, value.y), step(edge.z, value.z));
    public static float4 step(float4 edge, float4 value) => new(step(edge.x, value.x), step(edge.y, value.y), step(edge.z, value.z), step(edge.w, value.w));

    public static float smoothstep(float minimum, float maximum, float value)
    {
        var t = saturate((value - minimum) / (maximum - minimum));
        return t * t * (3.0f - 2.0f * t);
    }

    public static float2 smoothstep(float2 minimum, float2 maximum, float2 value)
    {
        var t = saturate((value - minimum) / (maximum - minimum));
        return t * t * (3.0f - 2.0f * t);
    }

    public static float3 smoothstep(float3 minimum, float3 maximum, float3 value)
    {
        var t = saturate((value - minimum) / (maximum - minimum));
        return t * t * (3.0f - 2.0f * t);
    }

    public static float4 smoothstep(float4 minimum, float4 maximum, float4 value)
    {
        var t = saturate((value - minimum) / (maximum - minimum));
        return t * t * (3.0f - 2.0f * t);
    }

    public static float dot(float2 left, float2 right) => left.x * right.x + left.y * right.y;
    public static float dot(float3 left, float3 right) => left.x * right.x + left.y * right.y + left.z * right.z;
    public static float dot(float4 left, float4 right) => left.x * right.x + left.y * right.y + left.z * right.z + left.w * right.w;

    public static float3 cross(float3 left, float3 right) =>
        new(
            left.y * right.z - left.z * right.y,
            left.z * right.x - left.x * right.z,
            left.x * right.y - left.y * right.x);

    public static float lengthsq(float2 value) => dot(value, value);
    public static float lengthsq(float3 value) => dot(value, value);
    public static float lengthsq(float4 value) => dot(value, value);

    public static float length(float2 value) => MathF.Sqrt(lengthsq(value));
    public static float length(float3 value) => MathF.Sqrt(lengthsq(value));
    public static float length(float4 value) => MathF.Sqrt(lengthsq(value));

    public static float distance(float2 left, float2 right) => length(left - right);
    public static float distance(float3 left, float3 right) => length(left - right);
    public static float distance(float4 left, float4 right) => length(left - right);

    // dxc lowers normalize(x) to x * Rsqrt(Dot(x, x)), and DXIL Rsqrt is 1 / sqrt(src): a zero vector is 0 * inf = NaN.
    public static float2 normalize(float2 value) => value * (1.0f / MathF.Sqrt(dot(value, value)));
    public static float3 normalize(float3 value) => value * (1.0f / MathF.Sqrt(dot(value, value)));
    public static float4 normalize(float4 value) => value * (1.0f / MathF.Sqrt(dot(value, value)));
    public static quaternion normalize(quaternion value)
    {
        var length = MathF.Sqrt(value.x * value.x + value.y * value.y + value.z * value.z + value.w * value.w);
        if (length <= 1.0e-20f)
        {
            return quaternion.identity;
        }

        var invLength = 1.0f / length;
        return new quaternion(value.x * invLength, value.y * invLength, value.z * invLength, value.w * invLength);
    }

    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 float2 rotate(float2 value, float radians)
    {
        var s = sin(radians);
        var c = cos(radians);
        return new float2(value.x * c - value.y * s, value.x * s + value.y * c);
    }

    public static float2 rotate_degrees(float2 value, float degrees) => rotate(value, radians(degrees));

    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)
    {
        closest = segment_closest(point, start, end);
        return distance(point, closest);
    }

    /// <summary>The distance from <paramref name="point"/> to the segment <paramref name="start"/>-<paramref name="end"/>; a degenerate segment is its start.</summary>
    public static float distance_to_segment(float2 point, float2 start, float2 end) =>
        distance(point, segment_closest(point, start, end));

    /// <summary>The distance from <paramref name="point"/> to the segment <paramref name="start"/>-<paramref name="end"/>; a degenerate segment is its start.</summary>
    public static float distance_to_segment(float3 point, float3 start, float3 end) =>
        distance(point, segment_closest(point, start, end));

    /// <summary>The point of the segment <paramref name="start"/>-<paramref name="end"/> nearest <paramref name="point_value"/>; a degenerate segment is its start.</summary>
    public static float2 segment_closest(float2 point_value, float2 start, float2 end)
    {
        var segment = end - start;
        var segmentLengthSq = lengthsq(segment);
        if (segmentLengthSq <= 0.0f)
            return start;

        var t = saturate(dot(point_value - start, segment) / segmentLengthSq);
        return start + segment * t;
    }

    /// <summary>The point of the segment <paramref name="start"/>-<paramref name="end"/> nearest <paramref name="point_value"/>; a degenerate segment is its start.</summary>
    public static float3 segment_closest(float3 point_value, float3 start, float3 end)
    {
        var segment = end - start;
        var segmentLengthSq = lengthsq(segment);
        if (segmentLengthSq <= 0.0f)
            return start;

        var t = saturate(dot(point_value - start, segment) / segmentLengthSq);
        return start + segment * t;
    }

    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 uint asuint(float value) => (uint)BitConverter.SingleToInt32Bits(value);
    public static float asfloat(uint value) => BitConverter.Int32BitsToSingle((int)value);

    // PCG hashes from Jarzynski and Olano, "Hash Functions for GPU Rendering" (JCGT 9(3), 2020). Integer-only, so
    // dxc and C# agree bit for bit. pcg3d/pcg4d carry uint bit patterns in int3/int4; float overloads hash IEEE bits.
    public static uint pcg(uint value)
    {
        var state = value * 747796405u + 2891336453u;
        var word = ((state >> (int)((state >> 28) + 4)) ^ state) * 277803737u;
        return (word >> 22) ^ word;
    }

    public static int3 pcg3d(int3 value)
    {
        uint x = (uint)value.x * 1664525u + 1013904223u, y = (uint)value.y * 1664525u + 1013904223u, z = (uint)value.z * 1664525u + 1013904223u;
        x += y * z; y += z * x; z += x * y;
        x ^= x >> 16; y ^= y >> 16; z ^= z >> 16;
        x += y * z; y += z * x; z += x * y;
        return new int3((int)x, (int)y, (int)z);
    }

    // The paper hashes 2 -> N inputs through pcg3d with the unused component held constant.
    public static int3 pcg3d(float2 value) => pcg3d(new float3(value, 0.0f));
    public static int3 pcg3d(float3 value) => pcg3d(new int3((int)asuint(value.x), (int)asuint(value.y), (int)asuint(value.z)));

    public static int4 pcg4d(int4 value)
    {
        uint x = (uint)value.x * 1664525u + 1013904223u, y = (uint)value.y * 1664525u + 1013904223u;
        uint z = (uint)value.z * 1664525u + 1013904223u, w = (uint)value.w * 1664525u + 1013904223u;
        x += y * w; y += z * x; z += x * y; w += y * z;
        x ^= x >> 16; y ^= y >> 16; z ^= z >> 16; w ^= w >> 16;
        x += y * w; y += z * x; z += x * y; w += y * z;
        return new int4((int)x, (int)y, (int)z, (int)w);
    }

    public static int4 pcg4d(float4 value) =>
        pcg4d(new int4((int)asuint(value.x), (int)asuint(value.y), (int)asuint(value.z), (int)asuint(value.w)));

    public static float log(float value) => MathF.Log(value);
    public static float2 log(float2 value) => new(log(value.x), log(value.y));
    public static float3 log(float3 value) => new(log(value.x), log(value.y), log(value.z));
    public static float4 log(float4 value) => new(log(value.x), log(value.y), log(value.z), log(value.w));

    public static float2 pow(float2 value, float2 power) => new(pow(value.x, power.x), pow(value.y, power.y));
    public static float3 pow(float3 value, float3 power) => new(pow(value.x, power.x), pow(value.y, power.y), pow(value.z, power.z));
    public static float4 pow(float4 value, float4 power) => new(pow(value.x, power.x), pow(value.y, power.y), pow(value.z, power.z), pow(value.w, power.w));

    public static float2 sqrt(float2 value) => new(sqrt(value.x), sqrt(value.y));
    public static float3 sqrt(float3 value) => new(sqrt(value.x), sqrt(value.y), sqrt(value.z));
    public static float4 sqrt(float4 value) => new(sqrt(value.x), sqrt(value.y), sqrt(value.z), sqrt(value.w));

    public static float2 sin(float2 value) => new(sin(value.x), sin(value.y));
    public static float3 sin(float3 value) => new(sin(value.x), sin(value.y), sin(value.z));
    public static float4 sin(float4 value) => new(sin(value.x), sin(value.y), sin(value.z), sin(value.w));

    public static float2 cos(float2 value) => new(cos(value.x), cos(value.y));
    public static float3 cos(float3 value) => new(cos(value.x), cos(value.y), cos(value.z));
    public static float4 cos(float4 value) => new(cos(value.x), cos(value.y), cos(value.z), cos(value.w));

    public static float2 acos(float2 value) => new(acos(value.x), acos(value.y));
    public static float3 acos(float3 value) => new(acos(value.x), acos(value.y), acos(value.z));
    public static float4 acos(float4 value) => new(acos(value.x), acos(value.y), acos(value.z), acos(value.w));

    public static float2 atan2(float2 y, float2 x) => new(atan2(y.x, x.x), atan2(y.y, x.y));
    public static float3 atan2(float3 y, float3 x) => new(atan2(y.x, x.x), atan2(y.y, x.y), atan2(y.z, x.z));
    public static float4 atan2(float4 y, float4 x) => new(atan2(y.x, x.x), atan2(y.y, x.y), atan2(y.z, x.z), atan2(y.w, x.w));

    // DXIL FMin/FMax take double overloads with the same rules as float (see min/max above).
    public static double min(double left, double right) => left < right || double.IsNaN(right) ? left : right;
    public static double max(double left, double right) => left >= right || double.IsNaN(right) ? left : right;

    public static bool any(bool2 value) => value.x || value.y;
    public static bool any(bool3 value) => value.x || value.y || value.z;
    public static bool any(bool4 value) => value.x || value.y || value.z || value.w;
    public static bool all(bool2 value) => value.x && value.y;
    public static bool all(bool3 value) => value.x && value.y && value.z;
    public static bool all(bool4 value) => value.x && value.y && value.z && value.w;

    // HLSL any/all on numeric vectors test components against zero (NaN != 0 counts as true).
    public static bool any(float2 value) => value.x != 0.0f || value.y != 0.0f;
    public static bool any(float3 value) => value.x != 0.0f || value.y != 0.0f || value.z != 0.0f;
    public static bool any(float4 value) => value.x != 0.0f || value.y != 0.0f || value.z != 0.0f || value.w != 0.0f;
    public static bool all(float2 value) => value.x != 0.0f && value.y != 0.0f;
    public static bool all(float3 value) => value.x != 0.0f && value.y != 0.0f && value.z != 0.0f;
    public static bool all(float4 value) => value.x != 0.0f && value.y != 0.0f && value.z != 0.0f && value.w != 0.0f;
    public static bool any(int2 value) => value.x != 0 || value.y != 0;
    public static bool any(int3 value) => value.x != 0 || value.y != 0 || value.z != 0;
    public static bool any(int4 value) => value.x != 0 || value.y != 0 || value.z != 0 || value.w != 0;
    public static bool all(int2 value) => value.x != 0 && value.y != 0;
    public static bool all(int3 value) => value.x != 0 && value.y != 0 && value.z != 0;
    public static bool all(int4 value) => value.x != 0 && value.y != 0 && value.z != 0 && value.w != 0;

    // HLSL 2021 argument order: select(condition, whenTrue, whenFalse).
    // Unity.Mathematics' select(falseValue, trueValue, condition) is the reverse.
    public static float select(bool condition, float whenTrue, float whenFalse) => condition ? whenTrue : whenFalse;
    public static float2 select(bool2 condition, float2 whenTrue, float2 whenFalse) =>
        new(condition.x ? whenTrue.x : whenFalse.x, condition.y ? whenTrue.y : whenFalse.y);
    public static float3 select(bool3 condition, float3 whenTrue, float3 whenFalse) =>
        new(condition.x ? whenTrue.x : whenFalse.x, condition.y ? whenTrue.y : whenFalse.y, condition.z ? whenTrue.z : whenFalse.z);
    public static float4 select(bool4 condition, float4 whenTrue, float4 whenFalse) =>
        new(condition.x ? whenTrue.x : whenFalse.x, condition.y ? whenTrue.y : whenFalse.y, condition.z ? whenTrue.z : whenFalse.z, condition.w ? whenTrue.w : whenFalse.w);

    // HLSL mul: mul(m, v) treats v as a column vector, mul(v, m) as a row vector.
    public static float mul(float left, float right) => left * right;

    public static float2 mul(float2x2 m, float2 v) =>
        new(m._m00 * v.x + m._m01 * v.y, m._m10 * v.x + m._m11 * v.y);

    public static float2 mul(float2 v, float2x2 m) =>
        new(v.x * m._m00 + v.y * m._m10, v.x * m._m01 + v.y * m._m11);

    public static float2x2 mul(float2x2 a, float2x2 b) =>
        new(
            a._m00 * b._m00 + a._m01 * b._m10, a._m00 * b._m01 + a._m01 * b._m11,
            a._m10 * b._m00 + a._m11 * b._m10, a._m10 * b._m01 + a._m11 * b._m11);

    public static float3 mul(float3x3 m, float3 v) =>
        new(
            m._m00 * v.x + m._m01 * v.y + m._m02 * v.z,
            m._m10 * v.x + m._m11 * v.y + m._m12 * v.z,
            m._m20 * v.x + m._m21 * v.y + m._m22 * v.z);

    public static float3 mul(float3 v, float3x3 m) =>
        new(
            v.x * m._m00 + v.y * m._m10 + v.z * m._m20,
            v.x * m._m01 + v.y * m._m11 + v.z * m._m21,
            v.x * m._m02 + v.y * m._m12 + v.z * m._m22);

    public static float3x3 mul(float3x3 a, float3x3 b) =>
        new(
            a._m00 * b._m00 + a._m01 * b._m10 + a._m02 * b._m20,
            a._m00 * b._m01 + a._m01 * b._m11 + a._m02 * b._m21,
            a._m00 * b._m02 + a._m01 * b._m12 + a._m02 * b._m22,
            a._m10 * b._m00 + a._m11 * b._m10 + a._m12 * b._m20,
            a._m10 * b._m01 + a._m11 * b._m11 + a._m12 * b._m21,
            a._m10 * b._m02 + a._m11 * b._m12 + a._m12 * b._m22,
            a._m20 * b._m00 + a._m21 * b._m10 + a._m22 * b._m20,
            a._m20 * b._m01 + a._m21 * b._m11 + a._m22 * b._m21,
            a._m20 * b._m02 + a._m21 * b._m12 + a._m22 * b._m22);

    // Abramowitz & Stegun 7.1.26, a single-precision rational/exponential fit with a stated
    // max absolute error of 1.5e-7. erf is odd, so only x >= 0 is fit and the sign is restored.
    // Measured against an independent reference: worst absolute error 6.621e-7, inside the
    // 2e-6 bar this package's tests pin.
    public static float erf(float value)
    {
        const float p = 0.3275911f;
        const float a1 = 0.254829592f;
        const float a2 = -0.284496736f;
        const float a3 = 1.421413741f;
        const float a4 = -1.453152027f;
        const float a5 = 1.061405429f;

        var s = sign(value);
        var x = abs(value);
        var t = 1.0f / (1.0f + p * x);
        var poly = ((((a5 * t + a4) * t + a3) * t + a2) * t + a1) * t;
        var result = 1.0f - poly * exp(-x * x);
        return s * result;
    }

    // Giles, "Approximating the erfinv function" (GPU Computing Gems, 2010): a single-precision
    // minimax polynomial in two ranges of w = -ln((1-x)(1+x)), central and tail.
    // Domain is [-1, 1]. Outside it, including +-infinity, the result is NaN (there is no real
    // inverse). At the edges the mathematical limit is infinite: erfinv(1) = +infinity and
    // erfinv(-1) = -infinity. The unguarded polynomial hits log(0) = -infinity at those edges and
    // the tail branch's leading coefficient is negative, so it delivers the wrong sign; guard the
    // edges explicitly instead of trusting the fit through the singularity.
    // Measured against an independent reference: worst absolute error 5.066e-7, and the
    // erf(erfinv(y)) round trip 6.109e-7 worst case, both inside the 2e-6 bar this package's
    // tests pin.
    public static float erfinv(float value)
    {
        if (value >= 1.0f) return value > 1.0f ? float.NaN : float.PositiveInfinity;
        if (value <= -1.0f) return value < -1.0f ? float.NaN : float.NegativeInfinity;

        var w = -log((1.0f - value) * (1.0f + value));
        float p;
        if (w < 5.0f)
        {
            w -= 2.5f;
            p = 2.81022636e-08f;
            p = 3.43273939e-07f + p * w;
            p = -3.5233877e-06f + p * w;
            p = -4.39150654e-06f + p * w;
            p = 0.00021858087f + p * w;
            p = -0.00125372503f + p * w;
            p = -0.00417768164f + p * w;
            p = 0.246640727f + p * w;
            p = 1.50140941f + p * w;
        }
        else
        {
            w = sqrt(w) - 3.0f;
            p = -0.000200214257f;
            p = 0.000100950558f + p * w;
            p = 0.00134934322f + p * w;
            p = -0.00367342844f + p * w;
            p = 0.00573950773f + p * w;
            p = -0.0076224613f + p * w;
            p = 0.00943887047f + p * w;
            p = 1.00167406f + p * w;
            p = 2.83297682f + p * w;
        }

        return p * value;
    }

    // Ashima Arts / Ian McEwan 3D simplex noise (MIT), mirrored by cultmath_snoise(float3)
    // in shaders/CultMath.hlsl with the same float32 evaluation order. Follows stegu/webgl-noise
    // src/noise3D.glsl at 22434e04d7 (kernel radius^2 0.5, scale 105, permute (34x+10)x); see
    // THIRD-PARTY-NOTICES.md. Continuous across simplex-cell boundaries
    // (NoiseGradTests.SnoiseAndItsGradientAreContinuousAcrossSimplexCellBoundaries). Ties in x0 resolve
    // in Gustavson's simplexnoise1234 order (x over y, y over z, x over z), a total order, so the four
    // corners are always a simplex; upstream's step is cyclic and degenerates where x0's components are
    // equal (NoiseGradTests.SnoiseIsContinuousOnTheSimplexDiagonal).
    // SNOISE_LIPSCHITZ (math.Interval.cs) is measured from this kernel: a change here must re-run
    // NoiseBoundTests.MeasureLipschitz and re-pin it.
    public static float snoise(float3 value)
    {
        const float cx = 1.0f / 6.0f;
        const float cy = 1.0f / 3.0f;
        var i = floor(value + dot(value, new float3(cy, cy, cy)));
        var x0 = value - i + dot(i, new float3(cx, cx, cx));
        var g = new float3(step(x0.y, x0.x), step(x0.z, x0.y), 1.0f - step(x0.z, x0.x));
        var l = 1.0f - g;
        var lzxy = new float3(l.z, l.x, l.y);
        var i1 = min(g, lzxy);
        var i2 = max(g, lzxy);
        var x1 = x0 - i1 + cx;
        var x2 = x0 - i2 + cy;
        var x3 = x0 - 0.5f;

        i = snoise_mod289(i);
        var p = snoise_permute(snoise_permute(snoise_permute(
            i.z + new float4(0.0f, i1.z, i2.z, 1.0f)) + i.y + new float4(0.0f, i1.y, i2.y, 1.0f)) + i.x + new float4(0.0f, i1.x, i2.x, 1.0f));

        const float n = 0.142857142857f;
        var ns = new float3(2.0f * n, 0.5f * n - 1.0f, n);
        var j = p - 49.0f * floor(p * ns.z * ns.z);
        var xs = floor(j * ns.z);
        var ys = floor(j - 7.0f * xs);
        var x = xs * ns.x + ns.y;
        var y = ys * ns.x + ns.y;
        var h = 1.0f - abs(x) - abs(y);

        var b0 = new float4(x.x, x.y, y.x, y.y);
        var b1 = new float4(x.z, x.w, y.z, y.w);
        var s0 = floor(b0) * 2.0f + 1.0f;
        var s1 = floor(b1) * 2.0f + 1.0f;
        var sh = -step(h, 0.0f);
        var a0 = new float4(b0.x, b0.z, b0.y, b0.w) + new float4(s0.x, s0.z, s0.y, s0.w) * new float4(sh.x, sh.x, sh.y, sh.y);
        var a1 = new float4(b1.x, b1.z, b1.y, b1.w) + new float4(s1.x, s1.z, s1.y, s1.w) * new float4(sh.z, sh.z, sh.w, sh.w);

        var p0 = new float3(a0.x, a0.y, h.x);
        var p1 = new float3(a0.z, a0.w, h.y);
        var p2 = new float3(a1.x, a1.y, h.z);
        var p3 = new float3(a1.z, a1.w, h.w);
        var norm = 1.79284291400159f - 0.85373472095314f * new float4(dot(p0, p0), dot(p1, p1), dot(p2, p2), dot(p3, p3));
        p0 *= norm.x;
        p1 *= norm.y;
        p2 *= norm.z;
        p3 *= norm.w;

        var m = max(0.5f - new float4(dot(x0, x0), dot(x1, x1), dot(x2, x2), dot(x3, x3)), 0.0f);
        m = m * m;
        return 105.0f * dot(m * m, new float4(dot(p0, x0), dot(p1, x1), dot(p2, x2), dot(p3, x3)));
    }

    // Ashima Arts / Ian McEwan 2D simplex noise, kept component-explicit so
    // the C# and HLSL mirrors preserve the same float32 evaluation order.
    public static float snoise(float2 value)
    {
        var c = new float4(
            0.211324865405187f,
            0.366025403784439f,
            -0.577350269189626f,
            0.024390243902439f);
        var i = floor(value + dot(value, new float2(c.y, c.y)));
        var x0 = value - i + dot(i, new float2(c.x, c.x));
        var i1 = x0.x > x0.y ? new float2(1.0f, 0.0f) : new float2(0.0f, 1.0f);
        var x12 = new float4(
            x0.x + c.x - i1.x,
            x0.y + c.x - i1.y,
            x0.x + c.z,
            x0.y + c.z);

        i = snoise_mod289(i);
        var p = snoise_permute(snoise_permute(i.y + new float3(0.0f, i1.y, 1.0f)) + i.x + new float3(0.0f, i1.x, 1.0f));
        var m = max(0.5f - new float3(
            dot(x0, x0),
            x12.x * x12.x + x12.y * x12.y,
            x12.z * x12.z + x12.w * x12.w), 0.0f);
        m = m * m;
        m = m * m;

        var x = 2.0f * frac(p * c.w) - 1.0f;
        var h = abs(x) - 0.5f;
        var ox = floor(x + 0.5f);
        var a0 = x - ox;
        m *= 1.79284291400159f - 0.85373472095314f * (a0 * a0 + h * h);

        var g = new float3(
            a0.x * x0.x + h.x * x0.y,
            a0.y * x12.x + h.y * x12.y,
            a0.z * x12.z + h.z * x12.w);
        return 130.0f * dot(m, g);
    }

    private static float2 snoise_mod289(float2 value) => value - floor(value * (1.0f / 289.0f)) * 289.0f;
    private static float3 snoise_mod289(float3 value) => value - floor(value * (1.0f / 289.0f)) * 289.0f;
    private static float4 snoise_mod289(float4 value) => value - floor(value * (1.0f / 289.0f)) * 289.0f;
    private static float3 snoise_permute(float3 value) => snoise_mod289(((value * 34.0f) + 10.0f) * value);
    private static float4 snoise_permute(float4 value) => snoise_mod289(((value * 34.0f) + 10.0f) * value);

    // Analytic gradient of snoise(float3) (invariant 8), following the differentiation in
    // webgl-noise src/noise3Dgrad.glsl (Ashima Arts / McEwan, MIT), which uses the same kernel
    // radius^2 0.5, permutation (34x+10)x and scale 105 as snoise(float3) (stegu/webgl-noise
    // 22434e04d7). Every local up through p0..p3 and x0..x3 below is snoise(float3)'s own derivation
    // verbatim, so .w is that same field (NoiseGradTests.SnoiseGradValueIsBitEqualToSnoise pins the
    // two bit-equal over a random sweep). Per corner, m0 = max(0.5 - dot(x,x), 0), so dm0/dx = -2x,
    // and d(m0^4 * dot(p,x))/dx = 4*m0^3*(-2x)*dot(p,x) + m0^4*p = -8*m0^3*dot(p,x)*x + m0^4*p (p is
    // that corner's own gradient constant, independent of x); summed over the four corners and
    // scaled by the same 105 as the value.
    public static float4 snoise_grad(float3 value)
    {
        const float cx = 1.0f / 6.0f;
        const float cy = 1.0f / 3.0f;
        var i = floor(value + dot(value, new float3(cy, cy, cy)));
        var x0 = value - i + dot(i, new float3(cx, cx, cx));
        var g = new float3(step(x0.y, x0.x), step(x0.z, x0.y), 1.0f - step(x0.z, x0.x));
        var l = 1.0f - g;
        var lzxy = new float3(l.z, l.x, l.y);
        var i1 = min(g, lzxy);
        var i2 = max(g, lzxy);
        var x1 = x0 - i1 + cx;
        var x2 = x0 - i2 + cy;
        var x3 = x0 - 0.5f;

        i = snoise_mod289(i);
        var p = snoise_permute(snoise_permute(snoise_permute(
            i.z + new float4(0.0f, i1.z, i2.z, 1.0f)) + i.y + new float4(0.0f, i1.y, i2.y, 1.0f)) + i.x + new float4(0.0f, i1.x, i2.x, 1.0f));

        const float n = 0.142857142857f;
        var ns = new float3(2.0f * n, 0.5f * n - 1.0f, n);
        var j = p - 49.0f * floor(p * ns.z * ns.z);
        var xs = floor(j * ns.z);
        var ys = floor(j - 7.0f * xs);
        var x = xs * ns.x + ns.y;
        var y = ys * ns.x + ns.y;
        var h = 1.0f - abs(x) - abs(y);

        var b0 = new float4(x.x, x.y, y.x, y.y);
        var b1 = new float4(x.z, x.w, y.z, y.w);
        var s0 = floor(b0) * 2.0f + 1.0f;
        var s1 = floor(b1) * 2.0f + 1.0f;
        var sh = -step(h, 0.0f);
        var a0 = new float4(b0.x, b0.z, b0.y, b0.w) + new float4(s0.x, s0.z, s0.y, s0.w) * new float4(sh.x, sh.x, sh.y, sh.y);
        var a1 = new float4(b1.x, b1.z, b1.y, b1.w) + new float4(s1.x, s1.z, s1.y, s1.w) * new float4(sh.z, sh.z, sh.w, sh.w);

        var p0 = new float3(a0.x, a0.y, h.x);
        var p1 = new float3(a0.z, a0.w, h.y);
        var p2 = new float3(a1.x, a1.y, h.z);
        var p3 = new float3(a1.z, a1.w, h.w);
        var norm = 1.79284291400159f - 0.85373472095314f * new float4(dot(p0, p0), dot(p1, p1), dot(p2, p2), dot(p3, p3));
        p0 *= norm.x;
        p1 *= norm.y;
        p2 *= norm.z;
        p3 *= norm.w;

        var m0 = max(0.5f - new float4(dot(x0, x0), dot(x1, x1), dot(x2, x2), dot(x3, x3)), 0.0f);
        var m2 = m0 * m0;
        var m3 = m2 * m0;
        var m4 = m2 * m2;

        var px = new float4(dot(p0, x0), dot(p1, x1), dot(p2, x2), dot(p3, x3));
        var value2 = 105.0f * dot(m4, px);
        var grad = -8.0f * m3.x * x0 * px.x + m4.x * p0
            + -8.0f * m3.y * x1 * px.y + m4.y * p1
            + -8.0f * m3.z * x2 * px.z + m4.z * p2
            + -8.0f * m3.w * x3 * px.w + m4.w * p3;
        grad *= 105.0f;

        return new float4(grad, value2);
    }

    // Octave sum of snoise_grad (invariant 8). Each octave samples snoise_grad at p scaled by the
    // running frequency, so by the chain rule its gradient scales by that same frequency
    // (d/dp[n(f*p)] = f * grad_n(f*p)); amplitude and frequency both compound per octave (gain,
    // lacunarity respectively), not applied once for the whole sum, so with 1 octave this is exactly
    // snoise_grad (NoiseGradTests.FbmGradWithOneOctaveEqualsSnoiseGrad). octaves is clamped to
    // [0, 16]: HlslSourceCompatibilityTests' bit-parity oracle drives every mirrored int parameter
    // across the full int32 range, 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 this mirror.
    public static float4 fbm_grad(float3 p, int octaves, float lacunarity, float gain)
    {
        octaves = clamp(octaves, 0, 16);
        var amplitude = 1.0f;
        var frequency = 1.0f;
        var value = 0.0f;
        var gradient = new float3(0.0f, 0.0f, 0.0f);
        for (var i = 0; i < octaves; i++)
        {
            var n = snoise_grad(p * frequency);
            value += amplitude * n.w;
            gradient += amplitude * frequency * new float3(n.x, n.y, n.z);
            frequency *= lacunarity;
            amplitude *= gain;
        }

        return new float4(gradient, value);
    }

    // Value-only fractional Brownian motion: the octave sum of snoise, octave i at frequency lacunarity^i and
    // amplitude gain^i. It is fbm_grad(p, octaves, lacunarity, gain).w without the gradient work (to float
    // rounding: snoise and snoise_grad.w order their arithmetic differently), so a loop that only wants the height
    // (a sphere trace) does not pay for the gradient. octaves clamps to [0, 16] as fbm_grad's does. The one-argument
    // form is the stock four octaves, lacunarity 2, gain 1/2: its range is within +-1.875 times snoise's.
    public static float fbm(float3 p, int octaves, float lacunarity, float gain)
    {
        octaves = clamp(octaves, 0, 16);
        var amplitude = 1.0f;
        var frequency = 1.0f;
        var value = 0.0f;
        for (var i = 0; i < octaves; i++)
        {
            value += amplitude * snoise(p * frequency);
            frequency *= lacunarity;
            amplitude *= gain;
        }

        return value;
    }

    public static float fbm(float3 value) => fbm(value, 4, 2.0f, 0.5f);

    // Ridged multifractal noise (Musgrave, Texturing and Modeling, ch. 16; the operator's dune term,
    // 2026-09-25): each octave folds snoise_grad about zero, value = Sum a_i*(1 - |n_i|), gradient =
    // -Sum a_i*f_i*sign(n_i)*grad(n_i), the same per-octave frequency scaling as fbm_grad's chain
    // rule plus the fold's sign flip (d|n|/dp = sign(n)*dn/dp). The crease at n_i = 0 is real and
    // deliberate, not a bug: NoiseGradTests excludes a band around it, the same way cellular excludes
    // its F1 = F2 seam. sign follows this codebase's HLSL-matching convention (design.md: sign
    // returns int, 0 at exactly 0), so the crease itself contributes no gradient rather than an
    // arbitrary one. Same octave clamp as fbm_grad, for the same bit-parity-oracle reason.
    public static float4 ridged_grad(float3 p, int octaves, float lacunarity, float gain)
    {
        octaves = clamp(octaves, 0, 16);
        var amplitude = 1.0f;
        var frequency = 1.0f;
        var value = 0.0f;
        var gradient = new float3(0.0f, 0.0f, 0.0f);
        for (var i = 0; i < octaves; i++)
        {
            var n = snoise_grad(p * frequency);
            var s = (float)sign(n.w);
            value += amplitude * (1.0f - abs(n.w));
            gradient += -amplitude * frequency * s * new float3(n.x, n.y, n.z);
            frequency *= lacunarity;
            amplitude *= gain;
        }

        return new float4(gradient, value);
    }

    // Inigo Quilez, "smooth minimum" (https://iquilezles.org/articles/smin/): the quadratic-polynomial
    // smin, carried to value-and-gradient form (invariant 8). h is affine in (b.w - a.w) inside the
    // smoothing band, so differentiating value = lerp(b.w, a.w, h) - k*h*(1-h) with respect to
    // position, the terms carrying dh/dp cancel exactly (dh/dp * [(a.w-b.w) - k + 2*k*h] and the
    // bracket is identically zero given how h is built from a.w-b.w). What is left is exactly
    // lerp(∇b, ∇a, h): the gradient blend is not an approximation, it is the analytic gradient.
    // Outside the band (|a.w - b.w| >= k), h saturates to 0 or 1 and this is exactly min(a, b) with
    // that input's own gradient, continuously (h and the correction term both reach the boundary at
    // the same value from the smooth side, so there is no kink to exclude here the way cellular's
    // F1 = F2 seam needs one).
    public static float4 smin_grad(float4 a, float4 b, float k)
    {
        var h = saturate(0.5f + 0.5f * (b.w - a.w) / k);
        var value = lerp(b.w, a.w, h) - k * h * (1.0f - h);
        var gradient = lerp(new float3(b.x, b.y, b.z), new float3(a.x, a.y, a.z), h);
        return new float4(gradient, value);
    }

    // Worley, "A Cellular Texture Basis Function" (SIGGRAPH 1996): F1/F2 and their analytic
    // gradients, searched over the jittered neighbourhood of feature points. Feature points are
    // selected with pcg3d over each neighbour's integer cell coordinate (design.md, "Integer
    // hashing"), never the float32 bit pattern of that coordinate and never the sin-based hash.
    // Hashing the integer cell coordinate directly is exact and runtime-independent, and it is what
    // the HLSL mirror does, so C# and shader stay bit-for-bit aligned. What pins this is the oracle
    // and mirror tests, not a uniformity measurement: reverting to hashing the float32 bit pattern
    // (or the sin-based hash) is caught by ProductionSearchMatchesTheBruteForceOracleExactly,
    // EveryMirrorFunctionMatchesCSharpMath, F1NeverExceedsF2, and the other tests that key off the
    // integer-hash fixtures. The
    // gradient of a distance field, del|p - c|, is the unit vector (p - c)/|p - c|; that is undefined exactly at
    // a feature point (F1 = 0 or F2 = 0), so cellular follows the repo's existing degenerate-normal
    // convention (GameCult.Geometry.CultGeometryIsoSurface.EmitOrientedTriangle: guard the zero-length
    // case and return the zero vector instead of the NaN a bare normalize would produce there).
    //
    // The search radius is 2 cells on every axis (5x5x5), not the classic 1-cell (3x3x3) radius,
    // because 3x3x3 is not exact: measured 9 of 400,000 points in [-50,50]^3 land on the wrong F2
    // (worst error 0.088), and past |p| of about 2^23 the float32 jitter can put two distinct cells'
    // features on the same value, which 3x3x3 has no chance of noticing either way. 5x5x5 is exact,
    // not merely safer, by the following argument:
    //
    // Let u = p - cell be the query's position inside its own cell (u in [0,1)^3 per component), and
    // let a candidate feature at integer cell offset d = (dx, dy, dz) sit at d + j for some jitter j
    // in [0,1)^3. Per axis, the displacement from p to that feature is (u_i - j_i) - d_i; since
    // u_i - j_i is always strictly inside (-1, 1) regardless of the actual u_i and j_i values,
    // |displacement_i| is always strictly greater than max(0, |d_i| - 1). So every cell with |d_i| >=
    // 3 on some axis is strictly farther than 2 from p (that axis alone already gives max(0, 3-1) =
    // 2). Now bound F1 and F2 from above using two real, always-present candidates: the query's own
    // cell (d = 0), whose distance is always strictly less than sqrt(3) (every axis strictly under
    // 1); and the near-side neighbour on whichever single axis's u_i sits closest to a cell
    // boundary. Let m = min(u_i, 1 - u_i) on that chosen axis (m is in [0, 1/2] by construction,
    // since it is the smaller of the two). Stepping one cell towards that nearer boundary puts the
    // query within (1 + m) of the far face on the offset axis and within (1 - m) of the near face on
    // each of the other two axes, so that neighbour's squared distance is at most
    // (1 + m)^2 + 2*(1 - m)^2 = 3 - 2m + 3m^2, which is <= 3 for every m in [0, 1/2] (3 at m = 0,
    // dipping to a minimum of 8/3 ~= 2.667 at m = 1/3, then rising back to 2.75 at m = 1/2), and
    // strictly below 3 because the jitter never reaches exactly 0 or
    // 1 (u_i - j_i stays strictly inside (-1, 1), so this bound, like the own-cell one, is never
    // tight). F2 is at most the larger of these two real candidates (adding candidates to a set can
    // only lower or hold its 2nd-smallest value), so F1 < sqrt(3) and F2 < sqrt(3) always: roughly
    // 1.732 < 2, margin 2 - sqrt(3) ~= 0.268 (Soul's adversarial check, cut 2a-i: the worst
    // realizable in-box F2 is exactly sqrt(3), reached only in the unattainable limit at a cell
    // corner). Every cell outside the 5x5x5 block is farther than 2, hence farther than either F1 or
    // F2, so the true global F1 and F2 are always inside it. (The matching argument for 3x3x3 only
    // reaches radius 1, whose excluded cells start at distance > 1, which is less than the 1.732
    // ceiling — the gap 3x3x3's failures live in.)
    //
    // An earlier version of this proof picked the near-side neighbour on a fixed axis rather than
    // whichever axis sits closest to its boundary, and that version is wrong (Soul, cut 2a-i):
    // u = (0.5, 0, 0) with the x-neighbour's jittered feature at (-1, 1 - 2^-24, 1 - 2^-24) gives a
    // distance of about 2.06, outside the radius-2 exclusion the proof needs. Choosing the axis by
    // proximity to its own boundary is what keeps m <= 1/2 and the bound at sqrt(3); x was the wrong
    // axis to fix there because u.x = 0.5 is the point on the whole cube farthest from any boundary.
    //
    // The loop below still visits every one of the 125 cells in a fixed order, but skips the
    // pcg3d/jitter/distance work for a cell whose distance cannot possibly beat the current f2: that
    // per-cell lower bound (the max(0, |d_i| - 1) term above, combined across axes) depends only on
    // the integer offset, never on p or the hash, so pruning a cell the unpruned loop would also
    // have rejected changes no output bit (CellularExactSearchTests pins the pruned/unpruned
    // equivalence and the exactness against a 7x7x7 brute-force oracle).
    public static CultCellular cellular(float3 p)
    {
        var cell = floor(p);
        var f1 = 1.0e30f;
        var f2 = 1.0e30f;
        var c1 = new float3(0.0f, 0.0f, 0.0f);
        var c2 = new float3(0.0f, 0.0f, 0.0f);
        var cellId = new float3(0.0f, 0.0f, 0.0f);

        for (var dz = -2; dz <= 2; dz++)
        {
            var lz = max(0.0f, abs((float)dz) - 1.0f);
            for (var dy = -2; dy <= 2; dy++)
            {
                var ly = max(0.0f, abs((float)dy) - 1.0f);
                for (var dx = -2; dx <= 2; dx++)
                {
                    var lx = max(0.0f, abs((float)dx) - 1.0f);
                    var lowerBound = sqrt(lx * lx + ly * ly + lz * lz);
                    if (lowerBound >= f2)
                        continue;

                    var neighbor = cell + new float3(dx, dy, dz);
                    var hash = pcg3d(int3(neighbor));
                    var jitter = new float3(cellular_unit(hash.x), cellular_unit(hash.y), cellular_unit(hash.z));
                    var feature = neighbor + jitter;
                    var d = distance(p, feature);

                    if (d < f1)
                    {
                        f2 = f1; c2 = c1;
                        f1 = d; c1 = feature; cellId = neighbor;
                    }
                    else if (d < f2)
                    {
                        f2 = d; c2 = feature;
                    }
                }
            }
        }

        // F1 = 0 is reachable (p sits exactly on a feature point). F2 = 0 needs two distinct cells'
        // pcg3d-jittered features to land on the same float32 value in every component: unreachable
        // at ordinary magnitudes, but confirmed reachable past roughly |p| = 2^23 (design.md,
        // "Precision domain"), so both share the same degenerate-gradient guard.
        var grad1 = f1 > 0.0f ? (p - c1) / f1 : new float3(0.0f, 0.0f, 0.0f);
        var grad2 = f2 > 0.0f ? (p - c2) / f2 : new float3(0.0f, 0.0f, 0.0f);

        // id comes from pcg4d over the winning cell's integer coordinate alone (same int-hashing
        // rule as the jitter above), never from the pcg3d hash that produced its jitter: sharing a
        // hash between id and jitter would make id a deterministic function of the jitter it is
        // supposed to be independent from, correlated across every cell (F4).
        var idHash = pcg4d(new int4(int3(cellId), 0));

        return new CultCellular(
            new float4(grad1, f1),
            new float4(grad2 - grad1, f2 - f1),
            cellular_unit(idHash.w));
    }

    // Maps a hash component's uint bit pattern to [0, 1) using the top 24 bits, so the result is
    // always strictly less than 1: a plain divide-by-2^32 rounds 0xFFFFFF80..0xFFFFFFFF up to
    // exactly 1.0f in float32, outside the documented range. (uint >> 8) has at most 24 significant
    // bits, which float32's 24-bit mantissa represents exactly, and the division is by an exact power
    // of two, so both the shift-and-convert and the divide are the same IEEE-754 round-to-nearest on
    // dxc and C#, bit-exact on both sides.
    internal static float cellular_unit(int bits) => ((uint)bits >> 8) * (1.0f / 16777216.0f);

    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);
    }
}