Add tests, complete splines, SIMD-accelerate transpose, add a few SIMD methods, add SIMD-based quantization, XorShift constructor, heuristics (incomplete), spline encoder, v256/v128 support to DCT output/source, JxlBlending (blending.cc, blending.h)

src/ImageSharp/Common/Helpers/SimdUtils.cs

@@ -1,4 +1,4 @@
-// Copyright (c) Six Labors.
+// Copyright (c) Six Labors.
 // Licensed under the Six Labors Split License.
 
 using System.Diagnostics;
@@ -71,6 +71,43 @@ internal static partial class SimdUtils
         return val_2p23_f32 | sign;
     }
 
+    /// <summary>
+    /// Estimates the reciprocal of this vector.
+    /// </summary>
+    /// <param name="v">The vector to get reciprocal estimate of.</param>
+    /// <returns>An estimated reciprocal of each element in the vector.</returns>
+    internal static Vector<float> ReciprocalEstimate(this Vector<float> v)
+    {
+        // TODO: System.Runtime.Intrinsics.Arm has Sve and Sve2
+        // support but is for evaluation purposes only; add SVE/SVE2
+        // support when possible
+        if (Avx512F.IsSupported && Vector<float>.Count == 16)
+        {
+            // x86
+            return Avx512F.Reciprocal14(v.AsVector512()).AsVector();
+        }
+        else if (Avx.IsSupported && Vector<float>.Count == 8)
+        {
+            // x86
+            return Avx.Reciprocal(v.AsVector256()).AsVector();
+        }
+        else if (AdvSimd.IsSupported && Vector<float>.Count == 4)
+        {
+            // ARM
+            return AdvSimd.ReciprocalEstimate(v.AsVector128()).AsVector();
+        }
+        else if (Sse.IsSupported && Vector<float>.Count == 4)
+        {
+            // x86
+            return Sse.Reciprocal(v.AsVector128()).AsVector();
+        }
+        else
+        {
+            // Exact reciprocal fallback (slower)
+            return Vector<float>.One / v;
+        }
+    }
+
     [Conditional("DEBUG")]
     private static void DebugVerifySpanInput(ReadOnlySpan<byte> source, ReadOnlySpan<float> dest, int shouldBeDivisibleBy)
     {

src/ImageSharp/Common/Helpers/Vector128Utilities.cs

@@ -899,4 +899,44 @@ internal static class Vector128_
 
         return (shuffledA & maskA) | (shuffledB & maskB);
     }
+
+    /// <summary>
+    /// Interleaves the lower half of the vector.
+    /// </summary>
+    /// <param name="a">First vector</param>
+    /// <param name="b">Second vector</param>
+    /// <returns>
+    ///   <c>{ a[0], b[0], a[1], b[1] }</c>
+    /// </returns>
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static Vector128<float> InterleaveLower(Vector128<float> a, Vector128<float> b)
+    {
+        Vector128<float> shuffledA = Vector128.Shuffle(a, Vector128.Create(0, 0, 1, 1));
+        Vector128<float> shuffledB = Vector128.Shuffle(b, Vector128.Create(0, 0, 1, 1));
+
+        Vector128<float> maskA = Vector128.Create(-1, 0, -1, 0).AsSingle();
+        Vector128<float> maskB = Vector128.Create(0, -1, 0, -1).AsSingle();
+
+        return (shuffledA & maskA) | (shuffledB & maskB);
+    }
+
+    /// <summary>
+    /// Interleaves the upper half of the vector.
+    /// </summary>
+    /// <param name="a">First vector</param>
+    /// <param name="b">Second vector</param>
+    /// <returns>
+    ///   <c>{ a[2], b[2], a[3], b[3] }</c>
+    /// </returns>
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static Vector128<float> InterleaveUpper(Vector128<float> a, Vector128<float> b)
+    {
+        Vector128<float> shuffledA = Vector128.Shuffle(a, Vector128.Create(2, 2, 3, 3));
+        Vector128<float> shuffledB = Vector128.Shuffle(b, Vector128.Create(2, 2, 3, 3));
+
+        Vector128<float> maskA = Vector128.Create(-1, 0, -1, 0).AsSingle();
+        Vector128<float> maskB = Vector128.Create(0, -1, 0, -1).AsSingle();
+
+        return (shuffledA & maskA) | (shuffledB & maskB);
+    }
 }

src/ImageSharp/Common/Helpers/Vector256Utilities.cs

@@ -563,4 +563,44 @@ internal static class Vector256_
 
         return (shuffledA & maskA) | (shuffledB & maskB);
     }
+
+    /// <summary>
+    /// Interleaves the lower half of the vector.
+    /// </summary>
+    /// <param name="a">First vector</param>
+    /// <param name="b">Second vector</param>
+    /// <returns>
+    ///   <c>{ a[0], b[0], a[1], b[1], a[2], b[2], a[3], b[3] }</c>
+    /// </returns>
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static Vector256<float> InterleaveLower(Vector256<float> a, Vector256<float> b)
+    {
+        Vector256<float> shuffledA = Vector256.Shuffle(a, Vector256.Create(0, 0, 1, 1, 2, 2, 3, 3));
+        Vector256<float> shuffledB = Vector256.Shuffle(b, Vector256.Create(0, 0, 1, 1, 2, 2, 3, 3));
+
+        Vector256<float> maskA = Vector256.Create(-1, 0, -1, 0, -1, 0, -1, 0).AsSingle();
+        Vector256<float> maskB = Vector256.Create(0, -1, 0, -1, 0, -1, 0, -1).AsSingle();
+
+        return (shuffledA & maskA) | (shuffledB & maskB);
+    }
+
+    /// <summary>
+    /// Interleaves the upper half of the vector.
+    /// </summary>
+    /// <param name="a">First vector</param>
+    /// <param name="b">Second vector</param>
+    /// <returns>
+    ///   <c>{ a[4], b[4], a[5], b[5], a[6], b[6], a[7], b[7] }</c>
+    /// </returns>
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static Vector256<float> InterleaveUpper(Vector256<float> a, Vector256<float> b)
+    {
+        Vector256<float> shuffledA = Vector256.Shuffle(a, Vector256.Create(4, 4, 5, 5, 6, 6, 7, 7));
+        Vector256<float> shuffledB = Vector256.Shuffle(b, Vector256.Create(4, 4, 5, 5, 6, 6, 7, 7));
+
+        Vector256<float> maskA = Vector256.Create(-1, 0, -1, 0, -1, 0, -1, 0).AsSingle();
+        Vector256<float> maskB = Vector256.Create(0, -1, 0, -1, 0, -1, 0, -1).AsSingle();
+
+        return (shuffledA & maskA) | (shuffledB & maskB);
+    }
 }

src/ImageSharp/Formats/Jxl/Memory/JxlImage3{T}.cs

@@ -1,6 +1,8 @@
 // Copyright (c) Six Labors.
 // Licensed under the Six Labors Split License.
 
+using System.Buffers;
+using System.Runtime.CompilerServices;
 using System.Runtime.InteropServices;
 
 namespace SixLabors.ImageSharp.Formats.Jxl.Memory;
@@ -9,6 +11,21 @@ namespace SixLabors.ImageSharp.Formats.Jxl.Memory;
 internal class JxlImage3<T> : IDisposable
     where T : unmanaged
 {
+    private sealed class TypeChangingMemoryManager<TTarget>(Memory<byte> memory) : MemoryManager<TTarget>
+        where TTarget : unmanaged
+    {
+        public override Span<TTarget> GetSpan() => MemoryMarshal.Cast<byte, TTarget>(memory.Span);
+
+        // we don't use these
+        public override MemoryHandle Pin(int elementIndex = 0) => throw new NotImplementedException();
+
+        public override void Unpin() => throw new NotImplementedException();
+
+        protected override void Dispose(bool disposing)
+        {
+        }
+    }
+
     private const int PlaneCount = 3;
 
     private JxlPlane<T>[] planes = new JxlPlane<T>[3];
@@ -46,6 +63,22 @@ internal class JxlImage3<T> : IDisposable
         return rowSpan;
     }
 
