伪随机数生成器“Mersenne Vortex” - SIMD 实现

我需要一个分布良好且周期长的高质量 HRNG。现存最流行的是梅森漩涡。

我从 C 的作者SFMT 1.5.1那里获取了 SIMD实现， 并用 C# (.NET 6) 重写了它。

抱歉，游戏开发者，这个版本的 Whirlwind 与 Unity 不兼容，不容易修改，我们必须等到 Unity 本身切换到 CoreCLR 和 .NET 6。:(

完成了以下任务：

• 处理器必须支持 SSE2，无回退
• 支持静态成员级别的线程安全
• 与公众成员的行为兼容System.Random
• 没有零种病
• 让它尽可能简单
• 面向 x64 应用程序

/// <summary>
/// Генератор псевдослучайных чисел на основе алгоритма Вихрь Мерсенна
/// Конфигурация: SFMT-19937:122-18-1-11-1:dfffffef-ddfecb7f-bffaffff-bffffff6
/// </summary>
public class MersenneTwister
{
    const int MersenneExponent = 19937;
    const int Length128 = MersenneExponent / 128 + 1;
    const int Length32 = Length128 * 4;

    private readonly Vector128<uint>[] state = new Vector128<uint>[Length128];
    private int index;

    [ThreadStatic]
    private static MersenneTwister _shared;

    /// <summary>
    /// Общий экземпляр генератора, создается отдельно для каждого потока
    /// </summary>
    public static MersenneTwister Shared => _shared ??= new(GetRandomSeed());

    /// <summary>
    /// Создает экземпляр генератора со случайным сидом
    /// </summary>
    public MersenneTwister() 
    {
        InitGenerator(Shared.NextUInt32());
    }

    /// <summary>
    /// Создает экземпляр генератора с указанным сидом
    /// </summary>
    /// <param name="seed">Сид для создания генератора</param>
    public MersenneTwister(int seed) : this((uint)seed) { }

    private MersenneTwister(uint seed)
    {
        InitGenerator(seed);
    }

    private static uint GetRandomSeed()
    {
        ReadOnlySpan<byte> bytes = Guid.NewGuid().ToByteArray();
        ReadOnlySpan<uint> hash = MemoryMarshal.Cast<byte, uint>(SHA256.HashData(bytes));
        uint result = 0;
        foreach (uint value in hash)
            result ^= value;
        return result;
    }

    private void InitGenerator(uint seed)
    {
        if (!Sse2.IsSupported)
            throw new InvalidOperationException("SSE2 не поддерживается на данном устройстве.");

        Span<uint> values = MemoryMarshal.Cast<Vector128<uint>, uint>(state);
        values[0] = seed;
        for (int i = 1; i < Length32; i++)
        {
            values[i] = (uint)(1812433253UL * (values[i - 1] ^ (values[i - 1] >> 30)) + (uint)i);
        }
        index = Length32;
        PeriodCertification(values);
    }

    private void PeriodCertification(Span<uint> values)
    {
        uint inner = 0;
        ReadOnlySpan<uint> parity = stackalloc uint[] { 0x00000001U, 0x00000000U, 0x00000000U, 0x13c9e684U };

        for (int i = 0; i < 4; i++)
            inner ^= values[i] & parity[i];
        for (int i = 16; i > 0; i >>= 1)
            inner ^= inner >> i;
        inner &= 1;

        if (inner == 1)
            return;

        for (int i = 0; i < 4; i++)
        {
            for (uint work = 1; work != 0; work <<= 1)
            {
                if ((work & parity[i]) != 0)
                {
                    values[i] ^= work;
                    return;
                }
            }
        }
    }

    private void UpdateState()
    {
        const int offset = 122;
        int i;
        Vector128<uint> r1 = state[^2];
        Vector128<uint> r2 = state[^1];
        for (i = 0; i < Length128 - offset; i++)
        {
            state[i] = DoRecursion(state[i], state[i + offset], r1, r2);
            r1 = r2;
            r2 = state[i];
        }
        for (; i < Length128; i++)
        {
            state[i] = DoRecursion(state[i], state[i + offset - Length128], r1, r2);
            r1 = r2;
            r2 = state[i];
        }
    }

    private Vector128<uint> DoRecursion(Vector128<uint> a, Vector128<uint> b, Vector128<uint> c, Vector128<uint> d)
    {
        Vector128<uint> z = Sse2.ShiftRightLogical128BitLane(c, 1);
        z = Sse2.Xor(z, a);
        Vector128<uint> v = Sse2.ShiftLeftLogical(d, 18);
        z = Sse2.Xor(z, v);
        Vector128<uint> x = Sse2.ShiftLeftLogical128BitLane(a, 1);
        z = Sse2.Xor(z, x);
        Vector128<uint> y = Sse2.ShiftRightLogical(b, 11);
        Vector128<uint> mask = Vector128.Create(0xdfffffefU, 0xddfecb7fU, 0xbffaffffU, 0xbffffff6U);
        y = Sse2.And(y, mask);
        return Sse2.Xor(z, y);
    }

    /// <summary>
    /// Возвращает число в диапазоне [0,ulong.MaxValue]
    /// </summary>
    public ulong NextUInt64()
    {
        if (index >= Length32)
        {
            UpdateState();
            index = 0;
        }
        Span<ulong> values = MemoryMarshal.Cast<Vector128<uint>, ulong>(state);
        ulong r = values[index / 2];
        index += 2;
        return r;
    }

    /// <summary>
    /// Возвращает число в диапазоне [0,uint.MaxValue]
    /// </summary>
    public uint NextUInt32()
    {
        if (index >= Length32)
        {
            UpdateState();
            index = 0;
        }
        Span<uint> values = MemoryMarshal.Cast<Vector128<uint>, uint>(state);
        return values[index++];
    }

    /// <summary>
    /// Возвращает число в диапазоне [0,int.MaxValue)
    /// </summary>
    public int Next()
    {
        return (int)(((ulong)int.MaxValue * NextUInt32()) >> 32);
    }

    /// <summary>
    /// Возвращает число в диапазоне [0,maxValue)
    /// </summary>
    public int Next(int maxValue)
    {
        if (maxValue <= 0)
            throw new ArgumentOutOfRangeException(nameof(maxValue), "Значение должно быть больше 0.");

        return (int)(((ulong)maxValue * NextUInt32()) >> 32);
    }

    /// <summary>
    /// Возвращает число в диапазоне [minValue,maxValue)
    /// </summary>
    public int Next(int minValue, int maxValue)
    {
        if (minValue < 0)
            throw new ArgumentOutOfRangeException(nameof(minValue), "Значение должно быть больше либо равно 0.");
        if (maxValue <= minValue)
            throw new ArgumentOutOfRangeException(nameof(maxValue), $"Значение должно быть больше, чем {nameof(minValue)}.");

        return (int)(((ulong)(maxValue - minValue) * NextUInt32()) >> 32) + minValue;
    }

    /// <summary>
    /// Возвращает число в диапазоне [0,long.MaxValue)
    /// </summary>
    public long NextInt64()
    {
        return (long)(((BigInteger)long.MaxValue * NextUInt64()) >> 64);
    }

    /// <summary>
    /// Возвращает число в диапазоне [0,maxValue)
    /// </summary>
    public long NextInt64(long maxValue)
    {
        if (maxValue <= 0)
            throw new ArgumentOutOfRangeException(nameof(maxValue), "Значение должно быть больше 0.");

        return (long)(((BigInteger)maxValue * NextUInt64()) >> 64);
    }

    /// <summary>
    /// Возвращает число в диапазоне [minValue,maxValue)
    /// </summary>
    public long NextInt64(long minValue, long maxValue)
    {
        if (minValue < 0)
            throw new ArgumentOutOfRangeException(nameof(minValue), "Значение должно быть больше либо равно 0.");
        if (maxValue <= minValue)
            throw new ArgumentOutOfRangeException(nameof(maxValue), $"Значение должно быть больше, чем {nameof(minValue)}.");

        return (long)(((BigInteger)(maxValue - minValue) * NextUInt64()) >> 64) + minValue;
    }

    /// <summary>
    /// Возвращает число в диапазоне [0,1)
    /// </summary>
    public double NextDouble()
    {
        return (NextUInt64() >> 11) / 9007199254740992.0; // 2^53
    }

    /// <summary>
    /// Возвращает число в диапазоне [0,1)
    /// </summary>
    public float NextSingle()
    {
        return (NextUInt32() >> 8) / 16777216.0f; // 2^24
    }
}

它很容易使用，如下所示：

static void Main(string[] args)
{
    MersenneTwister mt = new MersenneTwister();
    int randomNumber = mt.Next(10);
    Console.WriteLine(randomNumber);

    Console.ReadKey();
}

性能方面，不错，我就不给数字System.Random了.

它在同一个种子上产生相同的数字序列，就像SFMT 1.5.1在 C 中的原始开发一样，也就是说，这里的计算也一切正常。

从System.Random没有开始被继承，因为不想拉它的内部逻辑。而且由于他自己System.Random没有实现任何东西（他们可能已经为它附加了这么多年的接口），我没有看到联系他的意义。

请看看我在公共方法的转换中是否一切正常，以及我是否对种子太聪明了。或者可能需要调整以使其更快？