Fast shuffle by batching

Shuffle.cs

using System.Numerics;
using System.Runtime.InteropServices;
using System.Security.Cryptography;
using System.Linq;
using System;
using System.Runtime.CompilerServices;


// Test code
var rng = new Seiran128();
var array = Enumerable.Range(0, 32).ToArray();
rng.Shuffle(array.AsSpan());
Console.WriteLine(string.Join(", ", array));


public interface IRandom
{
    public ulong Next();
}

public sealed class Seiran128 : IRandom
{
    private ulong state0, state1;

    public Seiran128()
    {
        var buffer = (stackalloc ulong[2]);
        RandomNumberGenerator.Fill(MemoryMarshal.AsBytes(buffer));

        state0 = buffer[0];
        state1 = buffer[1];
    }

    public ulong Next()
    {
        ulong s0 = state0, s1 = state1;

        ulong result = BitOperations.RotateLeft((s0 + s1) * 9, 29) + s0;

        state0 = s0 ^ BitOperations.RotateLeft(s1, 29);
        state1 = s0 ^ s1 << 9;

        return result;
    }
}

public static class RandomExtension
{
    [MethodImpl(MethodImplOptions.AggressiveInlining)]
    private static void Swap<T>(Span<T> span, int indexA, int indexB)
    {
        ref var head = ref MemoryMarshal.GetReference(span);
        var t = Unsafe.Add(ref head, indexA);
        Unsafe.Add(ref head, indexA) = Unsafe.Add(ref head, indexB);
        Unsafe.Add(ref head, indexB) = t;
    }

    public static void Shuffle<TRng, TSpan>(this TRng rng, Span<TSpan> span)
        where TRng : IRandom
    {
        ulong bound = 2432902008176640000;  // 20!;
        int k = Math.Min(20, span.Length);


        for (int i = 1; i < span.Length;)
        {
            ulong r = rng.Next();

            // correct r
            ulong rbound = r * bound;
            if (rbound > 0 - bound)
            {
                ulong sum, carry;
                do
                {
                    ulong r2 = rng.Next();
                    ulong lohi = Math.BigMul(r2, bound, out ulong lolo);
                    sum = rbound + lohi;
                    carry = sum < rbound ? 1ul : 0ul;
                    rbound = lolo;
                } while (sum == ~0ul);
                r += carry;
            }


            // run fisher-yates shuffle
            for (int m = i; m < k; m++)
            {
                int index = (int)Math.BigMul(r, (ulong)(m + 1), out r);

                Swap(span, m, index);
            }


            // calculates next `n`
            i = k;
            bound = (ulong)(i + 1);
            for (k = i + 1; k < span.Length; k++)
            {
                if (Math.BigMul(bound, (ulong)(k + 1), out var newbound) == 0)
                {
                    bound = newbound;
                }
                else
                {
                    break;
                }
            }
        }
    }
}