+    // This method performs minor allocations!
+    public Memory<T> PlaneRowMemory(int plane, int row)
+    {
+        this.PlaneRowBoundsCheck(plane, row);
+
+        int rowOffset = row * this.planes[0].BytesPerRow;
+        Memory<byte> rowMemoryBytes = this.planes[plane].Bytes[rowOffset..];
+
+        // we have to allocate a utility class so we can reinterpret
+        // a Memory<T>.
+        // Unsafe.As is truly unsafe because, f.e. what if there are
+        // 400 bytes but T is 4 bytes? the length will remain as 400.
+        TypeChangingMemoryManager<T> reinterpreter = new(rowMemoryBytes);
+        return reinterpreter.Memory;
+    }
+
     public Span<T> PlaneRow(Rectangle rectangle, int c, int y)
     {
         DebugGuard.MustBeGreaterThanOrEqualTo(y + rectangle.Top, 0, nameof(y));

src/ImageSharp/Formats/Jxl/Processing/AcStrategy/JxlAcStrategyRow.cs

@@ -6,17 +6,13 @@ using System.Runtime.InteropServices;
 
 namespace SixLabors.ImageSharp.Formats.Jxl.Processing.AcStrategy;
 
-internal sealed class JxlAcStrategyRow
+internal readonly struct JxlAcStrategyRow(ReadOnlyMemory<byte> row)
 {
-    private readonly ReadOnlyMemory<byte> row;
-
-    public JxlAcStrategyRow(ReadOnlyMemory<byte> row) => this.row = row;
-
-    public JxlAcStrategy this[int x]
+    public readonly JxlAcStrategy this[int x]
     {
         get
         {
-            ReadOnlySpan<byte> span = this.row.Span;
+            ReadOnlySpan<byte> span = row.Span;
 
             DebugGuard.MustBeLessThan(x * 8, span.Length, "x overflows");
 

src/ImageSharp/Formats/Jxl/Processing/Blending/JxlAlphaBlendingInputLayer.cs

@@ -3,13 +3,13 @@
 
 namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Blending;
 
-internal sealed class JxlAlphaBlendingInputLayer
+internal ref struct JxlAlphaBlendingInputLayer(ReadOnlySpan<float> singleSpan)
 {
-    public ReadOnlyMemory<float> R { get; set; }
+    public ReadOnlySpan<float> R = singleSpan;
 
-    public ReadOnlyMemory<float> G { get; set; }
+    public ReadOnlySpan<float> G = singleSpan;
 
-    public ReadOnlyMemory<float> B { get; set; }
+    public ReadOnlySpan<float> B = singleSpan;
 
-    public ReadOnlyMemory<float> A { get; set; }
+    public ReadOnlySpan<float> A = singleSpan;
 }

src/ImageSharp/Formats/Jxl/Processing/Blending/JxlAlphaBlendingOutput.cs

@@ -3,13 +3,13 @@
 
 namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Blending;
 
-internal sealed class JxlAlphaBlendingOutput
+internal ref struct JxlAlphaBlendingOutput(Span<float> singleSpan)
 {
-    public Memory<float> R { get; set; }
+    public Span<float> R = singleSpan;
 
-    public Memory<float> G { get; set; }
+    public Span<float> G = singleSpan;
 
-    public Memory<float> B { get; set; }
+    public Span<float> B = singleSpan;
 
-    public Memory<float> A { get; set; }
+    public Span<float> A = singleSpan;
 }

src/ImageSharp/Formats/Jxl/Processing/Blending/JxlBlending.cs

@@ -0,0 +1,368 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using SixLabors.ImageSharp.Formats.Jxl.IO.FrameHeader;
+using SixLabors.ImageSharp.Formats.Jxl.IO.Metadata;
+using SixLabors.ImageSharp.Formats.Jxl.Memory.ImageTypes;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+using SixLabors.ImageSharp.Memory;
+
+namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Blending;
+
+internal static class JxlBlending
+{
+    public static bool NeedsBlending(JxlFrameHeader header)
+    {
+        if (header.FrameType is not JxlFrameType.RegularFrame and not JxlFrameType.SkipProgressive)
+        {
+            return false;
+        }
+
+        JxlBlendingInfo? blendingInfo = header.BlendingInfo;
+        if (blendingInfo is null)
+        {
+            return false;
+        }
+
+        bool replaceAll = blendingInfo.BlendMode == JxlBlendMode.Replace;
+
+        foreach (JxlBlendingInfo info in header.ExtraChannelBlendingInfo)
+        {
+            if (info.BlendMode != JxlBlendMode.Replace)
+            {
+                replaceAll = false;
+            }
+        }
+
+        if (!header.CustomSizeOrOrigin && replaceAll)
+        {
+            return false;
+        }
+
+        return true;
+    }
+
+    public static void PerformBlending(
+        Configuration configuration,
+        Buffer2D<float> bg,
+        Buffer2D<float> fg,
+        Buffer2D<float> output,
+        int x0,
+        int xsize,
+        JxlPatchBlending colorBlending,
+        Span<JxlPatchBlending> ecBlending,
+        List<JxlExtraChannelInfo> extraChannelInfo)
+    {
+        bool hasAlpha = extraChannelInfo.Any(x => x.Type == JxlExtraChannel.Alpha);
+
+        int numEc = extraChannelInfo.Count;
+        using JxlImageF tmp = new(configuration, xsize, 3 + numEc);
+
+        for (int i = 0; i < numEc; i++)
+        {
+            int i3 = 3 + i;
+
+            switch (ecBlending[i].Mode)
+            {
+                case JxlPatchBlendMode.Add:
+                {
+                    Span<float> row = tmp.GetRow(i3);
+                    for (int x = 0; x < xsize; x++)
+                    {
+                        row[x] = bg[i3, x + x0] + fg[i3, x + x0];
+                    }
+
+                    continue;
+                }
+
+                case JxlPatchBlendMode.BlendAbove:
+                {
+                    int alpha = ecBlending[i].AlphaChannel;
+                    bool isPremultiplied = extraChannelInfo[alpha].AlphaAssociated;
+
+                    Span<float> bgSpan3 = bg.DangerousGetRowSpan(i3)[x0..];
+                    Span<float> bgSpan3Alpha = bg.DangerousGetRowSpan(3 + alpha)[x0..];
+                    Span<float> fgSpan3 = fg.DangerousGetRowSpan(i3)[x0..];
+                    Span<float> fgSpan3Alpha = fg.DangerousGetRowSpan(3 + alpha)[x0..];
+
+                    JxlAlphaHelper.PerformAlphaBlending(
+                        bgSpan3,
+                        bgSpan3Alpha,
+                        fgSpan3,
+                        fgSpan3Alpha,
+                        tmp.GetRow(i3),
+                        xsize,
+                        isPremultiplied,
+                        ecBlending[i].Clamp);
+
+                    continue;
+                }
+
+                case JxlPatchBlendMode.BlendBelow:
+                {
+                    int alpha = ecBlending[i].AlphaChannel;
+                    bool isPremultiplied = extraChannelInfo[alpha].AlphaAssociated;
+
+                    Span<float> bgSpan3 = bg.DangerousGetRowSpan(i3)[x0..];
+                    Span<float> bgSpan3Alpha = bg.DangerousGetRowSpan(3 + alpha)[x0..];
+                    Span<float> fgSpan3 = fg.DangerousGetRowSpan(i3)[x0..];
+                    Span<float> fgSpan3Alpha = fg.DangerousGetRowSpan(3 + alpha)[x0..];
+
+                    JxlAlphaHelper.PerformAlphaBlending(
+                        bgSpan3,
+                        bgSpan3Alpha,
+                        fgSpan3,
+                        fgSpan3Alpha,
+                        tmp.GetRow(3 + i),
+                        xsize,
+                        isPremultiplied,
+                        ecBlending[i].Clamp);
+
+                    continue;
+                }
+
+                case JxlPatchBlendMode.AlphaWeightedAddAbove:
+                {
+                    int alpha = ecBlending[i].AlphaChannel;
+
+                    Span<float> bgSpan3 = bg.DangerousGetRowSpan(i3)[x0..];
+                    Span<float> bgSpan3Alpha = bg.DangerousGetRowSpan(3 + alpha)[x0..];
+                    Span<float> fgSpan3 = fg.DangerousGetRowSpan(i3)[x0..];
+
+                    JxlAlphaHelper.PerformAlphaWeightedAdd(
+                        bgSpan3,
+                        fgSpan3,
+                        bgSpan3Alpha,
+                        tmp.GetRow(3 + i),
+                        xsize,
+                        ecBlending[i].Clamp);
+
+                    continue;
+                }
+
+                case JxlPatchBlendMode.AlphaWeightedAddBelow:
+                {
+                    int alpha = ecBlending[i].AlphaChannel;
+
+                    Span<float> bgSpan3 = bg.DangerousGetRowSpan(i3)[x0..];
+                    Span<float> bgSpan3Alpha = bg.DangerousGetRowSpan(3 + alpha)[x0..];
+                    Span<float> fgSpan3 = fg.DangerousGetRowSpan(i3)[x0..];
+
+                    JxlAlphaHelper.PerformAlphaWeightedAdd(
+                        fgSpan3,
+                        bgSpan3,
+                        bgSpan3Alpha,
+                        tmp.GetRow(3 + i),
+                        xsize,
+                        ecBlending[i].Clamp);
+
+                    continue;
+                }
+
+                case JxlPatchBlendMode.Multiply:
+                {
+                    Span<float> bgSpan3 = bg.DangerousGetRowSpan(i3)[x0..];
+                    Span<float> fgSpan3 = fg.DangerousGetRowSpan(i3)[x0..];
+
+                    JxlAlphaHelper.PerformMultiplyBlending(
+                        bgSpan3,
+                        fgSpan3,
+                        tmp.GetRow(i3),
+                        xsize,
+                        ecBlending[i].Clamp);
+
+                    continue;
+                }
+
+                case JxlPatchBlendMode.Replace:
+                    if (xsize > 0)
+                    {
+                        Span<float> fgSpan3 = fg.DangerousGetRowSpan(i3)[x0..];
+                        fgSpan3.Slice(0, xsize).CopyTo(tmp.GetRow(i3));
+                    }
+
+                    continue;
+
+                case JxlPatchBlendMode.None:
+                    if (xsize > 0)
+                    {
+                        Span<float> bgSpan3 = bg.DangerousGetRowSpan(i3)[x0..];
+                        bgSpan3.Slice(0, xsize).CopyTo(tmp.GetRow(i3));
+                    }
+
+                    continue;
+            }
+        }
+
+        int colorBlendingAlpha = colorBlending.AlphaChannel;
+
+        void Add()
+        {
+            for (int p = 0; p < 3; p++)
+            {
+                Span<float> output = tmp.GetRow(p);
+                Span<float> bgSpan = bg.DangerousGetRowSpan(p);
+                Span<float> fgSpan = fg.DangerousGetRowSpan(p);
+
+                for (int x = 0; x < xsize; x++)
+                {
+                    int xPlusX0 = x + x0;
+
+                    output[x] = bgSpan[xPlusX0] + fgSpan[xPlusX0];
+                }
+            }
+        }
+
+        void BlendWeighted(Span<float> bottom, Span<float> top)
+        {
+            bool isPremultiplied = extraChannelInfo[colorBlendingAlpha].AlphaAssociated;
+
+            JxlAlphaHelper.PerformAlphaBlending(
+                new JxlAlphaBlendingInputLayer()
+                {
+                    R = bottom[x0..],
+                    G = bottom[(x0 + 1)..],
+                    B = bottom[(2 + x0)..],
+                    A = bottom[(3 + colorBlendingAlpha + x0)..]
+                },
+                new JxlAlphaBlendingInputLayer()
+                {
+                    R = top[x0..],
+                    G = top[(x0 + 1)..],
+                    B = top[(x0 + 2)..],
+                    A = top[(3 + colorBlendingAlpha + x0)..]
+                },
+                new JxlAlphaBlendingOutput()
+                {
+                    R = tmp.GetRow(0),
+                    G = tmp.GetRow(1),
+                    B = tmp.GetRow(2),
+                    A = tmp.GetRow(3)
+                },
+                xsize,
+                isPremultiplied,
+                colorBlending.Clamp);
+        }
+
+        void AddWeighted(Span<float> bottom, Span<float> top)
+        {
+            for (int c = 0; c < 3; c++)
+            {
+                JxlAlphaHelper.PerformAlphaWeightedAdd(bottom[(c + x0)..], top[(c + x0)..], top[(3 + colorBlendingAlpha + x0)..], tmp.GetRow(c), xsize, colorBlending.Clamp);
+            }
+        }
+
+        void Copy(Span<float> src)
+        {
+            for (int p = 0; p < 3; p++)
+            {
+                src.Slice(p + x0, xsize).CopyTo(tmp.GetRow(p));
+            }
+        }
+
+        switch (colorBlending.Mode)
+        {
+            case JxlPatchBlendMode.Add:
+            {
+                Add();
+                break;
+            }
+
+            case JxlPatchBlendMode.AlphaWeightedAddAbove:
+            {
+                if (hasAlpha)
+                {
+                    AddWeighted(bg.DangerousGetSingleSpan(), fg.DangerousGetSingleSpan());
+                }
+                else
+                {
+                    Add();
+                }
+
+                break;
+            }
+
+            case JxlPatchBlendMode.AlphaWeightedAddBelow:
+            {
+                if (hasAlpha)
+                {
+                    AddWeighted(fg.DangerousGetSingleSpan(), bg.DangerousGetSingleSpan());
+                }
+                else
+                {
+                    Add();
+                }
+
+                break;
+            }
+
+            case JxlPatchBlendMode.BlendAbove:
+            {
+                if (hasAlpha)
+                {
+                    BlendWeighted(bg.DangerousGetSingleSpan(), fg.DangerousGetSingleSpan());
+                }
+                else
+                {
+                    Copy(fg.DangerousGetSingleSpan());
+                }
+
+                break;
+            }
+
+            case JxlPatchBlendMode.BlendBelow:
+            {
+                if (hasAlpha)
+                {
+                    BlendWeighted(fg.DangerousGetSingleSpan(), bg.DangerousGetSingleSpan());
+                }
+                else
+                {
+                    Copy(fg.DangerousGetSingleSpan());
+                }
+
+                break;
+            }
+
+            case JxlPatchBlendMode.Multiply:
+            {
+                Span<float> bgSpan = bg.DangerousGetSingleSpan();
+                Span<float> fgSpan = fg.DangerousGetSingleSpan();
+
+                for (int p = 0; p < 3; p++)
+                {
+                    JxlAlphaHelper.PerformMultiplyBlending(
+                        bgSpan[(p + x0)..],
+                        fgSpan[(p + x0)..],
+                        tmp.GetRow(p),
+                        xsize,
+                        colorBlending.Clamp);
+                }
+
+                break;
+            }
+
+            case JxlPatchBlendMode.Replace:
+            {
+                Copy(fg.DangerousGetSingleSpan());
+                break;
+            }
+
+            case JxlPatchBlendMode.None:
+            {
+                Copy(bg.DangerousGetSingleSpan());
+                break;
+            }
+        }
+
+        if (xsize != 0)
+        {
+            Span<float> outputSpan = output.DangerousGetSingleSpan();
+
+            for (int i = 0; i < 3; i++)
+            {
+                tmp.GetRow(i).Slice(0, xsize).CopyTo(outputSpan[(i + x0)..]);
+            }
+        }
+    }
+}

src/ImageSharp/Formats/Jxl/Processing/Dct/JxlDctOutput.cs

@@ -3,6 +3,7 @@
 
 using System.Numerics;
 using System.Runtime.CompilerServices;
+using System.Runtime.Intrinsics;
 
 namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Dct;
 
@@ -49,4 +50,22 @@ internal ref struct JxlDctOutput(Span<float> data, int stride)
     /// <param name="index">The offset.</param>
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     public readonly void StorePart(Vector<float> value, int row, int index) => value.CopyTo(this.Address(row, index));
+
+    /// <summary>
+    /// Stores the vector into the data at the specified row and offset.
+    /// </summary>
+    /// <param name="value">The vector to write.</param>
+    /// <param name="row">The row index.</param>
+    /// <param name="index">The offset.</param>
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public readonly void StorePart256(Vector256<float> value, int row, int index) => value.CopyTo(this.Address(row, index));
+
+    /// <summary>
+    /// Stores the vector into the data at the specified row and offset.
+    /// </summary>
+    /// <param name="value">The vector to write.</param>
+    /// <param name="row">The row index.</param>
+    /// <param name="index">The offset.</param>
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public readonly void StorePart128(Vector128<float> value, int row, int index) => value.CopyTo(this.Address(row, index));
 }

src/ImageSharp/Formats/Jxl/Processing/Dct/JxlDctSource.cs

@@ -3,6 +3,7 @@
 
 using System.Numerics;
 using System.Runtime.CompilerServices;
+using System.Runtime.Intrinsics;
 
 namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Dct;
 
@@ -52,4 +53,24 @@ internal readonly ref struct JxlDctSource(Span<float> data, int stride)
     /// Vector at that row and offset.
     /// </returns>
     public Vector<float> LoadPart(int row, int i) => new(this.Address(row, i));
+
+    /// <summary>
+    /// Loads a vector at the specified row and offset.
+    /// </summary>
+    /// <param name="row">The row index.</param>
+    /// <param name="i">The offset.</param>
+    /// <returns>
+    /// Vector at that row and offset.
+    /// </returns>
+    public Vector256<float> LoadPart256(int row, int i) => Vector256.Create<float>(this.Address(row, i));
+
+    /// <summary>
+    /// Loads a vector at the specified row and offset.
+    /// </summary>
+    /// <param name="row">The row index.</param>
+    /// <param name="i">The offset.</param>
+    /// <returns>
+    /// Vector at that row and offset.
+    /// </returns>
+    public Vector128<float> LoadPart128(int row, int i) => Vector128.Create<float>(this.Address(row, i));
 }

src/ImageSharp/Formats/Jxl/Processing/Encoder/Ans/JxlToken.cs

@@ -0,0 +1,13 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Splines;
+
+namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Encoder.Ans;
+
+internal struct JxlToken(JxlSplineEntropyContext c, uint value)
+{
+    public bool IsLz77Length;
+    public JxlSplineEntropyContext Context = c;
+    public uint Value = value;
+}

src/ImageSharp/Formats/Jxl/Processing/Encoder/JxlHeuristics.cs

@@ -0,0 +1,170 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using System.Buffers;
+using System.Runtime.InteropServices;
+using SixLabors.ImageSharp.Formats.Jxl.Memory.ImageTypes;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.AcStrategy;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Primitives;
+
+namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Encoder;
+
+internal static class JxlHeuristics
+{
+    private static ReadOnlySpan<byte> SimpleContextMap =>
+    [
+        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
+    ];
+
+    public static void FindBestBlockEntropyModel(JxlCompressParameters cparameters, JxlImageI rqf, JxlAcStrategyImage acStrategy, JxlBlockContextMap blockCtxMap)
+    {
+        if (cparameters.DecodingSpeedTier >= 1)
+        {
+            SimpleContextMap.CopyTo(blockCtxMap.ContextMap.AsSpan());
+            blockCtxMap.ContextCount = 2;
+            blockCtxMap.DcContextCount = 1;
+            return;
+        }
+
+        if (cparameters.SpeedTier >= JxlSpeedTier.Falcon)
+        {
+            return;
+        }
+
+        int total = rqf.XSize * rqf.YSize;
+        int sizeForContextModel = (1 << 10) * cparameters.ButteraugliDistance;
+
+        if (total < sizeForContextModel)
+        {
+            return;
+        }
+
+        OccCounters counters = new(rqf, acStrategy);
+        int sizeForQfSplit = (1 << 13) * cparameters.ButteraugliDistance;
+        int numQfSegments = total < sizeForQfSplit ? 1 : 2;
+        List<uint> qft = blockCtxMap.QfThresholds;
+        qft.Clear();
+        int cumulativeSum = 0;
+        int next = 1;
+        int lastCut = 256;
+        int cut = total * next / numQfSegments;
+
+        for (int j = 0; j < 256; j++)
+        {
+            cumulativeSum += counters.QfCounts[j];
+
+            if (cumulativeSum > cut)
+            {
+                if (j != 0)
+                {
+                    qft.Add((uint)j);
+                }
+
+                lastCut = j;
+
+                while (cumulativeSum > cut)
+                {
+                    next++;
+                    cut = total * next / numQfSegments;
+                }
+            }
+            else if (next > qft.Count + 1)
+            {
+                if (j - 1 == lastCut && j != 0)
+                {
+                    qft.Add((uint)j);
+                }
+            }
+        }
+
+        int[]? pooledCounts = null;
+        int[]? pooledRemap = null;
+        int[]? pooledClusters = null;
+        int countsLength = JxlForwardCoefficientOrder.OrderCount * (qft.Count + 1);
+
+        Span<int> counts =
+            countsLength <= 128
+                ? stackalloc int[128].Slice(0, countsLength)
+                : pooledCounts = ArrayPool<int>.Shared.Rent(countsLength);
+
+        Span<int> remap =
+            countsLength <= 128
+                ? stackalloc int[128].Slice(0, countsLength)
+                : pooledRemap = ArrayPool<int>.Shared.Rent(countsLength);
+
+        Span<int> clusters =
+            countsLength <= 128
+                ? stackalloc int[128].Slice(0, countsLength)
+                : pooledClusters = ArrayPool<int>.Shared.Rent(countsLength);
+
+        int qftPos = 0;
+
+        for (int j = 0; j < 256; j++)
+        {
+            if (qftPos < qft.Count && j == qft[qftPos])
+            {
+                qftPos++;
+            }
+
+            for (int i = 0; i < JxlForwardCoefficientOrder.OrderCount; i++)
+            {
+                counts[qftPos + (i * (qft.Count + 1))] += counters.QfOrdCounts[i, j];
+            }
+        }
+
+        JxlSimdUtils.Iota(remap, 0);
+        remap.CopyTo(clusters);
+
+        int numClusters = Math.Clamp(total / sizeForContextModel / 2, 2, 9);
+        int numClustersChroma = Math.Clamp(total / sizeForContextModel / 3, 1, 5);
+
+        // TODO: method incomplete
+        // do not forget to ArrayPool<int>.Shared.Return pooledCounts, pooledRemap, pooledClusters if needed
+    }
+
+    private sealed class OccCounters : IDisposable
+    {
+        private readonly int[] qfCounts;
+        private readonly int[] dataForQfOrdCounts;
+        private readonly int[] ordCounts;
+
+        public OccCounters(JxlImageI rqf, JxlAcStrategyImage acStrategy)
+        {
+            this.qfCounts = ArrayPool<int>.Shared.Rent(256);
+            this.dataForQfOrdCounts = ArrayPool<int>.Shared.Rent(256 * JxlForwardCoefficientOrder.OrderCount);
+            this.ordCounts = ArrayPool<int>.Shared.Rent(JxlForwardCoefficientOrder.OrderCount);
+
+            this.QfOrdCounts = new(JxlForwardCoefficientOrder.OrderCount, 256, this.dataForQfOrdCounts);
+
+            for (int y = 0; y < rqf.YSize; y++)
+            {
+                Span<int> qfRow = rqf.GetRow(y);
+                JxlAcStrategyRow acsRow = acStrategy.GetRow(y);
+
+                for (int x = 0; x < rqf.XSize; x++)
+                {
+                    int ord = JxlCoefficientOrder.StrategyOrder[acsRow[x].RawStrategy];
+                    int qf = qfRow[x] - 1;
+                    this.qfCounts[qf]++;
+                    this.QfOrdCounts[ord, qf]++;
+                    this.ordCounts[ord]++;
+                }
+            }
+        }
+
+        public Span<int> QfCounts => this.qfCounts.AsSpan();
+
+        public DenseMatrix<int> QfOrdCounts { get; }
+
+        public Span<int> OrdCounts => this.ordCounts.AsSpan();
+
+        public void Dispose()
+        {
+            ArrayPool<int>.Shared.Return(this.qfCounts);
+            ArrayPool<int>.Shared.Return(this.dataForQfOrdCounts);
+            ArrayPool<int>.Shared.Return(this.ordCounts);
+        }
+    }
+}

src/ImageSharp/Formats/Jxl/Processing/Encoder/JxlSplineEncoder.cs

@@ -0,0 +1,83 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Encoder.Ans;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Encoder.AuxiliaryOutput;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Primitives;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Splines;
+
+namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Encoder;
+
+internal sealed class JxlSplineEncoder
+{
+    private static void Tokenize(JxlQuantizedSpline spline, List<JxlToken> tokens)
+    {
+        tokens.Add(new(JxlSplineEntropyContext.NumControlPoints, (uint)spline.ControlPoints.Length));
+
+        foreach (JxlControlPoint point in spline.ControlPoints.Span)
+        {
+            tokens.Add(new(JxlSplineEntropyContext.ControlPoints, JxlPackSigned.PackUnsigned(point.First)));
+            tokens.Add(new(JxlSplineEntropyContext.ControlPoints, JxlPackSigned.PackUnsigned(point.Second)));
+        }
+
+        void EncodeDCT(Span<int> dct)
+        {
+            for (int i = 0; i < 32; i++)
+            {
+                tokens.Add(new(JxlSplineEntropyContext.Dct, JxlPackSigned.PackUnsigned(dct[i])));
+            }
+        }
+
+        foreach (Span<int> dct in spline.ColorDct)
+        {
+            EncodeDCT(dct);
+        }
+
+        EncodeDCT(spline.SigmaDct);
+    }
+
+    public static void EncodeAllStartingPoints(Span<PointF> points, List<JxlToken> tokens)
+    {
+        long lastX = 0;
+        long lastY = 0;
+
+        for (int i = 0; i < points.Length; i++)
+        {
+            long x = (long)MathF.Round(points[i].X, MidpointRounding.AwayFromZero);
+            long y = (long)MathF.Round(points[i].Y, MidpointRounding.AwayFromZero);
+
+            if (i == 0)
+            {
+                tokens.Add(new(JxlSplineEntropyContext.StartingPosition, (uint)x));
+                tokens.Add(new(JxlSplineEntropyContext.StartingPosition, (uint)y));
+            }
+            else
+            {
+                tokens.Add(new(JxlSplineEntropyContext.StartingPosition, JxlPackSigned.PackUnsigned((int)(x - lastX))));
+                tokens.Add(new(JxlSplineEntropyContext.StartingPosition, JxlPackSigned.PackUnsigned((int)(y - lastY))));
+            }
+
+            lastX = x;
+            lastY = y;
+        }
+    }
+
+    public static void EncodeSplines(JxlSplines splines, JxlBitWriter writer, JxlLayerType layer, JxlHistogramParameters histogramParameters, JxlAuxiliaryOutput auxOut)
+    {
+        Span<JxlQuantizedSpline> quantizedSplines = splines.QuantizedSplines;
+        List<List<JxlToken>> tokens = [[]];
+        tokens[0].Add(new(JxlSplineEntropyContext.NumSplineContexts, (uint)(quantizedSplines.Length - 1)));
+
+        EncodeAllStartingPoints(splines.StartingPoints, tokens[0]);
+
+        tokens[0].Add(new(JxlSplineEntropyContext.QuantizationAdjustment, JxlPackSigned.PackUnsigned(splines.QuantizationAdjustment)));
+
+        foreach (JxlQuantizedSpline spline in quantizedSplines)
+        {
+            Tokenize(spline, tokens[0]);
+        }
+
+        _ = BuildAndEncodeHistograms(writer, histogramParameters, JxlSplineEntropyContext.NumSplineContexts, tokens, out JxlEntropyEncodingData codes, writer, layer, auxOut);
+        WriteTokens(tokens[0], codes, 0, writer, layer, auxOut);
+    }
+}

src/ImageSharp/Formats/Jxl/Processing/JxlCoefficientOrder.cs

@@ -50,7 +50,7 @@ internal static class JxlCoefficientOrder
         return Math.Min(token, PermutationContexts - 1u);
     }
 
-    public static bool ReadPermutation(int skip, int size, Span<int> order, JxlBitReader bitReader, JxlAnsSymbolReader reader, Span<byte> contextMap)
+    public static bool ReadPermutation(int skip, int size, Span<uint> order, JxlBitReader bitReader, JxlAnsSymbolReader reader, Span<byte> contextMap)
     {
         DebugGuard.MustBeLessThanOrEqualTo(size, 65536, nameof(size));
 

src/ImageSharp/Formats/Jxl/Processing/JxlSimdUtils.StoreInterleaved.Generated.cs

@@ -10,124 +10,139 @@ namespace SixLabors.ImageSharp.Formats.Jxl.Processing;
 
 internal static partial class JxlSimdUtils
 {
-    public static void StoreInterleaved<T>(Vector<T> v1, Vector<T> v2, ref T memory)
+    public static unsafe void StoreInterleaved<T>(Vector<T> v1, Vector<T> v2, ref T memory)
+        where T : unmanaged
     {
-        v1.StoreUnsafe(ref Unsafe.Add(ref memory, 0));
-        v2.StoreUnsafe(ref Unsafe.Add(ref memory, 1));
+        v1.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 0));
+        v2.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 1));
     }
 
-    public static void StoreInterleaved<T>(Vector<T> v1, Vector<T> v2, Vector<T> v3, ref T memory)
+    public static unsafe void StoreInterleaved<T>(Vector<T> v1, Vector<T> v2, Vector<T> v3, ref T memory)
+        where T : unmanaged
     {
-        v1.StoreUnsafe(ref Unsafe.Add(ref memory, 0));
-        v2.StoreUnsafe(ref Unsafe.Add(ref memory, 1));
-        v3.StoreUnsafe(ref Unsafe.Add(ref memory, 2));
+        v1.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 0));
+        v2.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 1));
+        v3.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 2));
     }
 
-    public static void StoreInterleaved<T>(Vector<T> v1, Vector<T> v2, Vector<T> v3, Vector<T> v4, ref T memory)
+    public static unsafe void StoreInterleaved<T>(Vector<T> v1, Vector<T> v2, Vector<T> v3, Vector<T> v4, ref T memory)
+        where T : unmanaged
     {
-        v1.StoreUnsafe(ref Unsafe.Add(ref memory, 0));
-        v2.StoreUnsafe(ref Unsafe.Add(ref memory, 1));
-        v3.StoreUnsafe(ref Unsafe.Add(ref memory, 2));
-        v4.StoreUnsafe(ref Unsafe.Add(ref memory, 3));
+        v1.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 0));
+        v2.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 1));
+        v3.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 2));
+        v4.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 3));
     }
 
-    public static void StoreInterleaved<T>(Vector<T> v1, Vector<T> v2, Vector<T> v3, Vector<T> v4, Vector<T> v5, ref T memory)
+    public static unsafe void StoreInterleaved<T>(Vector<T> v1, Vector<T> v2, Vector<T> v3, Vector<T> v4, Vector<T> v5, ref T memory)
+        where T : unmanaged
     {
-        v1.StoreUnsafe(ref Unsafe.Add(ref memory, 0));
-        v2.StoreUnsafe(ref Unsafe.Add(ref memory, 1));
-        v3.StoreUnsafe(ref Unsafe.Add(ref memory, 2));
-        v4.StoreUnsafe(ref Unsafe.Add(ref memory, 3));
-        v5.StoreUnsafe(ref Unsafe.Add(ref memory, 4));
+        v1.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 0));
+        v2.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 1));
+        v3.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 2));
+        v4.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 3));
+        v5.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 4));
     }
 
-    public static void StoreInterleaved<T>(Vector<T> v1, Vector<T> v2, Vector<T> v3, Vector<T> v4, Vector<T> v5, Vector<T> v6, ref T memory)
+    public static unsafe void StoreInterleaved<T>(Vector<T> v1, Vector<T> v2, Vector<T> v3, Vector<T> v4, Vector<T> v5, Vector<T> v6, ref T memory)
+        where T : unmanaged
     {
-        v1.StoreUnsafe(ref Unsafe.Add(ref memory, 0));
-        v2.StoreUnsafe(ref Unsafe.Add(ref memory, 1));
-        v3.StoreUnsafe(ref Unsafe.Add(ref memory, 2));
-        v4.StoreUnsafe(ref Unsafe.Add(ref memory, 3));
-        v5.StoreUnsafe(ref Unsafe.Add(ref memory, 4));
-        v6.StoreUnsafe(ref Unsafe.Add(ref memory, 5));
+        v1.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 0));
+        v2.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 1));
+        v3.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 2));
+        v4.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 3));
+        v5.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 4));
+        v6.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 5));
     }
 
-    public static void StoreInterleaved<T>(Vector128<T> v1, Vector128<T> v2, ref T memory)
+    public static unsafe void StoreInterleaved<T>(Vector128<T> v1, Vector128<T> v2, ref T memory)
+        where T : unmanaged
     {
-        v1.StoreUnsafe(ref Unsafe.Add(ref memory, 0));
-        v2.StoreUnsafe(ref Unsafe.Add(ref memory, 1));
+        v1.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 0));
+        v2.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 1));
     }
 
-    public static void StoreInterleaved<T>(Vector128<T> v1, Vector128<T> v2, Vector128<T> v3, ref T memory)
+    public static unsafe void StoreInterleaved<T>(Vector128<T> v1, Vector128<T> v2, Vector128<T> v3, ref T memory)
+        where T : unmanaged
     {
-        v1.StoreUnsafe(ref Unsafe.Add(ref memory, 0));
-        v2.StoreUnsafe(ref Unsafe.Add(ref memory, 1));
-        v3.StoreUnsafe(ref Unsafe.Add(ref memory, 2));
+        v1.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 0));
+        v2.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 1));
+        v3.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 2));
     }
 
-    public static void StoreInterleaved<T>(Vector128<T> v1, Vector128<T> v2, Vector128<T> v3, Vector128<T> v4, ref T memory)
+    public static unsafe void StoreInterleaved<T>(Vector128<T> v1, Vector128<T> v2, Vector128<T> v3, Vector128<T> v4, ref T memory)
+        where T : unmanaged
     {
-        v1.StoreUnsafe(ref Unsafe.Add(ref memory, 0));
-        v2.StoreUnsafe(ref Unsafe.Add(ref memory, 1));
-        v3.StoreUnsafe(ref Unsafe.Add(ref memory, 2));
-        v4.StoreUnsafe(ref Unsafe.Add(ref memory, 3));
+        v1.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 0));
+        v2.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 1));
+        v3.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 2));
+        v4.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 3));
     }
 
-    public static void StoreInterleaved<T>(Vector128<T> v1, Vector128<T> v2, Vector128<T> v3, Vector128<T> v4, Vector128<T> v5, ref T memory)
+    public static unsafe void StoreInterleaved<T>(Vector128<T> v1, Vector128<T> v2, Vector128<T> v3, Vector128<T> v4, Vector128<T> v5, ref T memory)
+        where T : unmanaged
     {
-        v1.StoreUnsafe(ref Unsafe.Add(ref memory, 0));
-        v2.StoreUnsafe(ref Unsafe.Add(ref memory, 1));
-        v3.StoreUnsafe(ref Unsafe.Add(ref memory, 2));
-        v4.StoreUnsafe(ref Unsafe.Add(ref memory, 3));
-        v5.StoreUnsafe(ref Unsafe.Add(ref memory, 4));
+        v1.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 0));
+        v2.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 1));
+        v3.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 2));
+        v4.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 3));
+        v5.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 4));
     }
 
-    public static void StoreInterleaved<T>(Vector128<T> v1, Vector128<T> v2, Vector128<T> v3, Vector128<T> v4, Vector128<T> v5, Vector128<T> v6, ref T memory)
+    public static unsafe void StoreInterleaved<T>(Vector128<T> v1, Vector128<T> v2, Vector128<T> v3, Vector128<T> v4, Vector128<T> v5, Vector128<T> v6, ref T memory)
+        where T : unmanaged
     {
-        v1.StoreUnsafe(ref Unsafe.Add(ref memory, 0));
-        v2.StoreUnsafe(ref Unsafe.Add(ref memory, 1));
-        v3.StoreUnsafe(ref Unsafe.Add(ref memory, 2));
-        v4.StoreUnsafe(ref Unsafe.Add(ref memory, 3));
-        v5.StoreUnsafe(ref Unsafe.Add(ref memory, 4));
-        v6.StoreUnsafe(ref Unsafe.Add(ref memory, 5));
+        v1.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 0));
+        v2.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 1));
+        v3.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 2));
+        v4.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 3));
+        v5.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 4));
+        v6.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 5));
     }
 
-    public static void StoreInterleaved<T>(Vector256<T> v1, Vector256<T> v2, ref T memory)
+    public static unsafe void StoreInterleaved<T>(Vector256<T> v1, Vector256<T> v2, ref T memory)
+        where T : unmanaged
     {
-        v1.StoreUnsafe(ref Unsafe.Add(ref memory, 0));
-        v2.StoreUnsafe(ref Unsafe.Add(ref memory, 1));
+        v1.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 0));
+        v2.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 1));
     }
 
-    public static void StoreInterleaved<T>(Vector256<T> v1, Vector256<T> v2, Vector256<T> v3, ref T memory)
+    public static unsafe void StoreInterleaved<T>(Vector256<T> v1, Vector256<T> v2, Vector256<T> v3, ref T memory)
+        where T : unmanaged
     {
-        v1.StoreUnsafe(ref Unsafe.Add(ref memory, 0));
-        v2.StoreUnsafe(ref Unsafe.Add(ref memory, 1));
-        v3.StoreUnsafe(ref Unsafe.Add(ref memory, 2));
+        v1.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 0));
+        v2.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 1));
+        v3.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 2));
     }
 
-    public static void StoreInterleaved<T>(Vector256<T> v1, Vector256<T> v2, Vector256<T> v3, Vector256<T> v4, ref T memory)
+    public static unsafe void StoreInterleaved<T>(Vector256<T> v1, Vector256<T> v2, Vector256<T> v3, Vector256<T> v4, ref T memory)
+        where T : unmanaged
     {
-        v1.StoreUnsafe(ref Unsafe.Add(ref memory, 0));
-        v2.StoreUnsafe(ref Unsafe.Add(ref memory, 1));
-        v3.StoreUnsafe(ref Unsafe.Add(ref memory, 2));
-        v4.StoreUnsafe(ref Unsafe.Add(ref memory, 3));
+        v1.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 0));
+        v2.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 1));
+        v3.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 2));
+        v4.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 3));
     }
 
-    public static void StoreInterleaved<T>(Vector256<T> v1, Vector256<T> v2, Vector256<T> v3, Vector256<T> v4, Vector256<T> v5, ref T memory)
+    public static unsafe void StoreInterleaved<T>(Vector256<T> v1, Vector256<T> v2, Vector256<T> v3, Vector256<T> v4, Vector256<T> v5, ref T memory)
+        where T : unmanaged
     {
-        v1.StoreUnsafe(ref Unsafe.Add(ref memory, 0));
-        v2.StoreUnsafe(ref Unsafe.Add(ref memory, 1));
-        v3.StoreUnsafe(ref Unsafe.Add(ref memory, 2));
-        v4.StoreUnsafe(ref Unsafe.Add(ref memory, 3));
-        v5.StoreUnsafe(ref Unsafe.Add(ref memory, 4));
+        v1.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 0));
+        v2.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 1));
+        v3.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 2));
+        v4.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 3));
+        v5.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 4));
     }
 
-    public static void StoreInterleaved<T>(Vector256<T> v1, Vector256<T> v2, Vector256<T> v3, Vector256<T> v4, Vector256<T> v5, Vector256<T> v6, ref T memory)
+    public static unsafe void StoreInterleaved<T>(Vector256<T> v1, Vector256<T> v2, Vector256<T> v3, Vector256<T> v4, Vector256<T> v5, Vector256<T> v6, ref T memory)
+        where T : unmanaged
     {
-        v1.StoreUnsafe(ref Unsafe.Add(ref memory, 0));
-        v2.StoreUnsafe(ref Unsafe.Add(ref memory, 1));
-        v3.StoreUnsafe(ref Unsafe.Add(ref memory, 2));
-        v4.StoreUnsafe(ref Unsafe.Add(ref memory, 3));
-        v5.StoreUnsafe(ref Unsafe.Add(ref memory, 4));
-        v6.StoreUnsafe(ref Unsafe.Add(ref memory, 5));
+        v1.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 0));
+        v2.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 1));
+        v3.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 2));
+        v4.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 3));
+        v5.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 4));
+        v6.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * 5));
     }
 
 }

src/ImageSharp/Formats/Jxl/Processing/JxlSimdUtils.StoreInterleaved.tt

@@ -33,10 +33,11 @@ internal static partial class JxlSimdUtils
             }
             string inlineParameters = string.Join(", ", vectorParameters) + ", ";
 #>
-    public static void StoreInterleaved<T>(<#= inlineParameters #>ref T memory)
+    public static unsafe void StoreInterleaved<T>(<#= inlineParameters #>ref T memory)
+        where T : unmanaged
     {
 <# for (int j = 0; j < i; j++) { #>
-        v<#= j + 1 #>.StoreUnsafe(ref Unsafe.Add(ref memory, <#= j #>));
+        v<#= j + 1 #>.StoreUnsafe(ref Unsafe.Add(ref memory, Vector<T>.Count * <#= j #>));
 <# } #>
     }
 

src/ImageSharp/Formats/Jxl/Processing/JxlSimdUtils.cs

@@ -15,10 +15,12 @@ namespace SixLabors.ImageSharp.Formats.Jxl.Processing;
 internal static partial class JxlSimdUtils
 {
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static Vector256<int> ConcatLowerLower(Vector256<int> a, Vector256<int> b) => Vector256.Create(a.GetLower(), b.GetLower());
+    public static Vector256<T> ConcatLowerLower<T>(Vector256<T> a, Vector256<T> b)
+        where T : unmanaged => Vector256.Create(a.GetLower(), b.GetLower());
 
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
-    private static Vector256<int> ConcatUpperUpper(Vector256<int> a, Vector256<int> b) => Vector256.Create(a.GetUpper(), b.GetUpper());
+    public static Vector256<T> ConcatUpperUpper<T>(Vector256<T> a, Vector256<T> b)
+        where T : unmanaged => Vector256.Create(a.GetUpper(), b.GetUpper());
 
     public static void Transpose8x8Block(Span<int> fromSpan, Span<int> toSpan, int stride)
     {
@@ -109,4 +111,330 @@ internal static partial class JxlSimdUtils
         Vector<float> vec = Vector.Log2(@base) * exponent;
         return vec * vec;
     }
+
+    /// <summary>
+    ///   <para>
+    ///     Fills the span so its first value is equal to <paramref name="start"/>
+    ///     and subsequent values increment by one. For example, with start=5,
+    ///     the span's values will be:
+    ///     <code>
+    ///       { start, start+1, start+2, start+3, start+4, ... to the end of the span }
+    ///     </code>
+    ///   </para>
+    ///   <para>
+    ///     or, more precisely:
+    ///     <code>
+    ///       { 5, 6, 7, 8, 9, 10, 11, ... to the end of the span }
+    ///     </code>
+    ///   </para>
+    ///   <seealso href="https://en.cppreference.com/cpp/algorithm/iota" />
+    /// </summary>
+    /// <param name="span">The span where the values are filled.</param>
+    /// <param name="start">Initial value.</param>
+    public static void Iota(Span<int> span, int start)
+    {
+        ref int spanRef = ref MemoryMarshal.GetReference(span);
+
+        // Using fixed-size vectors so we can construct an
+        // incrementMask more easily.
+        if (Vector512.IsHardwareAccelerated)
+        {
+            Vector512<int> incrementMask = Vector512.Create(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16);
+            Vector512<int> v = Vector512.Create(start) + (incrementMask - Vector512<int>.One);
+
+            if ((span.Length % Vector512<int>.Count) == 0)
+            {
+                // Aligned length
+                for (int i = 0; i < span.Length; i += Vector512<int>.Count)
+                {
+                    v.StoreUnsafe(ref Unsafe.Add(ref spanRef, i));
+                    v += incrementMask;
+                }
+            }
+            else
+            {
+                // We will need a scalar remainder
+                int vectorLength = span.Length - (span.Length % Vector512<int>.Count);
+
+                int i;
+                for (i = 0; i < vectorLength; i += Vector512<int>.Count)
+                {
+                    v.StoreUnsafe(ref Unsafe.Add(ref spanRef, i));
+                    v += incrementMask;
+                }
+
+                int val = v.ToScalar();
+                for (; i < span.Length; i++)
+                {
+                    Unsafe.Add(ref spanRef, i) = val++;
+                }
+            }
+        }
+        else if (Vector256.IsHardwareAccelerated)
+        {
+            Vector256<int> incrementMask = Vector256.Create(1, 2, 3, 4, 5, 6, 7, 8);
+            Vector256<int> v = Vector256.Create(start) + (incrementMask - Vector256<int>.One);
+
+            if ((span.Length % Vector256<int>.Count) == 0)
+            {
+                // Aligned length
+                for (int i = 0; i < span.Length; i += Vector256<int>.Count)
+                {
+                    v.StoreUnsafe(ref Unsafe.Add(ref spanRef, i));
+                    v += incrementMask;
+                }
+            }
+            else
+            {
+                // We will need a scalar remainder
+                int vectorLength = span.Length - (span.Length % Vector256<int>.Count);
+
+                int i;
+                for (i = 0; i < vectorLength; i += Vector256<int>.Count)
+                {
+                    v.StoreUnsafe(ref Unsafe.Add(ref spanRef, i));
+                    v += incrementMask;
+                }
+
+                int val = v.ToScalar();
+                for (; i < span.Length; i++)
+                {
+                    Unsafe.Add(ref spanRef, i) = val++;
+                }
+            }
+        }
+        else if (Vector128.IsHardwareAccelerated)
+        {
+            Vector128<int> incrementMask = Vector128.Create(1, 2, 3, 4);
+            Vector128<int> v = Vector128.Create(start) + (incrementMask - Vector128<int>.One);
+
+            if ((span.Length % Vector128<int>.Count) == 0)
+            {
+                // Aligned length
+                for (int i = 0; i < span.Length; i += Vector128<int>.Count)
+                {
+                    v.StoreUnsafe(ref Unsafe.Add(ref spanRef, i));
+                    v += incrementMask;
+                }
+            }
+            else
+            {
+                // We will need a scalar remainder
+                int vectorLength = span.Length - (span.Length % Vector128<int>.Count);
+
+                int i;
+                for (i = 0; i < vectorLength; i += Vector128<int>.Count)
+                {
+                    v.StoreUnsafe(ref Unsafe.Add(ref spanRef, i));
+                    v += incrementMask;
+                }
+
+                int val = v.ToScalar();
+                for (; i < span.Length; i++)
+                {
+                    Unsafe.Add(ref spanRef, i) = val++;
+                }
+            }
+        }
+        else if (Vector64.IsHardwareAccelerated)
+        {
+            Vector64<int> incrementMask = Vector64.Create(1, 2);
+            Vector64<int> v = Vector64.Create(start) + (incrementMask - Vector64<int>.One);
+
+            if ((span.Length % Vector64<int>.Count) == 0)
+            {
+                // Aligned length
+                for (int i = 0; i < span.Length; i += Vector64<int>.Count)
+                {
+                    v.StoreUnsafe(ref Unsafe.Add(ref spanRef, i));
+                    v += incrementMask;
+                }
+            }
+            else
+            {
+                // We will need a scalar remainder
+                int vectorLength = span.Length - (span.Length % Vector64<int>.Count);
+
+                int i;
+                for (i = 0; i < vectorLength; i += Vector64<int>.Count)
+                {
+                    v.StoreUnsafe(ref Unsafe.Add(ref spanRef, i));
+                    v += incrementMask;
+                }
+
+                int val = v.ToScalar();
+                for (; i < span.Length; i++)
+                {
+                    Unsafe.Add(ref spanRef, i) = val++;
+                }
+            }
+        }
+        else
+        {
+            // No SIMD
+            int value = start;
+            spanRef = value;
+            value++;
+            for (int i = 1; i < span.Length; i++)
+            {
+                Unsafe.Add(ref spanRef, i) = value++;
+            }
+        }
+    }
+
+    /// <summary>
+    ///   <para>
+    ///     Fills the span so its first value is equal to <paramref name="start"/>
+    ///     and subsequent values increment by one. For example, with start=5,
+    ///     the span's values will be:
+    ///     <code>
+    ///       { start, start+1, start+2, start+3, start+4, ... to the end of the span }
+    ///     </code>
+    ///   </para>
+    ///   <para>
+    ///     or, more precisely:
+    ///     <code>
+    ///       { 5, 6, 7, 8, 9, 10, 11, ... to the end of the span }
+    ///     </code>
+    ///   </para>
+    ///   <seealso href="https://en.cppreference.com/cpp/algorithm/iota" />
+    /// </summary>
+    /// <param name="span">The span where the values are filled.</param>
+    /// <param name="start">Initial value.</param>
+    public static void Iota<T>(Span<T> span, T start)
+        where T : unmanaged, INumber<T>
+    {
+        // Slightly slower than the int variant
+        ref T spanRef = ref MemoryMarshal.GetReference(span);
+
+        if (Vector<T>.IsSupported && Vector.IsHardwareAccelerated)
+        {
+            Vector<T> incrementMask = IotaMask<T>.IncrementMask;
+            Vector<T> v = Vector.Create(start) + (incrementMask - Vector<T>.One);
+
+            if ((span.Length % Vector<T>.Count) == 0)
+            {
+                // Aligned length
+                for (int i = 0; i < span.Length; i += Vector<T>.Count)
+                {
+                    v.StoreUnsafe(ref Unsafe.Add(ref spanRef, i));
+                    v += incrementMask;
+                }
+            }
+            else
+            {
+                // Remainder needed
+                int vectorLength = span.Length - (span.Length % Vector<T>.Count);
+
+                int i;
+                for (i = 0; i < vectorLength; i += Vector<T>.Count)
+                {
+                    v.StoreUnsafe(ref Unsafe.Add(ref spanRef, i));
+                    v += incrementMask;
+                }
+
+                T val = v.ToScalar();
+                for (; i < span.Length; i++)
+                {
+                    Unsafe.Add(ref spanRef, i) = val++;
+                }
+            }
+        }
+        else
+        {
+            // Scalar (slow)
+            T value = start;
+            spanRef = value;
+            value++;
+            for (int i = 1; i < span.Length; i++)
+            {
+                Unsafe.Add(ref spanRef, i) = value++;
+            }
+        }
+    }
+
+    public static Vector<T> Iota<T>(T start)
+        where T : unmanaged, INumber<T>
+        => IotaMask<T>.IncrementMask + Vector.Create(start);
+
+    /// <summary>
+    /// Vectorized floating-point error function (precise approximate).
+    /// </summary>
+    /// <param name="x">Vector to compute error of.</param>
+    /// <returns>Vector whose each item is an error (similar to std::erf).</returns>
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static Vector<float> FastErff(Vector<float> x)
+    {
+        Vector<float> zero = Vector<float>.Zero;
+        Vector<float> one = Vector<float>.One;
+
+        Vector<int> xle0 = Vector.LessThanOrEqual(x, zero);
+        Vector<float> absx = Vector.Abs(x);
+
+        Vector<float> denom1 = (absx * new Vector<float>(0.0777394369f)) + new Vector<float>(0.000205260015f);
+        Vector<float> denom2 = (denom1 * absx) + new Vector<float>(0.232120216f);
+        Vector<float> denom3 = (denom2 * absx) + new Vector<float>(0.277820801f);
+        Vector<float> denom4 = (denom3 * absx) + one;
+        Vector<float> denom5 = denom4 * denom4;
+        Vector<float> invDenom5 = one / denom5;
+        Vector<float> result = one - Vector.Multiply(invDenom5, invDenom5);
+
+        // Change sign if x <= 0.
+        return Vector.ConditionalSelect(xle0, -result, result);
+    }
+
+    /// <summary>
+    /// Scalar floating-point error function (precise approximate).
+    /// </summary>
+    /// <param name="x">Value to compute error of.</param>
+    /// <returns>A scalar error value (similar to std::erf).</returns>
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public static float FastErff(float x)
+    {
+        float zero = 0.0f;
+        float one = 1.0f;
+
+        bool xle0 = x <= zero;
+        float absx = MathF.Abs(x);
+
+        float denom1 = (absx * 0.0777394369f) + 0.000205260015f;
+        float denom2 = (denom1 * absx) + 0.232120216f;
+        float denom3 = (denom2 * absx) + 0.277820801f;
+        float denom4 = (denom3 * absx) + one;
+        float denom5 = denom4 * denom4;
+        float invDenom5 = one / denom5;
+        float result = one - (invDenom5 * invDenom5);
+
+        // Change sign if x <= 0.
+        return xle0 ? -result : result;
+    }
+
+    /// <summary>
+    /// Incrementing values to compute the Iota function.
+    /// </summary>
+    /// <remarks>
+    /// Creating a Vector&lt;T&gt; incrementing values would be
+    /// slow as Vector&lt;T&gt; is not a fixed-size vector, leaving
+    /// no other option but a slow loop. This class caches these
+    /// vectors for significantly better performance, though still
+    /// not as fast as an int variant.
+    /// </remarks>
+    /// <typeparam name="T">Type of the vector.</typeparam>
+    private static class IotaMask<T>
+        where T : unmanaged, INumber<T>
+    {
+        public static readonly Vector<T> IncrementMask;
+
+        static IotaMask()
+        {
+            Span<T> values = stackalloc T[Vector<T>.Count];
+
+            for (int i = 0; i < Vector<T>.Count; i++)
+            {
+                values[i] = T.CreateSaturating(i + 1);
+            }
+
+            IncrementMask = Vector.Create<T>(values);
+        }
+    }
 }

src/ImageSharp/Formats/Jxl/Processing/JxlTranspose.cs

@@ -1,6 +1,8 @@
 // Copyright (c) Six Labors.
 // Licensed under the Six Labors Split License.
 
+using System.Runtime.Intrinsics;
+using SixLabors.ImageSharp.Common.Helpers;
 using SixLabors.ImageSharp.Formats.Jxl.Processing.Dct;
 
 namespace SixLabors.ImageSharp.Formats.Jxl.Processing;
@@ -10,9 +12,29 @@ namespace SixLabors.ImageSharp.Formats.Jxl.Processing;
 /// </summary>
 internal static class JxlTranspose
 {
-    // TODO: SIMD
     public static void Transpose(int r, int c, JxlDctSource from, JxlDctOutput to)
     {
+        if (Vector256.IsHardwareAccelerated)
+        {
+            if (((r | c) & 7) == 0) // equivalent to (r % 8 == 0 && c % 8 == 0); micro-optimization, reduces one branch
+            {
+                // we can use SIMD
+                TransposeSimd256(r, c, from, to, r, c);
+                return;
+            }
+        }
+        else if (Vector128.IsHardwareAccelerated)
+        {
+            if (((r | c) & 3) == 0) // equivalent to (r % 4 == 0 && c % 4 == 0); micro-optimization, reduces one branch
+            {
+                // we can use SIMD
+                TransposeSimd128(r, c, from, to, r, c);
+                return;
+            }
+        }
+
+        // fallback: can't use SIMD (block size isn't aligned or
+        // there's no v128/v256 support)
         for (int n = 0; n < r; n++)
         {
             for (int m = 0; m < c; m++)
@@ -21,4 +43,93 @@ internal static class JxlTranspose
             }
         }
     }
+
+    private static void TransposeSimd256(int rowsOr0, int colsOr0, JxlDctSource from, JxlDctOutput to, int rowsP, int colsP)
+    {
+        int rows = rowsOr0 == 0 ? rowsP : rowsOr0;
+        int cols = colsOr0 == 0 ? colsP : colsOr0;
+
+        for (int n = 0; n < rows; n += 8)
+        {
+            for (int m = 0; m < cols; m += 8)
+            {
+                Vector256<float> i0 = from.LoadPart256(n, m);
+                Vector256<float> i1 = from.LoadPart256(n + 1, m);
+                Vector256<float> i2 = from.LoadPart256(n + 2, m);
+                Vector256<float> i3 = from.LoadPart256(n + 3, m);
+                Vector256<float> i4 = from.LoadPart256(n + 4, m);
+                Vector256<float> i5 = from.LoadPart256(n + 5, m);
+                Vector256<float> i6 = from.LoadPart256(n + 6, m);
+                Vector256<float> i7 = from.LoadPart256(n + 7, m);
+
+                Vector256<float> q0 = Vector256_.InterleaveLower(i0, i2);
+                Vector256<float> q1 = Vector256_.InterleaveLower(i1, i3);
+                Vector256<float> q2 = Vector256_.InterleaveUpper(i0, i2);
+                Vector256<float> q3 = Vector256_.InterleaveUpper(i1, i3);
+                Vector256<float> q4 = Vector256_.InterleaveLower(i4, i6);
+                Vector256<float> q5 = Vector256_.InterleaveLower(i5, i7);
+                Vector256<float> q6 = Vector256_.InterleaveUpper(i4, i6);
+                Vector256<float> q7 = Vector256_.InterleaveUpper(i5, i7);
+
+                Vector256<float> r0 = Vector256_.InterleaveLower(q0, q1);
+                Vector256<float> r1 = Vector256_.InterleaveUpper(q0, q1);
+                Vector256<float> r2 = Vector256_.InterleaveLower(q2, q3);
+                Vector256<float> r3 = Vector256_.InterleaveUpper(q2, q3);
+                Vector256<float> r4 = Vector256_.InterleaveLower(q4, q5);
+                Vector256<float> r5 = Vector256_.InterleaveUpper(q4, q5);
+                Vector256<float> r6 = Vector256_.InterleaveLower(q6, q7);
+                Vector256<float> r7 = Vector256_.InterleaveUpper(q6, q7);
+
+                i0 = JxlSimdUtils.ConcatLowerLower(r4, r0);
+                i1 = JxlSimdUtils.ConcatLowerLower(r5, r1);
+                i2 = JxlSimdUtils.ConcatLowerLower(r6, r2);
+                i3 = JxlSimdUtils.ConcatLowerLower(r7, r3);
+                i4 = JxlSimdUtils.ConcatUpperUpper(r4, r0);
+                i5 = JxlSimdUtils.ConcatUpperUpper(r5, r1);
+                i6 = JxlSimdUtils.ConcatUpperUpper(r6, r2);
+                i7 = JxlSimdUtils.ConcatUpperUpper(r7, r3);
+
+                to.StorePart256(i0, m, n);
+                to.StorePart256(i1, m + 1, n);
+                to.StorePart256(i2, m + 2, n);
+                to.StorePart256(i3, m + 3, n);
+                to.StorePart256(i4, m + 4, n);
+                to.StorePart256(i5, m + 5, n);
+                to.StorePart256(i6, m + 6, n);
+                to.StorePart256(i7, m + 7, n);
+            }
+        }
+    }
+
+    private static void TransposeSimd128(int rowsOr0, int colsOr0, JxlDctSource from, JxlDctOutput to, int rowsP, int colsP)
+    {
+        int rows = rowsOr0 == 0 ? rowsP : rowsOr0;
+        int cols = colsOr0 == 0 ? colsP : colsOr0;
+
+        for (int n = 0; n < rows; n += 4)
+        {
+            for (int m = 0; m < cols; m += 4)
+            {
+                Vector128<float> p0 = from.LoadPart128(n, m);
+                Vector128<float> p1 = from.LoadPart128(n + 1, m);
+                Vector128<float> p2 = from.LoadPart128(n + 2, m);
+                Vector128<float> p3 = from.LoadPart128(n + 3, m);
+
+                Vector128<float> q0 = Vector128_.InterleaveLower(p0, p2);
+                Vector128<float> q1 = Vector128_.InterleaveLower(p1, p3);
+                Vector128<float> q2 = Vector128_.InterleaveUpper(p0, p2);
+                Vector128<float> q3 = Vector128_.InterleaveUpper(p1, p3);
+
+                Vector128<float> r0 = Vector128_.InterleaveLower(q0, q1);
+                Vector128<float> r1 = Vector128_.InterleaveUpper(q0, q1);
+                Vector128<float> r2 = Vector128_.InterleaveLower(q2, q3);
+                Vector128<float> r3 = Vector128_.InterleaveUpper(q2, q3);
+
+                to.StorePart128(r0, m, n);
+                to.StorePart128(r1, m + 1, n);
+                to.StorePart128(r2, m + 2, n);
+                to.StorePart128(r3, m + 3, n);
+            }
+        }
+    }
 }

src/ImageSharp/Formats/Jxl/Processing/Primitives/JxlLehmerCode.cs

@@ -10,16 +10,16 @@ internal static class JxlLehmerCode
     [MethodImpl(MethodImplOptions.AggressiveInlining)]
     private static int ValueOfLowest1Bit(int n) => n & -n;
 
-    public static bool ComputeLehmerCode(ReadOnlySpan<int> permutation, Span<uint> temp, int n, Span<uint> code)
+    public static bool ComputeLehmerCode(ReadOnlySpan<uint> permutation, Span<uint> temp, int n, Span<uint> code)
     {
         temp[(n + 1)..].Clear();
 
         for (int idx = 0; idx < n; idx++)
         {
-            int s = permutation[idx];
+            uint s = permutation[idx];
 
             uint penalty = 0u;
-            uint i = (uint)s + 1u;
+            uint i = s + 1u;
 
             while (i != 0u)
             {
@@ -32,8 +32,8 @@ internal static class JxlLehmerCode
                 return false;
             }
 
-            code[idx] = (uint)s - penalty;
-            i = (uint)s + 1u;
+            code[idx] = s - penalty;
+            i = s + 1u;
 
             while (i < n + 1u)
             {
@@ -45,7 +45,7 @@ internal static class JxlLehmerCode
         return true;
     }
 
-    public static bool DecodeLehmerCode(ReadOnlySpan<uint> code, Span<uint> temp, int n, Span<int> permutation)
+    public static bool DecodeLehmerCode(ReadOnlySpan<uint> code, Span<uint> temp, int n, Span<uint> permutation)
     {
         if (n == 0)
         {
@@ -91,7 +91,7 @@ internal static class JxlLehmerCode
                 }
             }
 
-            permutation[i] = next;
+            permutation[i] = unchecked((uint)next);
 
             next++;
             while (next <= paddedN)

src/ImageSharp/Formats/Jxl/Processing/Primitives/JxlXorShift.cs

@@ -12,7 +12,7 @@ internal sealed class JxlXorShift
     private readonly ulong[] s0 = new ulong[8];
     private readonly ulong[] s1 = new ulong[8];
 
-    public void XorShift128Plus(ulong seed)
+    public JxlXorShift(ulong seed)
     {
         this.s0[0] = SplitMix64(seed + 0x9E3779B97F4A7C15L);
         this.s1[0] = SplitMix64(this.s0[0]);

src/ImageSharp/Formats/Jxl/Processing/Quantization/JxlQuantizerSimd.cs

@@ -0,0 +1,28 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using System.Numerics;
+
+namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Quantization;
+
+/// <summary>
+/// SIMD utilities used by the quantizer.
+/// </summary>
+internal static class JxlQuantizerSimd
+{
+    public static Vector<float> AdjustQuantBias(int c, Vector<int> quantI, Span<float> biases)
+    {
+        Vector<float> quant = quantI.As<int, float>();
+        Vector<float> constSign = Vector.Create(int.MinValue).As<int, float>();
+        Vector<float> sign = quant & constSign;
+        Vector<float> absoluteQuant = Vector.AndNot(constSign, quant);
+
+        Vector<int> is01 = Vector.LessThan(absoluteQuant, Vector.Create(1.125f));
+        Vector<int> not0 = Vector.GreaterThan(absoluteQuant, Vector<float>.One);
+
+        Vector<float> oneBias = Vector.ConditionalSelect(not0, Vector.Create(biases[c]) ^ sign, Vector<float>.Zero);
+        Vector<float> bias = -(Vector.Create(biases[3]) * quant.ReciprocalEstimate()) + quant;
+
+        return Vector.ConditionalSelect(is01, oneBias, bias);
+    }
+}

src/ImageSharp/Formats/Jxl/Processing/Splines/Dct32.cs

@@ -0,0 +1,15 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using System.Runtime.CompilerServices;
+
+namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Splines;
+
+/// <summary>
+/// Storage for 32 DCT coefficients (floating-point).
+/// </summary>
+[InlineArray(32)]
+internal struct Dct32
+{
+    private float first;
+}

src/ImageSharp/Formats/Jxl/Processing/Splines/JxlQuantizedSpline.cs

@@ -224,7 +224,7 @@ internal sealed class JxlQuantizedSpline : IDisposable
 
         for (int i = 0; i < 32; i++)
         {
-            float inverseDctFactor = (i == 0) ? Sqrt05 : 1.0f;
+            float inverseDctFactor = (i == 0) ? JxlDctScales.Sqrt05 : 1.0f;
             result.SigmaDct[i] = this.SigmaDct[i] * inverseDctFactor * ChannelWeight[3] * inverseQuant;
             float weightF = MathF.Ceiling(inverseQuant * MathF.Abs(this.SigmaDct[i]));
             long weight = (long)Math.Min(weightLimit, Math.Max(1.0f, weightF));
@@ -240,7 +240,7 @@ internal sealed class JxlQuantizedSpline : IDisposable
         return true;
     }
 
-    public bool Decode(
+    public bool TryDecode(
         Configuration configuration,
         Span<byte> contextMap,
         JxlAnsSymbolReader decoder,

src/ImageSharp/Formats/Jxl/Processing/Splines/JxlSpline.cs

@@ -11,9 +11,9 @@ internal sealed class JxlSpline : IDisposable
 
     public Memory<PointF> ControlPoints { get; private set; }
 
-    public JxlDct32[] ColorDct { get; set; } = [];
+    public Dct32[] ColorDct { get; set; } = [];
 
-    public JxlDct32 SigmaDct { get; set; }
+    public Dct32 SigmaDct { get; set; }
 
     public void ClearControlPoints()
     {

src/ImageSharp/Formats/Jxl/Processing/Splines/JxlSplineUtils.cs

@@ -0,0 +1,394 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using System.Numerics;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
+using System.Runtime.Intrinsics;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Dct;
+
+namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Splines;
+
+internal static class JxlSplineUtils
+{
+    public const float DesiredRenderingDistance = 1f;
+
+    private static ReadOnlySpan<float> ContinuousIDCTMultipliers =>
+    [
+        MathF.PI / 32 * 0,  MathF.PI / 32 * 1,  MathF.PI / 32 * 2,  MathF.PI / 32 * 3,  MathF.PI / 32 * 4,
+        MathF.PI / 32 * 5,  MathF.PI / 32 * 6,  MathF.PI / 32 * 7,  MathF.PI / 32 * 8,  MathF.PI / 32 * 9,
+        MathF.PI / 32 * 10, MathF.PI / 32 * 11, MathF.PI / 32 * 12, MathF.PI / 32 * 13, MathF.PI / 32 * 14,
+        MathF.PI / 32 * 15, MathF.PI / 32 * 16, MathF.PI / 32 * 17, MathF.PI / 32 * 18, MathF.PI / 32 * 19,
+        MathF.PI / 32 * 20, MathF.PI / 32 * 21, MathF.PI / 32 * 22, MathF.PI / 32 * 23, MathF.PI / 32 * 24,
+        MathF.PI / 32 * 25, MathF.PI / 32 * 26, MathF.PI / 32 * 27, MathF.PI / 32 * 28, MathF.PI / 32 * 29,
+        MathF.PI / 32 * 30, MathF.PI / 32 * 31,
+    ];
+
+    public static float ContinuousInverseDCT(in Dct32 dct, float t)
+    {
+        ref float multipliers = ref MemoryMarshal.GetReference(ContinuousIDCTMultipliers);
+        ReadOnlySpan<float> dctData = dct;
+        ref float dctDataRef = ref MemoryMarshal.GetReference(dctData);
+
+        if (Vector<float>.Count <= 32 && Vector.IsHardwareAccelerated)
+        {
+            Vector<float> result = Vector<float>.Zero;
+            Vector<float> tandhalf = Vector.Create(t + 0.5f);
+
+            for (int i = 0; i < 32; i += Vector<float>.Count)
+            {
+                Vector<float> cosArg = Vector.LoadUnsafe(ref Unsafe.Add(ref multipliers, i)) * tandhalf;
+                Vector<float> cos = Vector.Cos(cosArg);
+                Vector<float> localRes = Vector.LoadUnsafe(ref Unsafe.Add(ref dctDataRef, i)) * cos;
+                result = (Vector.Create(JxlDctScales.Sqrt2) * localRes) + result;
+            }
+
+            return Vector.Sum(result);
+        }
+
+        // Might have SIMD support but Vector<T> > 32 (e.g. on some CPUs),
+        // so let's try different fixed-size vectors first.
+        else if (Vector512.IsHardwareAccelerated)
+        {
+            Vector512<float> result = Vector512<float>.Zero;
+            Vector512<float> tandhalf = Vector512.Create(t + 0.5f);
+
+            for (int i = 0; i < 32; i += Vector512<float>.Count)
+            {
+                Vector512<float> cosArg = Vector512.LoadUnsafe(ref Unsafe.Add(ref multipliers, i)) * tandhalf;
+                Vector512<float> cos = Vector512.Cos(cosArg);
+                Vector512<float> localRes = Vector512.LoadUnsafe(ref Unsafe.Add(ref dctDataRef, i)) * cos;
+                result = (Vector512.Create(JxlDctScales.Sqrt2) * localRes) + result;
+            }
+
+            return Vector512.Sum(result);
+        }
+        else if (Vector256.IsHardwareAccelerated)
+        {
+            Vector256<float> result = Vector256<float>.Zero;
+            Vector256<float> tandhalf = Vector256.Create(t + 0.5f);
+
+            for (int i = 0; i < 32; i += Vector256<float>.Count)
+            {
+                Vector256<float> cosArg = Vector256.LoadUnsafe(ref Unsafe.Add(ref multipliers, i)) * tandhalf;
+                Vector256<float> cos = Vector256.Cos(cosArg);
+                Vector256<float> localRes = Vector256.LoadUnsafe(ref Unsafe.Add(ref dctDataRef, i)) * cos;
+                result = (Vector256.Create(JxlDctScales.Sqrt2) * localRes) + result;
+            }
+
+            return Vector256.Sum(result);
+        }
+        else if (Vector128.IsHardwareAccelerated)
+        {
+            Vector128<float> result = Vector128<float>.Zero;
+            Vector128<float> tandhalf = Vector128.Create(t + 0.5f);
+
+            for (int i = 0; i < 32; i += Vector128<float>.Count)
+            {
+                Vector128<float> cosArg = Vector128.LoadUnsafe(ref Unsafe.Add(ref multipliers, i)) * tandhalf;
+                Vector128<float> cos = Vector128.Cos(cosArg);
+                Vector128<float> localRes = Vector128.LoadUnsafe(ref Unsafe.Add(ref dctDataRef, i)) * cos;
+                result = (Vector128.Create(JxlDctScales.Sqrt2) * localRes) + result;
+            }
+
+            return Vector128.Sum(result);
+        }
+        else if (Vector64.IsHardwareAccelerated)
+        {
+            Vector64<float> result = Vector64<float>.Zero;
+            Vector64<float> tandhalf = Vector64.Create(t + 0.5f);
+
+            for (int i = 0; i < 32; i += Vector64<float>.Count)
+            {
+                Vector64<float> cosArg = Vector64.LoadUnsafe(ref Unsafe.Add(ref multipliers, i)) * tandhalf;
+                Vector64<float> cos = Vector64.Cos(cosArg);
+                Vector64<float> localRes = Vector64.LoadUnsafe(ref Unsafe.Add(ref dctDataRef, i)) * cos;
+                result = (Vector64.Create(JxlDctScales.Sqrt2) * localRes) + result;
+            }
+
+            return Vector64.Sum(result);
+        }
+        else
+        {
+            // Scalar fallback.
+            float result = 0f;
+            float tandhalf = t + 0.5f;
+
+            for (int i = 0; i < 32; i++)
+            {
+                float cosArg = Unsafe.Add(ref multipliers, i) * tandhalf;
+                float cos = MathF.Cos(cosArg);
+                float localRes = Unsafe.Add(ref dctDataRef, i) * cos;
+                result = (JxlDctScales.Sqrt2 * localRes) + result;
+            }
+
+            return result;
+        }
+    }
+
+    // SIMD version
+    private static void DrawSegmentPacked(ref JxlSplineSegment segment, bool add, int y, int x, int x0, InlineArray3<Memory<float>> rows)
+    {
+        Vector<float> inverseSigma = Vector.Create(segment.InverseSigma);
+        Vector<float> half = Vector.Create(0.5f);
+        Vector<float> oneOver2s2 = Vector.Create(0.353553391f);
+        Vector<float> sigmaOver4TimesIntensity = Vector.Create(segment.SigmaOver4TimesIntensity);
+
+        Vector<float> dx = JxlSimdUtils.Iota(x + x0).As<int, float>() - Vector.Create(segment.Center.X);
+        Vector<float> dy = Vector.Create(y - segment.Center.Y);
+
+        Vector<float> sqd = (dx * dx) + (dy * dy);
+        Vector<float> distance = Vector.SquareRoot(sqd);
+
+        Vector<float> oneDimensionalFactor =
+            JxlSimdUtils.FastErff(((distance * half) + oneOver2s2) * inverseSigma)
+            - JxlSimdUtils.FastErff(((distance * half) - oneOver2s2) * inverseSigma);
+
+        Vector<float> localIntensity = sigmaOver4TimesIntensity * (oneDimensionalFactor * oneDimensionalFactor);
+
+        for (int c = 0; c < 3; c++)
+        {
+            Span<float> currRow = rows[c].Span;
+            ref float currRowRef = ref MemoryMarshal.GetReference(currRow);
+
+            // TODO: move the add branch outside the loop and duplicate the
+            // loops twice? this removes the branch
+            Vector<float> cm = Vector.Create(add ? segment.Color[c] : -segment.Color[c]);
+
+            Vector<float> @in = Vector.LoadUnsafe(ref Unsafe.Add(ref currRowRef, x));
+            ((cm * localIntensity) + @in).StoreUnsafe(ref Unsafe.Add(ref currRowRef, x));
+        }
+    }
+
+    // Scalar version (for remaining items left to process)
+    private static void DrawSegmentScalar(ref JxlSplineSegment segment, bool add, int y, int x, int x0, InlineArray3<Memory<float>> rows)
+    {
+        float inverseSigma = segment.InverseSigma;
+        float half = 0.5f;
+        float oneOver2s2 = 0.353553391f;
+        float sigmaOver4TimesIntensity = segment.SigmaOver4TimesIntensity;
+
+        float dx = (x + x0) - segment.Center.X;
+        float dy = y - segment.Center.Y;
+
+        float sqd = (dx * dx) + (dy * dy);
+        float distance = MathF.Sqrt(sqd);
+
+        float oneDimensionalFactor =
+            JxlSimdUtils.FastErff(((distance * half) + oneOver2s2) * inverseSigma)
+            - JxlSimdUtils.FastErff(((distance * half) - oneOver2s2) * inverseSigma);
+
+        float localIntensity = sigmaOver4TimesIntensity * (oneDimensionalFactor * oneDimensionalFactor);
+
+        for (int c = 0; c < 3; c++)
+        {
+            Span<float> currRow = rows[c].Span;
+            float cm = add ? segment.Color[c] : -segment.Color[c];
+            currRow[x] = (cm * localIntensity) + currRow[x];
+        }
+    }
+
+    public static void DrawSegment(ref JxlSplineSegment segment, bool add, int y, int x0, int x1, InlineArray3<Memory<float>> rows)
+    {
+        int start = (int)MathF.Round(segment.Center.X - segment.MaximumDistance, MidpointRounding.AwayFromZero);
+        int end = (int)MathF.Round(segment.Center.X + segment.MaximumDistance, MidpointRounding.AwayFromZero);
+
+        if (end < x0 || start >= x1)
+        {
+            return;  // span does not intersect scan
+        }
+
+        int spanX0 = Math.Max(x0, start) - x0;
+        int spanX1 = Math.Min(x1, end + 1) - x0;
+
+        int x = spanX0;
+        for (; x + Vector<float>.Count <= spanX1; x += Vector<float>.Count)
+        {
+            DrawSegmentPacked(ref segment, add, y, x, x0, rows);
+        }
+
+        for (; x < spanX1; ++x)
+        {
+            DrawSegmentScalar(ref segment, add, y, x, x0, rows);
+        }
+    }
+
+    public static void ComputeSegments(int imageYSize, PointF center, float intensity, InlineArray3<float> color, float sigma, List<JxlSplineSegment> segments, List<JxlSplineSegmentSpan> segmentSpans)
+    {
+        if (!(float.IsFinite(sigma) && sigma != 0.0f && float.IsFinite(1.0f / sigma) && float.IsFinite(intensity)))
+        {
+            return;
+        }
+
+        // This is about 30% faster, but for higher precision
+        // one can change this to 5 instead.
+        const float distanceExp = 3f;
+
+        float maxColor = MathF.Max(0.01f, MathF.Abs(color[0] * intensity));
+        maxColor = MathF.Max(maxColor, MathF.Abs(color[1] * intensity));
+        maxColor = MathF.Max(maxColor, MathF.Abs(color[2] * intensity));
+
+        float maximumDistance = MathF.Sqrt(-2.0f * sigma * sigma * ((MathF.Log(0.1f) * distanceExp) - MathF.Log(maxColor)));
+
+        int y0 = (int)MathF.Round(center.Y - maximumDistance, MidpointRounding.AwayFromZero);
+        y0 = Math.Max(y0, 0);
+
+        int y1 = (int)MathF.Round(center.Y + maximumDistance, MidpointRounding.AwayFromZero) + 1;
+        y1 = Math.Min(y1, imageYSize);
+
+        if (y1 <= y0)
+        {
+            return;
+        }
+
+        JxlSplineSegment segment = new()
+        {
+            Center = center,
+            InverseSigma = 1.0f / sigma,
+            SigmaOver4TimesIntensity = 0.25f * sigma * intensity,
+            MaximumDistance = maximumDistance,
+            Color = color
+        };
+
+        segments.Add(segment);
+        segmentSpans.Add(new JxlSplineSegmentSpan(y0, y1));
+    }
+
+    public static void DrawSegments(Memory<float> rowX, Memory<float> rowY, Memory<float> rowB, int y, int x0, int x1, bool add, Span<JxlSplineSegment> segments, Span<int> segmentIndices, Span<int> segmentYStart)
+    {
+        InlineArray3<Memory<float>> rows = default;
+        rows[0] = rowX;
+        rows[1] = rowY;
+        rows[2] = rowB;
+
+        for (int i = segmentYStart[y]; i < segmentYStart[y + 1]; i++)
+        {
+            DrawSegment(ref segments[segmentIndices[i]], add, y, x0, x1, rows);
+        }
+    }
+
+    public static void SegmentsFromPoints(int imageYSize, JxlSpline spline, List<(PointF Point, float Multiplier)> pointsToDraw, float arcLength, List<JxlSplineSegment> segments, List<JxlSplineSegmentSpan> segmentsSpans)
+    {
+        float inverseArcLength = 1.0f / arcLength;
+        int k = 0;
+
+        foreach ((PointF point, float multiplier) in pointsToDraw)
+        {
+            float progressAlongArc = MathF.Min(1.0f, (k++ * DesiredRenderingDistance) * inverseArcLength);
+
+            InlineArray3<float> color = default;
+            color[0] = ContinuousInverseDCT(spline.ColorDct[0], (32 - 1) * progressAlongArc);
+            color[1] = ContinuousInverseDCT(spline.ColorDct[1], (32 - 1) * progressAlongArc);
+            color[2] = ContinuousInverseDCT(spline.ColorDct[2], (32 - 1) * progressAlongArc);
+
+            float sigma = ContinuousInverseDCT(spline.SigmaDct, (32 - 1) * progressAlongArc);
+            ComputeSegments(imageYSize, point, multiplier, color, sigma, segments, segmentsSpans);
+        }
+    }
+
+    public static void DrawCentripetalCatmullRomSpline(Span<PointF> points, List<PointF> result)
+    {
+        if (points.Length == 0)
+        {
+            return;
+        }
+
+        if (points.Length == 1)
+        {
+            result.Add(points[0]);
+            return;
+        }
+
+        List<PointF> pointsCopy = [];
+        for (int i = 0; i < points.Length; i++)
+        {
+            pointsCopy.Add(points[i]);
+        }
+
+        const int numPoints = 16;
+        pointsCopy.Insert(0, pointsCopy[0] + (pointsCopy[0] - pointsCopy[1]));
+        pointsCopy.Add(pointsCopy[^1] + (pointsCopy[^1] - pointsCopy[^2]));
+
+        for (int start = 0; start < pointsCopy.Count - 3; start++)
+        {
+            Span<PointF> p = CollectionsMarshal.AsSpan(pointsCopy)[start..];
+            result.Add(p[1]);
+
+            InlineArray3<float> d = default;
+            InlineArray4<float> t = default;
+
+            for (int k = 0; k < 3; ++k)
+            {
+                d[k] = MathF.Sqrt(JxlMath.Hypot(p[k + 1].X - p[k].X, p[k + 1].Y - p[k].Y));
+                t[k + 1] = t[k] + d[k];
+            }
+
+            for (int i = 1; i < numPoints; ++i)
+            {
+                float tt = d[0] + (((float)i / numPoints) * d[1]);
+                InlineArray3<PointF> a = default;
+
+                for (int k = 0; k < 3; ++k)
+                {
+                    a[k] = p[k] + (((tt - t[k]) / d[k]) * (p[k + 1] - p[k]));
+                }
+
+                InlineArray3<PointF> b = default;
+
+                for (int k = 0; k < 2; ++k)
+                {
+                    b[k] = a[k] + (((tt - t[k]) / (d[k] + d[k + 1])) * (a[k + 1] - a[k]));
+                }
+
+                result.Add(b[0] + (((tt - t[1]) / d[1]) * (b[1] - b[0])));
+            }
+        }
+
+        result.Add(pointsCopy[^2]);
+    }
+
+    public static void ForEachEquallySpacedPoint(Span<PointF> points, Action<PointF, float> functor)
+    {
+        PointF current = points[0];
+        functor(current, DesiredRenderingDistance);
+
+        ref PointF next = ref points[0];
+        ref PointF end = ref points[^1]; // last
+
+        while (!Unsafe.AreSame(ref next, ref end))
+        {
+            ref PointF previous = ref current;
+            float arcLengthFromPrevious = 0f;
+
+            while (true)
+            {
+                if (next == end)
+                {
+                    functor(previous, arcLengthFromPrevious);
+                    return;
+                }
+
+                float arcLengthToNext = MathF.Sqrt(SquaredNorm(next - previous));
+
+                if (arcLengthFromPrevious + arcLengthToNext >= DesiredRenderingDistance)
+                {
+                    current = previous + (((DesiredRenderingDistance - arcLengthFromPrevious) / arcLengthToNext) * (next - previous));
+                    functor(current, DesiredRenderingDistance);
+                    break;
+                }
+
+                arcLengthFromPrevious += arcLengthToNext;
+                previous = ref next;
+                next = ref Unsafe.Add(ref next, 1);
+            }
+        }
+    }
+
+    private static float SquaredNorm(PointF pointF)
+    {
+        float x = pointF.X;
+        float y = pointF.Y;
+
+        return (x * x) + (y * y);
+    }
+}

src/ImageSharp/Formats/Jxl/Processing/Splines/JxlSplines.cs

@@ -0,0 +1,289 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using System.Buffers;
+using System.Runtime.InteropServices;
+using SixLabors.ImageSharp.Formats.Jxl.Memory.ImageTypes;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Decoder;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Primitives;
+
+namespace SixLabors.ImageSharp.Formats.Jxl.Processing.Splines;
+
+internal sealed class JxlSplines
+{
+    private List<JxlQuantizedSpline> splinesStorage = [];
+    private readonly List<PointF> startingPointsStorage = [];
+    private JxlSplineDataView data = new();
+    private readonly List<JxlSplineSegment> segments = [];
+    private IMemoryOwner<byte> segmentIndices = EmptyMemoryOwner.Instance;
+    private IMemoryOwner<byte> segmentYStart = EmptyMemoryOwner.Instance;
+
+    public bool HasAny => this.data.HasAny;
+
+    public Span<JxlQuantizedSpline> QuantizedSplines => CollectionsMarshal.AsSpan(this.data.Splines);
+
+    public Span<PointF> StartingPoints => CollectionsMarshal.AsSpan(this.data.StartingPoints);
+
+    public int QuantizationAdjustment { get; private set; }
+
+    public void SetData(JxlSplineDataView data)
+    {
+        this.Clear();
+        this.data = data;
+    }
+
+    public void Clear()
+    {
+        this.QuantizationAdjustment = 0;
+        this.splinesStorage.Clear();
+        this.startingPointsStorage.Clear();
+        this.data = new();
+        this.segments.Clear();
+        this.segmentIndices = EmptyMemoryOwner.Instance;
+        this.segmentYStart = EmptyMemoryOwner.Instance;
+    }
+
+    public void Decode(Configuration configuration, JxlBitReader reader, int numPixels)
+    {
+        List<byte> contextMap = [];
+
+        JxlAnsReader.DecodeHistograms(reader, NumSplineContexts, out JxlAnsCode code, contextMap);
+        JxlAnsSymbolReader decoder = new(code, reader);
+
+        int numSplines = decoder.ReadHybridUint(NumSplinesContext, reader, contextMap);
+        int maxControlPoints = Math.Min(MaxNumControlPoints, numPixels / MaxNumControlPointsPerPixelRatio);
+
+        if (numSplines > maxControlPoints || numSplines + 1 > maxControlPoints)
+        {
+            throw new InvalidOperationException("Too many splines: " + numSplines);
+        }
+
+        numSplines++;
+
+        DecodeAllStartingPoints(this.startingPointsStorage, reader, decoder, contextMap, numSplines);
+
+        this.QuantizationAdjustment = JxlPackSigned.UnpackSigned(decoder.ReadHybridUint(QuantizationAdjustmentContext, reader, contextMap));
+        this.splinesStorage = new List<JxlQuantizedSpline>(numSplines);
+
+        int numControlPoints = numSplines;
+
+        for (int i = 0; i < numSplines; ++i)
+        {
+            JxlQuantizedSpline spline = new();
+            if (!spline.TryDecode(
+                configuration,
+                CollectionsMarshal.AsSpan(contextMap),
+                decoder,
+                reader,
+                maxControlPoints,
+                ref numControlPoints))
+            {
+                throw new InvalidOperationException("Could not decode quantized spline. Index of the quantized spline: " + i);
+            }
+
+            this.splinesStorage.Add(spline);
+        }
+
+        if (!decoder.CheckAnsFinalState())
+        {
+            throw new InvalidOperationException("Not ANS final state");
+        }
+
+        this.data = new JxlSplineDataView()
+        {
+            Splines = this.splinesStorage,
+            StartingPoints = this.startingPointsStorage
+        };
+
+        if (!this.HasAny)
+        {
+            throw new InvalidOperationException("Decoded splines but got none");
+        }
+    }
+
+    public void InitializeDrawCache(Configuration configuration, int imageXSize, int imageYSize, JxlColorCorrelation colorCorrelation)
+    {
+        this.segments.Clear();
+        this.segmentIndices = EmptyMemoryOwner.Instance;
+        this.segmentYStart = EmptyMemoryOwner.Instance;
+
+        List<JxlSplineSegmentSpan> segmentsSpans = [];
+        List<PointF> intermediatePoints = [];
+        List<JxlSpline> splines = [];
+        long totalEstimatedAreaReached = 0;
+
+        for (int i = 0; i < this.data.Splines.Count; i++)
+        {
+            JxlSpline spline = new();
+
+            if (!this.data.Splines[i].Dequantize(
+                configuration,
+                this.data.StartingPoints[i],
+                this.QuantizationAdjustment,
+                colorCorrelation.YToXRatio(0),
+                colorCorrelation.YToBRatio(0),
+                imageXSize * imageYSize,
+                ref totalEstimatedAreaReached,
+                spline))
+            {
+                throw new InvalidOperationException("Could not dequantize a quantized spline");
+            }
+
+            if (AdjacentFind(spline.ControlPoints.Span) != spline.ControlPoints.Length - 1)
+            {
+                throw new InvalidOperationException("Identical successive control points in spline " + i);
+            }
+
+            splines.Add(spline);
+        }
+
+#if JPEG_XL_THROW_ON_LARGE_SPLINE_AREA
+        if (totalEstimatedAreaReached > Math.Min((8 * imageXSize * imageYSize) + (1 << 25), 1 << 30))
+        {
+            throw new InvalidOperationException("Total spline area is too large");
+        }
+#endif
+
+        foreach (JxlSpline spline in splines)
+        {
+            List<(PointF Point, float Multiplier)> pointsToDraw = [];
+
+            void AddPoint(PointF point, float multiplier) => pointsToDraw.Add((point, multiplier));
+
+            intermediatePoints.Clear();
+
+            JxlSplineUtils.DrawCentripetalCatmullRomSpline(spline.ControlPoints.Span, intermediatePoints);
+            JxlSplineUtils.ForEachEquallySpacedPoint(CollectionsMarshal.AsSpan(intermediatePoints), AddPoint);
+
+            float arcLength = ((pointsToDraw.Count - 2) * JxlSplineUtils.DesiredRenderingDistance) + pointsToDraw[^1].Multiplier;
+            if (arcLength <= 0f)
+            {
+                // This spline wouldn't have any effect.
+                continue;
+            }
+
+            JxlSplineUtils.SegmentsFromPoints(imageYSize, spline, pointsToDraw, arcLength, this.segments, segmentsSpans);
+        }
+
+        int segmentYStartNumBytes = (imageYSize + 2) * 4;
+        this.segmentYStart = configuration.MemoryAllocator.Allocate<byte>(segmentYStartNumBytes);
+
+        Span<byte> segmentYStart = this.segmentYStart.Memory.Span;
+        segmentYStart.Clear();
+
+        Span<byte> population = segmentYStart[1..];
+
+        foreach (JxlSplineSegmentSpan segmentSpan in segmentsSpans)
+        {
+            population[segmentSpan.StartInclusive]++;
+            population[segmentSpan.EndInclusive]--;
+        }
+
+        int total = 0;
+        int coverage = 0;
+
+        for (int y = 0; y < imageYSize; y++)
+        {
+            if (population[y] < 0)
+            {
+                if (coverage < -population[y])
+                {
+                    throw new InvalidOperationException("Coverage is invalid");
+                }
+            }
+
+            coverage += population[y];
+            population[y] = (byte)total;
+            total += coverage;
+        }
+
+        this.segmentIndices = configuration.MemoryAllocator.Allocate<byte>(total * 4);
+        Span<int> segmentIndices = MemoryMarshal.Cast<byte, int>(this.segmentIndices.Memory.Span);
+
+        for (int i = 0; i < this.segments.Count; i++)
+        {
+            JxlSplineSegmentSpan segmentSpan = segmentsSpans[i];
+
+            for (int y = segmentSpan.StartInclusive; y < segmentSpan.EndInclusive; y++)
+            {
+                segmentIndices[population[y]++] = i;
+            }
+        }
+    }
+
+    private static int AdjacentFind<T>(Span<T> span)
+        where T : IEquatable<T>
+    {
+        for (int i = 0; i < span.Length - 1; i++)
+        {
+            if (span[i].Equals(span[i + 1]))
+            {
+                return i;
+            }
+        }
+
+        return span.Length;
+    }
+
+    public void AddTo(JxlImage3F opsin, Rectangle opsinRect) => this.Apply(add: true, opsin, opsinRect);
+
+    public void AddToRow(Memory<float> rowX, Memory<float> rowY, Memory<float> rowB, int y, int x0, int x1)
+        => this.ApplyToRow(add: true, rowX, rowY, rowB, y, x0, x1);
+
+    public void SubtractFrom(JxlImage3F opsin) => this.Apply(add: false, opsin, opsin.GetRectangle());
+
+    private void ApplyToRow(bool add, Memory<float> rowX, Memory<float> rowY, Memory<float> rowB, int y, int x0, int x1)
+    {
+        if (this.segments.Count == 0)
+        {
+            return;
+        }
+
+        JxlSplineUtils.DrawSegments(
+            rowX,
+            rowY,
+            rowB,
+            y,
+            x0,
+            x1,
+            add,
+            CollectionsMarshal.AsSpan(this.segments),
+            MemoryMarshal.Cast<byte, int>(this.segmentIndices.Memory.Span),
+            MemoryMarshal.Cast<byte, int>(this.segmentYStart.Memory.Span));
+    }
+
+    private void Apply(bool add, JxlImage3F opsin, Rectangle opsinRect)
+    {
+        if (this.segments.Count == 0)
+        {
+            return;
+        }
+
+        int y0 = RectangleUtils.Y0(in opsinRect);
+        int x0 = RectangleUtils.X0(in opsinRect);
+        int x1 = RectangleUtils.X1(in opsinRect);
+
+        for (int y = 0; y < opsinRect.Height; y++)
+        {
+            this.ApplyToRow(
+                add,
+                opsin.PlaneRowMemory(0, y0 + y)[x0..],
+                opsin.PlaneRowMemory(1, y0 + y)[x0..],
+                opsin.PlaneRowMemory(2, y0 + y)[x0..],
+                y0 + y,
+                x0,
+                x1);
+        }
+    }
+
+    private sealed class EmptyMemoryOwner : IMemoryOwner<byte>
+    {
+        public static readonly EmptyMemoryOwner Instance = new();
+
+        public Memory<byte> Memory => Memory<byte>.Empty;
+
+        public void Dispose()
+        {
+        }
+    }
+}

tests/ImageSharp.Tests/Formats/Jxl/Processing/AnsCommonTests.cs

@@ -0,0 +1,48 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using SixLabors.ImageSharp.Formats.Jxl.IO.Entropy;
+
+namespace SixLabors.ImageSharp.Tests.Formats.Jxl.Processing;
+
+public class AnsCommonTests
+{
+    private static void VerifyAliasDistribution(Span<int> distribution, uint logRange)
+    {
+        const int logAlphaSize = 8;
+
+        Span<JxlAnsEntry> table = stackalloc JxlAnsEntry[1 << logAlphaSize];
+        bool success = JxlAnsHelper.InitAliasTable(distribution, logRange, logAlphaSize, table);
+        Assert.True(success);
+
+        uint range = 1u << (int)logRange;
+        List<int>[] offsets = new List<int>[distribution.Length];
+
+        for (int i = 0; i < range; i++)
+        {
+            JxlAnsSymbol s = JxlAnsHelper.Lookup(table, i, JxlAnsConstants.AnsLogTableSize - 8, (1 << (JxlAnsConstants.AnsLogTableSize - 8)) - 1);
+
+            offsets[s.Value] ??= [];
+            offsets[s.Value].Add(s.Offset);
+        }
+
+        for (int i = 0; i < distribution.Length; i++)
+        {
+            Assert.Equal(distribution[i], offsets[i].Count);
+            offsets[i].Sort();
+
+            for (int j = 0; j < offsets[i].Count; j++)
+            {
+                Assert.Equal(offsets[i][j], j);
+            }
+        }
+    }
+
+    [Fact]
+    public void AliasDistributionSmoke()
+    {
+        VerifyAliasDistribution([JxlAnsConstants.AnsTableSize / 2, JxlAnsConstants.AnsTableSize / 2], JxlAnsConstants.AnsLogTableSize);
+        VerifyAliasDistribution([JxlAnsConstants.AnsTableSize], JxlAnsConstants.AnsLogTableSize);
+        VerifyAliasDistribution([0, 0, 0, JxlAnsConstants.AnsTableSize, 0], JxlAnsConstants.AnsLogTableSize);
+    }
+}

tests/ImageSharp.Tests/Formats/Jxl/Processing/Encoder/GammaCorrectionTests.cs

@@ -0,0 +1,37 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Encoder;
+
+namespace SixLabors.ImageSharp.Tests.Formats.Jxl.Processing.Encoder;
+
+public class GammaCorrectionTests
+{
+    [Fact]
+    public void TestLinearToSRgbEdgeCases()
+    {
+        Assert.Equal(0, JxlGammaCorrect.LinearToSRgb8Direct(0.0));
+
+        Assert.True(new TolerantMath(2E-5).AreEqual(0, JxlGammaCorrect.LinearToSRgb8Direct(1E-6)));
+
+        Assert.Equal(0, JxlGammaCorrect.LinearToSRgb8Direct(-1E-6));
+        Assert.Equal(0, JxlGammaCorrect.LinearToSRgb8Direct(-1E6));
+
+        Assert.True(new TolerantMath(1E-5).AreEqual(1, JxlGammaCorrect.LinearToSRgb8Direct(1 - 1E-6)));
+
+        Assert.Equal(1, JxlGammaCorrect.LinearToSRgb8Direct(1 + 1E-6));
+        Assert.Equal(1, JxlGammaCorrect.LinearToSRgb8Direct(1E6));
+    }
+
+    [Fact]
+    public void TestRoundTrip()
+    {
+        for (double linear = 0.0; linear <= 1.0; linear += 1E-7)
+        {
+            double srgb = JxlGammaCorrect.LinearToSRgb8Direct(linear);
+            double linear2 = JxlGammaCorrect.SRgb8ToLinearDirect(srgb);
+
+            Assert.True(Math.Abs(linear - linear2) < 2E-13, $"Linear = {linear}, Linear2 = {linear2}");
+        }
+    }
+}

tests/ImageSharp.Tests/Formats/Jxl/Processing/Encoder/Noise/PhotonNoiseTests.cs

@@ -0,0 +1,43 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using System.Diagnostics.CodeAnalysis;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Encoder.Noise;
+
+namespace SixLabors.ImageSharp.Tests.Formats.Jxl.Processing.Encoder.Noise;
+
+/// <summary>
+/// Tests for the simulation of photon noise in the JPEG XL encoder routines.
+/// </summary>
+public class PhotonNoiseTests
+{
+    [Fact]
+    public void TestPhotonNoiseEncoder()
+    {
+        ApproximateFloatComparer comparer = new(1e-6f);
+
+        Assert.Equal(
+            JxlPhotonNoise.SimulatePhotonNoise(xSize: 6000, ySize: 4000, iso: 100).Lookup,
+            [0.00259652f, 0.0139648f, 0.00681551f, 0.00632582f,
+             0.00694917f, 0.00803922f, 0.00934574f, 0.0107607f],
+            comparer);
+
+        Assert.Equal(
+            JxlPhotonNoise.SimulatePhotonNoise(xSize: 6000, ySize: 4000, iso: 800).Lookup,
+            [0.02077220f, 0.0420923f, 0.01820690f, 0.01439020f,
+             0.01293670f, 0.01254030f, 0.01277390f, 0.0134161f],
+            comparer);
+
+        Assert.Equal(
+            JxlPhotonNoise.SimulatePhotonNoise(xSize: 6000, ySize: 4000, iso: 6400).Lookup,
+            [0.1661770f, 0.1691120f, 0.05309080f, 0.03963960f,
+             0.03357410f, 0.03001650f, 0.02776740f, 0.0263478f],
+            comparer);
+
+        Assert.Equal(
+            JxlPhotonNoise.SimulatePhotonNoise(xSize: 4000, ySize: 3000, iso: 6400).Lookup,
+            [0.0830886f, 0.1008720f, 0.0367748f, 0.0280305f, 0.0240236f,
+             0.0218040f, 0.0205771f, 0.0200058f],
+            comparer);
+    }
+}

tests/ImageSharp.Tests/Formats/Jxl/Processing/Primitives/LehmerCodeTests.cs

@@ -0,0 +1,88 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using SixLabors.ImageSharp.Formats.Jxl.Processing;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Primitives;
+
+namespace SixLabors.ImageSharp.Tests.Formats.Jxl.Processing.Primitives;
+
+public class LehmerCodeTests
+{
+    private sealed class WorkingSet(int maxN)
+    {
+        public int PaddedN { get; } = maxN << JxlMath.CeilLog2Nonzero(maxN + 1);
+
+        public uint[] Permutation { get; } = new uint[maxN];
+
+        public uint[] Temporary { get; } = new uint[maxN];
+
+        public uint[] LehmerCodes { get; } = new uint[maxN];
+
+        public uint[] Decoded { get; } = new uint[maxN];
+    }
+
+    private static void RoundTrip(int n, WorkingSet ws)
+    {
+        Assert.NotEqual(0, n);
+        int paddedN = 1 << JxlMath.CeilLog2Nonzero(n);
+
+        Rng rng = new(((ulong)n * 65537) + 13);
+        Assert.True(n < 1 << (sizeof(uint) * 8));
+
+        Span<uint> permutationsSpan = ws.Permutation.AsSpan();
+        JxlSimdUtils.Iota(permutationsSpan[..n], 0u);
+
+        for (int rep = 0; rep < 3; rep++)
+        {
+            rng.Shuffle(permutationsSpan[..n]);
+
+            Assert.True(
+                JxlLehmerCode.ComputeLehmerCode(permutationsSpan, ws.Temporary, n, ws.LehmerCodes),
+                "Could not compute Lehmer code");
+
+            ws.Temporary.AsSpan()[..(paddedN * 4)].Clear();
+
+            Assert.True(
+                JxlLehmerCode.DecodeLehmerCode(ws.LehmerCodes.AsSpan(), ws.Temporary.AsSpan(), n, ws.Decoded.AsSpan()),
+                "Could not decode Lehmer code");
+
+            for (int i = 0; i < n; ++i)
+            {
+                Assert.Equal(permutationsSpan[i], ws.Decoded[i]);
+            }
+        }
+    }
+
+    private static void RoundTripSizeRange(int begin, int end)
+    {
+        Assert.NotEqual(0, begin);
+        List<WorkingSet> workingSets = [];
+
+        int numThreads = Environment.ProcessorCount;
+
+        // initialization
+        for (int i = 0; i < numThreads; i++)
+        {
+            workingSets.Add(new WorkingSet(end - 1));
+        }
+
+        // loop
+        Parallel.For(
+            begin,
+            end,
+            () => new WorkingSet(end - 1),
+            (n, _, workingSet) =>
+            {
+                RoundTrip(n, workingSet);
+                return workingSet;
+            },
+            _ => { });
+    }
+
+    [Fact]
+    public void TestLehmerCodes()
+    {
+        RoundTripSizeRange(1, 1026);
+        RoundTripSizeRange(65536, 65540);
+    }
+}

tests/ImageSharp.Tests/Formats/Jxl/Processing/Primitives/XorShiftTests.cs

@@ -0,0 +1,396 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using System.Diagnostics;
+using System.Numerics;
+using System.Runtime.InteropServices;
+using SixLabors.ImageSharp.Formats.Jxl.Processing.Primitives;
+using static SixLabors.ImageSharp.Tests.TestImages;
+
+namespace SixLabors.ImageSharp.Tests.Formats.Jxl.Processing.Primitives;
+
+public class XorShiftTests
+{
+    private const int Vectors = 64;
+
+    private static readonly ulong[][] ExpectedVectors =
+    [
+        [0x6E901576D477CBB1uL, 0xE9E53789195DA2A2uL, 0xB681F6DDA5E0AE99uL,
+     0x8EFD18CE21FD6896uL, 0xA898A80DF75CF532uL, 0x50CEB2C9E2DE7E32uL,
+     0x3CA7C2FEB25C0DD0uL, 0xA4D0866B80B4D836uL],
+    [0x8CD6A1E6233D3A26uL, 0x3D4603ADE98B112DuL, 0xDC427AF674019E36uL,
+     0xE28B4D230705AC53uL, 0x7297E9BBA88783DDuL, 0x34D3D23CFCD9B41AuL,
+     0x5A223615ADBE96B8uL, 0xE5EB529027CFBD01uL],
+    [
+    0xC1894CF00DFAC6A2uL, 0x18EDF8AE9085E404uL, 0x8E936625296B4CCDuL,
+     0x31971EF3A14A899BuL, 0xBE87535FCE0BF26AuL, 0x576F7A752BC6649FuL,
+     0xA44CBADCE0C6B937uL, 0x3DBA819BB17A353AuL],
+    [
+    0x27CE38DFCC1C5EB6uL, 0x920BEB5606340256uL, 0x3986CBC40C9AFC2CuL,
+     0xE22BCB3EEB1E191EuL, 0x6E1FCDD3602A8FBAuL, 0x052CB044E5415A29uL,
+     0x46266646EFB9ECD7uL, 0x8F44914618D29335uL],
+    [
+    0xDD30AEDF72A362C5uL, 0xBC1D824E16BB98F4uL, 0x9EA6009C2AA3D2F1uL,
+     0xF65C0FBBE17AF081uL, 0x22424D06A8738991uL, 0x8A62763F2B7611D2uL,
+     0x2F3E89F722637939uL, 0x84D338BEF50AFD50uL],
+    [
+    0x00F46494898E2B0BuL, 0x81239DC4FB8E8003uL, 0x414AD93EC5773FE7uL,
+     0x791473C450E4110FuL, 0x87F127BF68C959ACuL, 0x6429282D695EF67BuL,
+     0x661082E11546CBA8uL, 0x5815D53FA5436BFDuL],
+    [
+    0xB3DEADAB9BE6E0F9uL, 0xAA1B7B8F7CED0202uL, 0x4C5ED437699D279EuL,
+     0xA4471727F1CB39D3uL, 0xE439DA193F802F70uL, 0xF89401BB04FA6493uL,
+     0x3B08045A4FE898BAuL, 0x32137BFE98227950uL],
+    [
+    0xFBAE4A092897FEF3uL, 0x0639F6CE56E71C8EuL, 0xF0AD6465C07F0C1EuL,
+     0xFF8E28563361DCE5uL, 0xC2013DB7F86BC6B9uL, 0x8EFCC0503330102FuL,
+     0x3F6B767EA5C4DA40uL, 0xB9864B950B2232E1uL],
+    [
+    0x76EB58DE8E5EC22AuL, 0x9BBBF49A18B32F4FuL, 0xC8405F02B2B2FAB9uL,
+     0xC3E122A5F146BC34uL, 0xC90BB046660F5765uL, 0xB933981310DBECCFuL,
+     0x5A2A7BFC9126FD1CuL, 0x8BB388C94DF87901uL],
+    [
+    0x753EB89AD63EF3C3uL, 0xF24AAF40C89D65ADuL, 0x23F68931C1A6AA6DuL,
+     0xF47E79BF702C6DD0uL, 0xA3AD113244EE7EAEuL, 0xD42CBEA28F793DC3uL,
+     0xD896FCF1820F497CuL, 0x042B86D2818948C1uL],
+    [
+    0x8F2A4FC5A4265763uL, 0xEC499E6F95EAA10CuL, 0xE3786D4ECCD0DEB5uL,
+     0xC725C53D3AC4CC43uL, 0x065A4ACBBF83610EuL, 0x35C61C9FEF167129uL,
+     0x7B720AEAA7D70048uL, 0x14206B841377D039uL],
+    [
+    0xAD27D78BF96055F6uL, 0x5F43B20FF47ADCD4uL, 0xE184C2401E2BF71EuL,
+     0x30B263D78990045DuL, 0xC22F00EBFF9BA201uL, 0xAE7F86522B53A562uL,
+     0x2853312BC039F0A4uL, 0x868D619E6549C3C8uL],
+    [
+    0xFD5493D8AE9A8371uL, 0x773D5E224DF61B3BuL, 0x5377C54FBB1A8280uL,
+     0xCAD4DE3B8265CAFAuL, 0xCDF3F19C91EBD5F6uL, 0xC8EA0F182D73BD78uL,
+     0x220502D593433FF1uL, 0xB81205E612DC31B1uL],
+    [
+    0x8F32A39EAEDA4C70uL, 0x1D4B0914AA4DAC7FuL, 0x56EF1570F3A8B405uL,
+     0x29812CB17404A592uL, 0x97A2AAF69CAE90F2uL, 0x12BF5E02778BBFE5uL,
+     0x9D4B55AD42A05FD2uL, 0x06C2BAB5E6086620uL],
+    [
+    0x8DB4B9648302B253uL, 0xD756AD9E3AEA12C7uL, 0x68709B7F11D4B188uL,
+     0x7CC299DDCD707A4BuL, 0x97B860C370A7661DuL, 0xCECD314FC20E64F5uL,
+     0x55F412CDFB4C7EC3uL, 0x55EE97591193B525uL],
+    [
+    0xCF70F3ACA96E6254uL, 0x022FEDECA2E09F46uL, 0x686823DB60AE1ECFuL,
+     0xFD36190D3739830EuL, 0x74E1C09027F68120uL, 0xB5883A835C093842uL,
+     0x93E1EFB927E9E4E3uL, 0xB2721E249D7E5EBEuL],
+    [
+    0x69B6E21C44188CB8uL, 0x5D6CFB853655A7AAuL, 0x3E001A0B425A66DCuL,
+     0x8C57451103A5138FuL, 0x7BF8B4BE18EAB402uL, 0x494102EB8761A365uL,
+     0xB33796A9F6A81F0EuL, 0x10005AB3BCCFD960uL],
+    [
+    0xB2CF25740AE965DCuL, 0x6F7C1DF7EF53D670uL, 0x648DD6087AC2251EuL,
+     0x040955D9851D487DuL, 0xBD550FC7E21A7F66uL, 0x57408F484DEB3AB5uL,
+     0x481E24C150B506C1uL, 0x72C0C3EAF91A40D6uL],
+    [
+    0x1997A481858A5D39uL, 0x539718F4BEF50DC1uL, 0x2EC4DC4787E7E368uL,
+     0xFF1CE78879419845uL, 0xE219A93DD6F6DD30uL, 0x85328618D02FEC1AuL,
+     0xC86E02D969181B20uL, 0xEBEC8CD8BBA34E6EuL],
+    [
+    0x28B55088A16CE947uL, 0xDD25AC11E6350195uL, 0xBD1F176694257B1CuL,
+     0x09459CCF9FCC9402uL, 0xF8047341E386C4E4uL, 0x7E8E9A9AD984C6C0uL,
+     0xA4661E95062AA092uL, 0x70A9947005ED1152uL],
+    [
+    0x4C01CF75DBE98CCDuL, 0x0BA076CDFC7373B9uL, 0x6C5E7A004B57FB59uL,
+     0x336B82297FD3BC56uL, 0x7990C0BE74E8D60FuL, 0xF0275CC00EC5C8C8uL,
+     0x6CF29E682DFAD2E9uL, 0xFA4361524BD95D72uL],
+    [
+    0x631D2A19FF62F018uL, 0x41C43863B985B3FAuL, 0xE052B2267038EFD9uL,
+     0xE2A535FAC575F430uL, 0xE004EEA90B1FF5B8uL, 0x42DFE2CA692A1F26uL,
+     0x90FB0BFC9A189ECCuL, 0x4484102BD3536BD0uL],
+    [
+    0xD027134E9ACCA5A5uL, 0xBBAB4F966D476A9BuL, 0x713794A96E03D693uL,
+     0x9F6335E6B94CD44AuL, 0xC5090C80E7471617uL, 0x6D9C1B0C87B58E33uL,
+     0x1969CE82E31185A5uL, 0x2099B97E87754EBEuL],
+    [
+    0x60EBAF4ED934350FuL, 0xC26FBF0BA5E6ECFFuL, 0x9E54150F0312EC57uL,
+     0x0973B48364ED0041uL, 0x800A523241426CFCuL, 0x03AB5EC055F75989uL,
+     0x8CF315935DEEB40AuL, 0x83D3FC0190BD1409uL],
+    [
+    0x26D35394CF720A51uL, 0xCE9EAA15243CBAFEuL, 0xE2B45FBAF21B29E0uL,
+     0xDB92E98EDE73F9E0uL, 0x79B16F5101C26387uL, 0x1AC15959DE88C86FuL,
+     0x387633AEC6D6A580uL, 0xA6FC05807BFC5EB8uL],
+    [
+    0x2D26C8E47C6BADA9uL, 0x820E6EC832D52D73uL, 0xB8432C3E0ED0EE5BuL,
+     0x0F84B3C4063AAA87uL, 0xF393E4366854F651uL, 0x749E1B4D2366A567uL,
+     0x805EACA43480D004uL, 0x244EBF3AA54400A5uL],
+    [
+    0xBFDC3763AA79F75AuL, 0x9E3A74CC751F41DBuL, 0xF401302A149DBC55uL,
+     0x6B25F7973D7BF7BCuL, 0x13371D34FDBC3DAEuL, 0xC5E1998C8F484DCDuL,
+     0x7031B8AE5C364464uL, 0x3847F0C4F3DA2C25uL],
+    [
+    0x24C6387D2C0F1225uL, 0x77CCE960255C67A4uL, 0x21A0947E497B10EBuL,
+     0xBB5DB73A825A9D7EuL, 0x26294A41999E553DuL, 0x3953E0089F87D925uL,
+     0x3DAE6E5D4E5EAAFEuL, 0x74B545460341A7AAuL],
+    [
+    0x710E5EB08A7DB820uL, 0x7E43C4E77CAEA025uL, 0xD4C91529C8B060C1uL,
+     0x09AE26D8A7B0CA29uL, 0xAB9F356BB360A772uL, 0xB68834A25F19F6E9uL,
+     0x79B8D9894C5734E2uL, 0xC6847E7C8FFD265FuL],
+    [
+    0x10C4BCB06A5111E6uL, 0x57CB50955B6A2516uL, 0xEF53C87798B6995FuL,
+     0xAB38E15BBD8D0197uL, 0xA51C6106EFF73C93uL, 0x83D7F0E2270A7134uL,
+     0x0923FD330397FCE5uL, 0xF9DE54EDFE58FB45uL],
+    [
+    0x07D44833ACCD1A94uL, 0xAAD3C9E945E2F9F3uL, 0xABF4C879B876AA37uL,
+     0xF29C69A21B301619uL, 0x2DDCE959111C788BuL, 0x7CEDB48F8AC1729BuL,
+     0x93F3BA9A02B659BEuL, 0xF20A87FF17933CBEuL],
+    [
+    0x8E96EBE93180CFE6uL, 0x94CAA12873937079uL, 0x05F613D9380D4189uL,
+     0xBCAB40C1DC79F38AuL, 0x0AD8907B7C61D19EuL, 0x88534E189D103910uL,
+     0x2DB2FAABA160AB8FuL, 0xA070E7506B06F15CuL],
+    [
+    0x6FB1FCDAFFEF87A9uL, 0xE735CF25337A090DuL, 0x172C6EDCEFEF1825uL,
+     0x76957EA49EF0542DuL, 0x819BF4CD250F7C49uL, 0xD6FF23E4AD00C4D4uL,
+     0xE79673C1EC358FF0uL, 0xAC9C048144337938uL],
+    [
+    0x4C5387FF258B3AF4uL, 0xEDB68FAEC2CB1AA3uL, 0x02A624E67B4E1DA4uL,
+     0x5C44797A38E08AF2uL, 0x36546A70E9411B4BuL, 0x47C17B24D2FD9675uL,
+     0x101957AAA020CA26uL, 0x47A1619D4779F122uL],
+    [
+    0xF84B8BCDC92D9A3CuL, 0x951D7D2C74B3066BuL, 0x7AC287C06EDDD9B2uL,
+     0x4C38FC476608D38FuL, 0x224D793B19CB4BCDuL, 0x835A255899BF1A41uL,
+     0x4AD250E9F62DB4ABuL, 0xD9B44F4B58781096uL],
+    [
+    0xABBAF99A8EB5C6B8uL, 0xFB568E900D3A9F56uL, 0x11EDF63D23C5DF11uL,
+     0xA9C3011D3FA7C5A8uL, 0xAEDD3CF11AFFF725uL, 0xABCA472B5F1EDD6BuL,
+     0x0600B6BB5D879804uL, 0xDB4DE007F22191A0uL],
+    [
+    0xD76CC9EFF0CE9392uL, 0xF5E0A772B59BA49AuL, 0x7D1AE1ED0C1261B5uL,
+     0x79224A33B5EA4F4AuL, 0x6DD825D80C40EA60uL, 0x47FC8E747E51C953uL,
+     0x695C05F72888BF98uL, 0x1A012428440B9015uL],
+    [
+    0xD754DD61F9B772BFuL, 0xC4A2FCF4C0F9D4EBuL, 0x461167CDF67A24A2uL,
+     0x434748490EBCB9D4uL, 0x274DD9CDCA5781DEuL, 0x36BAC63BA9A85209uL,
+     0x30324DAFDA36B70FuL, 0x337570DB4FE6DAB3uL],
+    [
+    0xF46CBDD57C551546uL, 0x8E02507E676DA3E3uL, 0xD826245A8C15406DuL,
+     0xDFB38A5B71113B72uL, 0x5EA38454C95B16B5uL, 0x28C054FB87ABF3E1uL,
+     0xAA2724C0BA1A8096uL, 0xECA83EC980304F2FuL],
+    [
+    0x6AA76EC294EB3303uL, 0x42D4CDB2A8032E3BuL, 0x7999EDF75DCD8735uL,
+     0xB422BFFE696CCDCCuL, 0x8F721461FD7CCDFEuL, 0x148E1A5814FDE253uL,
+     0x4DC941F4375EF8FFuL, 0x27B2A9E0EB5B49CFuL],
+    [
+    0xCEA592EF9343EBE1uL, 0xF7D38B5FA7698903uL, 0x6CCBF352203FEAB6uL,
+     0x830F3095FCCDA9C5uL, 0xDBEEF4B81B81C8F4uL, 0x6D7EB9BCEECA5CF9uL,
+     0xC58ABB0FBE436C69uL, 0xE4B97E6DB2041A4BuL],
+    [
+    0x7E40FC772978AF14uL, 0xCDDA4BBAE28354A1uL, 0xE4F993B832C32613uL,
+     0xD3608093C68A4B35uL, 0x9A3B60E01BEE3699uL, 0x03BEF248F3288713uL,
+     0x70B9294318F3E9B4uL, 0x8D2ABB913B8610DEuL],
+    [
+    0x37F209128E7D8B2CuL, 0x81D2AB375BD874BCuL, 0xA716A1B7373F7408uL,
+     0x0CEE97BEC4706540uL, 0xA40C5FD9CDBC1512uL, 0x73CAF6C8918409E7uL,
+     0x45E11BCEDF0BBAA1uL, 0x612C612BFF6E6605uL],
+    [
+    0xF8ECB14A12D0F649uL, 0xDA683CD7C01BA1ACuL, 0xA2203F7510E124C1uL,
+     0x7F83E52E162F3C78uL, 0x77D2BB73456ACADBuL, 0x37FC34FC840BBA6FuL,
+     0x3076BC7D4C6EBC1FuL, 0x4F514123632B5FA9uL],
+    [
+    0x44D789DED935E884uL, 0xF8291591E09FEC9FuL, 0xD9CED2CF32A2E4B7uL,
+     0x95F70E1EB604904AuL, 0xDE438FE43C14F6ABuL, 0x4C8D23E4FAFCF8D8uL,
+     0xC716910A3067EB86uL, 0x3D6B7915315095D3uL],
+    [
+    0x3170FDBADAB92095uL, 0x8F1963933FC5650BuL, 0x72F94F00ABECFEABuL,
+     0x6E3AE826C6AAB4CEuL, 0xA677A2BF31068258uL, 0x9660CDC4F363AF10uL,
+     0xD81A15A152379EF1uL, 0x5D7D285E1080A3F9uL],
+    [
+    0xDAD5DDFF9A2249B3uL, 0x6F9721D926103FAEuL, 0x1418CBB83FFA349AuL,
+     0xE71A30AD48C012B2uL, 0xBE76376C63751132uL, 0x3496467ACA713AE6uL,
+     0x8D7EC01369F991A3uL, 0xD8C73A88B96B154EuL],
+    [
+    0x8B5D9C74AEB4833AuL, 0xF914FB3F867B912FuL, 0xB894EA034936B1DCuL,
+     0x8A16D21BE51C4F5BuL, 0x31FF048ED582D98EuL, 0xB95AB2F4DC65B820uL,
+     0x04082B9170561AF7uL, 0xA215610A5DC836FAuL],
+    [
+    0xB2ADE592C092FAACuL, 0x7A1E683BCBF13294uL, 0xC7A4DBF86858C096uL,
+     0x3A49940F97BFF316uL, 0xCAE5C06B82C46703uL, 0xC7F413A0F951E2BDuL,
+     0x6665E7BB10EB5916uL, 0x86F84A5A94EDE319uL],
+    [
+    0x4EA199D8FAA79CA3uL, 0xDFA26E5BF1981704uL, 0x0F5E081D37FA4E01uL,
+     0x9CB632F89CD675CDuL, 0x4A09DB89D48C0304uL, 0x88142742EA3C7672uL,
+     0xAC4F149E6D2E9BDBuL, 0x6D9E1C23F8B1C6C6uL],
+    [
+    0xD58BE47B92DEC0E9uL, 0x8E57573645E34328uL, 0x4CC094CCB5FB5126uL,
+     0x5F1D66AF6FB40E3CuL, 0x2BA15509132D3B00uL, 0x0D6545646120E567uL,
+     0x3CF680C45C223666uL, 0x96B28E32930179DAuL],
+    [
+    0x5900C45853AC7990uL, 0x61881E3E8B7FF169uL, 0x4DE5F835DF2230FFuL,
+     0x4427A9E7932F73FFuL, 0x9B641BAD379A8C8DuL, 0xDF271E5BF98F4E5CuL,
+     0xDFDA16DB830FF5EEuL, 0x371C7E7CFB89C0E9uL],
+    [
+    0x4410A8576247A250uL, 0x6AD2DA12B45AC0D9uL, 0x18DFC72AAC85EECCuL,
+     0x06FC8BB2A0EF25C8uL, 0xEB287619C85E6118uL, 0x19553ECA67F25A2CuL,
+     0x3B9557F1DCEC5BAAuL, 0x7BAD9E8B710D1079uL],
+    [
+    0x34F365D66BD22B28uL, 0xE6E124B9F10F835DuL, 0x0573C38ABF2B24DCuL,
+     0xD32E6AF10A0125AEuL, 0x383590ACEA979519uL, 0x8376ED7A39E28205uL,
+     0xF0B7F184DCBDA435uL, 0x062A203390E31794uL],
+    [
+    0xA2AFFD7E41918760uL, 0x7F90FC1BD0819C86uL, 0x5033C08E5A969533uL,
+     0x2707AF5C6D039590uL, 0x57BBD5980F17DF9CuL, 0xD3FE6E61D763268AuL,
+     0x9E0A0AE40F335A3BuL, 0x43CF4EB0A99613C5uL],
+    [
+    0xD4D2A397CE1A7C2EuL, 0x3DF7CE7CC3212DADuL, 0x0880F0D5D356C75AuL,
+     0xA8AFC44DD03B1346uL, 0x79263B46C13A29E0uL, 0x11071B3C0ED58E7AuL,
+     0xED46DC9F538406BFuL, 0x2C94974F2B94843DuL],
+    [
+    0xE246E13C39AB5D5EuL, 0xAC1018489D955B20uL, 0x8601B558771852B8uL,
+     0x110BD4C06DB40173uL, 0x738FC8A18CCA0EBBuL, 0x6673E09BE0EA76E5uL,
+     0x024BC7A0C7527877uL, 0x45E6B4652E2EC34EuL],
+    [
+    0xD1ED26A1A375CDC8uL, 0xAABC4E896A617CB8uL, 0x0A9C9E8E57D753C6uL,
+     0xA3774A75FEB4C30EuL, 0x30B816C01C93E49EuL, 0xF405BABC06D2408CuL,
+     0xCC0CE6B4CE788ABCuL, 0x75E7922D0447956CuL],
+    [
+    0xD07C1676A698BC95uL, 0x5F9AEA4840E2D860uL, 0xD5FC10D58BDF6F02uL,
+     0xF190A2AD4BC2EEA7uL, 0x0C24D11F51726931uL, 0xDB646899A16B6512uL,
+     0x7BC10670047B1DD8uL, 0x2413A5ABCD45F092uL],
+    [
+    0x4E66892190CFD923uL, 0xF10162440365EC8EuL, 0x158ACA5A6A2280AEuL,
+     0x0D60ED11C0224166uL, 0x7CD2E9A71B9D7488uL, 0x450D7289706AB2A3uL,
+     0x88FAE34EC9A0D7DCuL, 0x96FF9103575A97DAuL],
+    [
+    0x77990FAC6046C446uL, 0xB174B5FB30C76676uL, 0xE352CE3EB56CF82AuL,
+     0xC6039B6873A9A082uL, 0xE3F80F3AE333148AuL, 0xB853BA24BA3539B9uL,
+     0xE8863E52ECCB0C74uL, 0x309B4CC1092CC245uL],
+    [
+    0xBC2B70BEE8388D9FuL, 0xE48D92AE22216DCEuL, 0xF15F3BF3E2C15D8FuL,
+     0x1DD964D4812D8B24uL, 0xD56AF02FB4665E4CuL, 0x98002200595BD9A3uL,
+     0x049246D50BB8FA12uL, 0x1B542DF485B579B9uL],
+    [
+    0x2347409ADFA8E497uL, 0x36015C2211D62498uL, 0xE9F141F32EB82690uL,
+     0x1F839912D0449FB9uL, 0x4E4DCFFF2D02D97CuL, 0xF8A03AB4C0F625C9uL,
+     0x0605F575795DAC5CuL, 0x4746C9BEA0DDA6B1uL],
+    [
+    0xCA5BB519ECE7481BuL, 0xFD496155E55CA945uL, 0xF753B9DBB1515F81uL,
+     0x50549E8BAC0F70E7uL, 0x8614FB0271E21C60uL, 0x60C72947EB0F0070uL,
+     0xA6511C10AEE742B6uL, 0x48FB48F2CACCB43EuL]
+    ];
+
+    [Fact]
+    public void TestGolden()
+    {
+        JxlXorShift rng = new(12345);
+        Span<ulong> lanes = stackalloc ulong[JxlXorShift.Generators];
+
+        for (ulong vector = 0; vector < Vectors; vector++)
+        {
+            rng.Fill(lanes);
+
+            for (int i = 0; i < JxlXorShift.Generators; i++)
+            {
+                Assert.Equal(ExpectedVectors[(int)vector][i], lanes[i]);
+            }
+        }
+    }
+
+    [Fact]
+    public void TestSeedChanges()
+    {
+        Span<ulong> lanes = stackalloc ulong[JxlXorShift.Generators];
+        const int numberOfSeeds = 16384;
+
+        List<ulong> first = new(numberOfSeeds);
+        for (int seed = 0; seed < numberOfSeeds; seed++)
+        {
+            JxlXorShift xs128Plus = new((ulong)seed);
+            xs128Plus.Fill(lanes);
+            first.Add(lanes[0]);
+        }
+
+        Assert.Equal(numberOfSeeds, first.Count);
+        first.Sort();
+        first = [.. first.Distinct()];
+        Assert.Equal(numberOfSeeds, first.Count);
+    }
+
+    [Fact]
+    public void TestFloat()
+    {
+#if ALLOW_JPEGXL_SLOW_TESTS
+        const int seedMax = 4096;
+#else
+        const int seedMax = 256;
+#endif
+
+        Parallel.For(0, seedMax, seed =>
+        {
+            const int VecCap = 16;
+
+            JxlXorShift rng = new((ulong)seed);
+
+            Span<ulong> batch64 = stackalloc ulong[JxlXorShift.Generators];
+            Span<uint> batch32 = stackalloc uint[2 * JxlXorShift.Generators];
+
+            Span<float> lanes = stackalloc float[Vector<float>.Count];
+
+            int count = 0;
+            const int reps = 32000;
+            double sum = 0.0;
+            while (count < reps)
+            {
+                rng.Fill(batch64);
+
+                MemoryMarshal.Cast<ulong, uint>(batch64).CopyTo(batch32);
+
+                for (int i = 0; i < VecCap; i += Vector<float>.Count)
+                {
+                    Vector<uint> bits = new(batch32.Slice(i, Vector<uint>.Count));
+
+                    // (bits >> 9) | 0x3F800000
+                    Vector<uint> shifted = bits >> 9;
+                    Vector<uint> mantissa = shifted | new Vector<uint>(0x3F800000);
+
+                    Vector<float> rand12 = Vector.AsVectorSingle(mantissa);
+
+                    rand12.CopyTo(lanes);
+
+                    for (int j = 0; j < Vector<float>.Count; j++)
+                    {
+                        float lane = lanes[j];
+                        sum += lane;
+                        count++;
+
+                        Assert.True(lane < 2.0f);
+                        Assert.True(lane >= 1.0f);
+                    }
+                }
+            }
+        });
+    }
+
+    [Fact]
+    public void TestNotZero()
+    {
+#if ALLOW_JPEGXL_SLOW_TESTS
+        const int seedMax = 2000;
+#else
+        const int seedMax = 500;
+#endif
+
+        Parallel.For(0, seedMax, task =>
+        {
+            Span<ulong> lanes = stackalloc ulong[JxlXorShift.Generators];
+            JxlXorShift rng = new((ulong)task);
+            int numZero = 0;
+
+            for (int vectors = 0; vectors < 10000; vectors++)
+            {
+                rng.Fill(lanes);
+                for (int i = 0; i < lanes.Length; i++)
+                {
+                    if (lanes[i] == 0)
+                    {
+                        numZero++;
+                    }
+                }
+            }
+
+            Assert.True(numZero < 1, "There should not be any 0 values produced by the RNG");
+        });
+    }
+}

tests/ImageSharp.Tests/Formats/Jxl/README.md

@@ -0,0 +1,5 @@
+# JPEG XL Tests
+
+### Preprocessor directives
+Enable a ALLOW_JPEGXL_SLOW_TESTS preprocessor directive to allow
+slower tests that test the library more extensively.

tests/ImageSharp.Tests/Formats/Jxl/Rng.cs

@@ -0,0 +1,87 @@
+// Copyright (c) Six Labors.
+// Licensed under the Six Labors Split License.
+
+using System.Runtime.CompilerServices;
+using SixLabors.ImageSharp.Common.Helpers;
+
+namespace SixLabors.ImageSharp.Tests.Formats.Jxl;
+
+/// <summary>
+/// Deterministic random number generator used for compatibility
+/// with JPEG XL tests.
+/// </summary>
+internal struct Rng
+{
+    private ulong s0;
+    private ulong s1;
+
+    public Rng(ulong seed)
+    {
+        this.s0 = 0x94D049BB133111EBUL;
+        this.s1 = 0xBF58476D1CE4E5B9UL + seed;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public ulong Next()
+    {
+        ulong s1 = this.s0;
+        ulong s0 = this.s1;
+        ulong bits = s1 + s0;
+        this.s0 = s0;
+
+        s1 ^= s1 << 23;
+        s1 ^= s0 ^ (s1 >> 18) ^ (s0 >> 5);
+
+        this.s1 = s1;
+
+        return bits;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public long UniformI(long begin, long end) => (long)(this.Next() % (ulong)(end - begin)) + begin;
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public ulong UniformU(ulong begin, ulong end) => (this.Next() % (end - begin)) + begin;
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public float UniformF(float begin, float end)
+    {
+        uint u = (uint)(this.Next() >> (64 - 23)) | 0x3F800000u;
+        float f = BitConverter.UInt32BitsToSingle(u);
+
+        return ((end - begin) * (f - 1.0f)) + begin;
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public bool Bernoulli(float p) => this.UniformF(0, 1) < p;
+
+    internal readonly struct GeometricDistribution
+    {
+        public readonly float Value { get; }
+
+        public GeometricDistribution(float value) => this.Value = value;
+
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        public static GeometricDistribution Make(float p) => new(1.0f / MathF.Log(1.0f - p));
+    }
+
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
+    public uint Geometric(in GeometricDistribution dist)
+    {
+        float f = this.UniformF(0, 1);
+        float invLog1mp = dist.Value;
+
+        float log = MathF.Log(1.0f - f) * invLog1mp;
+
+        return (uint)log;
+    }
+
+    public void Shuffle<T>(Span<T> span)
+    {
+        for (nuint i = 0; i + 1 < (nuint)span.Length; i++)
+        {
+            nuint a = (nuint)this.UniformU(i, (nuint)span.Length);
+            RuntimeUtility.Swap(ref span[(int)a], ref span[(int)i]);
+        }
+    }
+}