Seiran128 sample implementation in C# (.NET 5)

Seiran.cs

using System;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Linq;
using System.Runtime.InteropServices;
using System.Security.Cryptography;

// This code is licensed under the terms of the MIT license.
// For details, see https://gist.github.com/andanteyk/d08ab296665b3fc68df58beff3ea39cb .

namespace AndanteSoft
{
    /// <summary>
    /// Seiran128 pseudorandom number generator.
    /// </summary>
    /// <seealso cref="https://github.com/andanteyk/prng-seiran"/>
    public sealed class Seiran
    {
        private ulong State0;
        private ulong State1;



        [ThreadStatic]
        private static Seiran StaticInstance = null;

        /// <summary>
        /// Gets a static instance of <see cref="Seiran"/>.
        /// </summary>
        /// <remarks>
        /// If you don't need reproducibility, we recommend using this.
        /// This instance is thread static.
        /// </remarks>
        public static Seiran I => StaticInstance ??= new Seiran();


        /// <summary>
        /// Initializes a new instance of the <see cref="Seiran"/> class using a cryptographically secure pseudo-random number generator.
        /// </summary>
        /// <remarks>
        /// Basically you should use this constructor.
        /// </remarks>
        public Seiran()
        {
            Span<ulong> state = stackalloc ulong[2];
            var stateBytes = MemoryMarshal.Cast<ulong, byte>(state);

            using (var csprng = new RNGCryptoServiceProvider())
            {
                do
                {
                    csprng.GetBytes(stateBytes);
                } while (state[0] == 0 && state[1] == 0);
            }

            State0 = state[0];
            State1 = state[1];
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="Seiran"/> class using the specified <paramref name="seed"/> value.
        /// </summary>
        /// <param name="seed">
        /// The seed value used to initialize the <see cref="Seiran"/> class.
        /// If the same value is specified, the same sequence of random numbers will be generated.
        /// </param>
        /// <remarks>
        /// This constructor is provided for compatibility with <see cref="System.Random"/>.
        /// You should use <see cref="Seiran()"/> if you want a completely random instance,
        /// and <see cref="Seiran(ReadOnlySpan{ulong})"/> if you need to reproduce a sequence of random numbers.
        /// </remarks>
        public Seiran(long seed)
        {
            ulong x = (ulong)seed;

            static (ulong value, ulong state) splitMix(ulong state)
            {
                ulong z = state += 0x9e3779b97f4a7c15;
                z = (z ^ z >> 30) * 0xbf58476d1ce4e5b9;
                z = (z ^ z >> 27) * 0x94d049bb133111eb;
                return (z ^ z >> 31, state);
            }

            (State0, x) = splitMix(x);
            (State1, _) = splitMix(x);
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="Seiran"/> class using the specified <paramref name="state"/>.
        /// </summary>
        /// <param name="state">
        /// Internal state obtained by <see cref="GetState"/>.
        /// </param>
        /// <exception cref="ArgumentOutOfRangeException">
        /// The length of <paramref name="state"/> is less than 2, or all elements of <paramref name="state"/> are 0.
        /// </exception>
        public Seiran(ReadOnlySpan<ulong> state)
        {
            if (state.Length < 2)
                throw new ArgumentOutOfRangeException(nameof(state));
            if (state[0] == 0 && state[1] == 0)
                throw new ArgumentOutOfRangeException(nameof(state));

            State0 = state[0];
            State1 = state[1];
        }

        /// <summary>
        /// Initialize a new instance of the <see cref="Seiran"/> class that will generate the same sequence of random numbers as the <paramref name="original"/> instance.
        /// </summary>
        /// <param name="original">Copy source instance.</param>
        /// <exception cref="ArgumentNullException"><paramref name="original"/> is null.</exception>
        public Seiran(Seiran original)
        {
            if (original == null)
                throw new ArgumentNullException(nameof(original));

            State0 = original.State0;
            State1 = original.State1;
        }

        /// <summary>
        /// Gets the internal state of this <see cref="Seiran"/> instance.
        /// </summary>
        /// <returns>Internal state. The element length is 2.</returns>
        /// <remarks>
        /// You can use <see cref="Seiran(ReadOnlySpan{ulong})"/> to restore the current state of this instance.
        /// </remarks>
        public ulong[] GetState() => new ulong[] { State0, State1 };



        /// <summary>
        /// Generates a 64-bit random number.
        /// </summary>
        /// <returns>[0, <see cref="ulong.MaxValue"/>]</returns>
        public ulong Next()
        {
            static ulong rotl(ulong x, int k) => x << k | x >> -k;

            ulong s0 = State0, s1 = State1;
            ulong result = rotl((s0 + s1) * 9, 29) + s0;

            State0 = s0 ^ rotl(s1, 29);
            State1 = s0 ^ s1 << 9;

            return result;
        }


        public uint NextUintStrict(uint maxExclusive)
        {
            var hi = Math.BigMul(Next(), maxExclusive, out var lo);
            if (lo < maxExclusive)
            {
                ulong lowerBound = (0ul - maxExclusive) % maxExclusive;

                while (lo < lowerBound)
                {
                    hi = Math.BigMul(Next(), maxExclusive, out lo);
                }
            }
            return (uint)hi;
        }


        /// <summary>
        /// Generates a non-negative random number.
        /// </summary>
        /// <returns>[0, <see cref="int.MaxValue"/>]</returns>
        public int NextInt() => (int)(Next() >> 33);

        /// <summary>
        /// Generates a random number smaller than the specified maximum.
        /// </summary>
        /// <param name="maxExclusive">
        /// An exclusive upper limit for the generated random numbers.
        /// <paramref name="maxExclusive"/> must be greater than 0.
        /// </param>
        /// <returns>[0, <paramref name="maxExclusive"/>)</returns>
        /// <exception cref="ArgumentOutOfRangeException">
        /// <paramref name="maxExclusive"/> is less than or equal to 0.
        /// </exception>
        public int NextInt(int maxExclusive)
        {
            if (maxExclusive <= 0)
                throw new ArgumentOutOfRangeException(nameof(maxExclusive));

            return (int)NextUintStrict((uint)maxExclusive);
        }

        /// <summary>
        /// Generates a random number in the specified range.
        /// </summary>
        /// <param name="minInclusive"> An inclusive lower limit for the generated random numbers. </param>
        /// <param name="maxExclusive"> An exclusive upper limit for the generated random numbers. <paramref name="maxExclusive"/> must be greater than <paramref name="minInclusive"/>. </param>
        /// <returns>[<paramref name="minInclusive"/>, <paramref name="maxExclusive"/>)</returns>
        /// <exception cref="ArgumentOutOfRangeException">
        /// <paramref name="minInclusive"/> is greater than or equal to <paramref name="maxExclusive"/>.
        /// </exception>
        public int NextInt(int minInclusive, int maxExclusive)
        {
            if (minInclusive >= maxExclusive)
                throw new ArgumentOutOfRangeException(nameof(minInclusive));

            return (int)NextUintStrict((uint)(maxExclusive - minInclusive)) + minInclusive;
        }


        /// <summary>
        /// Generates a floating-point random number greater than or equal to 0.0 and less than 1.0.
        /// </summary>
        /// <returns>[0.0, 1.0)</returns>
        /// <remarks>
        /// The actual minimum value greater than 0 is 1.1102230246251565E-16 and the actual maximum value is 0.99999999999999989.
        /// </remarks>
        public double NextDouble() => (Next() >> 11) * (1.0 / (1ul << 53));

        /// <summary>
        /// Generates a floating-point random number in the specified range.
        /// </summary>
        /// <param name="minInclusive"> An inclusive lower limit for the generated random numbers. </param>
        /// <param name="maxExclusive"> An exclusive upper limit for the generated random numbers. <paramref name="maxExclusive"/> must be greater than <paramref name="minInclusive"/>. </param>
        /// <returns>[<paramref name="minInclusive"/>, <paramref name="maxExclusive"/>)</returns>
        /// <exception cref="ArgumentOutOfRangeException">
        /// <paramref name="minInclusive"/> is greater than or equal to <paramref name="maxExclusive"/>.
        /// </exception>
        public double NextDouble(double minInclusive, double maxExclusive)
        {
            if (minInclusive >= maxExclusive)
                throw new ArgumentOutOfRangeException(nameof(minInclusive));

            double t = NextDouble();

            double ret = minInclusive * (1 - t) + maxExclusive * t;

            if (ret == maxExclusive && double.IsFinite(ret))
                ret = Math.BitDecrement(ret);
            return ret;
        }


        /// <summary>
        /// Fills the elements of the specified array with random numbers.
        /// </summary>
        /// <param name="bytes">An array to be filled with random numbers.</param>
        /// <exception cref="ArgumentNullException">
        /// <paramref name="bytes"/> is null.
        /// </exception>
        /// <remarks>
        /// The range of random numbers is [0x00, 0xff].
        /// </remarks>
        public void NextBytes(byte[] bytes)
        {
            if (bytes == null)
                throw new ArgumentNullException(nameof(bytes));

            NextBytes(bytes.AsSpan());
        }

        /// <summary>
        /// Fills the elements of the specified span with random numbers.
        /// </summary>
        /// <param name="buffer">A span to be filled with random numbers.</param>
        /// <remarks>
        /// The range of random numbers is [0x00, 0xff].
        /// If you want to apply to another type, use <see cref="MemoryMarshal.Cast{TFrom, TTo}(Span{TFrom})"/>.
        /// </remarks>
        public void NextBytes(Span<byte> buffer)
        {
            var ulongSpan = MemoryMarshal.Cast<byte, ulong>(buffer);
            for (int i = 0; i < ulongSpan.Length; i++)
                ulongSpan[i] = Next();

            if ((buffer.Length & 7) != 0)
            {
                ulong r = Next();
                for (int i = buffer.Length & ~7; i < buffer.Length; i++)
                {
                    buffer[i] = (byte)r;
                    r >>= 8;
                }
            }
        }


        /// <summary>
        /// Generates independent pairs of random numbers that follow a normal distribution.
        /// </summary>
        /// <returns>Independent pairs of random numbers that follow a normal distribution with mean 0.0 and standard deviation 1.0.</returns>
        public (double x, double y) NextGaussian()
        {
            double u, v, s;
            do
            {
                u = NextDouble(-1, 1);
                v = NextDouble(-1, 1);
                s = u * u + v * v;
            } while (s >= 1 || s == 0);

            s = Math.Sqrt(-2 * Math.Log(s) / s);
            return (u * s, v * s);
        }



        /// <summary>
        /// Shuffles the order of the elements.
        /// </summary>
        /// <typeparam name="T">The type of the elements of <paramref name="list"/>.</typeparam>
        /// <param name="list">List to be shuffled.</param>
        /// <exception cref="ArgumentNullException"> <paramref name="list"/> is null. </exception>
        public void ShuffleList<T>(IList<T> list)
        {
            if (list == null)
                throw new ArgumentNullException(nameof(list));

            for (int i = list.Count - 1; i > 0; i--)
            {
                int r = NextInt(i + 1);
                (list[i], list[r]) = (list[r], list[i]);
            }
        }

        /// <summary>
        /// Shuffles the order of the elements.
        /// </summary>
        /// <typeparam name="T">The type of the elements of <paramref name="span"/>.</typeparam>
        /// <param name="span">List to be shuffled.</param>
        public void ShuffleSpan<T>(Span<T> span)
        {
            for (int i = span.Length - 1; i > 0; i--)
            {
                int r = NextInt(i + 1);
                (span[i], span[r]) = (span[r], span[i]);
            }
        }

        /// <summary>
        /// Shuffles the order of the elements.
        /// </summary>
        /// <typeparam name="T">The type of the elements of <paramref name="source"/>.</typeparam>
        /// <param name="source">Sequence to be shuffled.</param>
        /// <returns>A sequence in which the order of the elements is randomly changed.</returns>
        /// <exception cref="ArgumentNullException"> <paramref name="source"/> is null. </exception>
        /// <remarks>
        /// This method is implemented using deferred execution.
        /// The <paramref name="source"/> is not evaluated and shuffled until GetEnumerator() is called.
        /// </remarks>
        public IEnumerable<T> Shuffle<T>(IEnumerable<T> source)
        {
            if (source == null)
                throw new ArgumentNullException(nameof(source));

            return new ShuffleIterator<T>(this, source);
        }

        public sealed class ShuffleIterator<T> : IEnumerable<T>, IEnumerator<T>
        {
            private readonly IEnumerable<T> Source;
            private readonly Seiran RandomSource;
            private int CurrentIndex;
            private T[] Elements;

            internal ShuffleIterator(Seiran randomSource, IEnumerable<T> source)
            {
                RandomSource = randomSource ?? throw new ArgumentNullException(nameof(randomSource));
                Source = source ?? throw new ArgumentNullException(nameof(source));
                CurrentIndex = -1;
            }

            public T Current
            {
                get
                {
                    if ((uint)CurrentIndex < Elements?.Length)
                        return Elements[CurrentIndex];
                    else
                        return default;
                }
            }
            object IEnumerator.Current => Current;

            public void Dispose() { }

            public bool MoveNext()
            {
                if (CurrentIndex >= Elements.Length)
                    return false;

                return ++CurrentIndex < Elements.Length;
            }

            public void Reset() => throw new NotSupportedException();

            public IEnumerator<T> GetEnumerator()
            {
                if (Elements == null)
                {
                    Elements = Source.ToArray();
                    RandomSource.ShuffleSpan(Elements.AsSpan());
                }

                if (CurrentIndex == -1)
                    return this;

                var clone = MemberwiseClone() as ShuffleIterator<T>;
                clone.CurrentIndex = -1;
                return clone;
            }
            IEnumerator IEnumerable.GetEnumerator() => GetEnumerator();
        }



        private void Jump(ReadOnlySpan<ulong> jumpPolynomial)
        {
            ulong t0 = 0, t1 = 0;

            for (int i = 0; i < 2; i++)
            {
                for (int b = 0; b < 64; b++)
                {
                    if (((jumpPolynomial[i] >> b) & 1) != 0)
                    {
                        t0 ^= State0;
                        t1 ^= State1;
                    }
                    Next();
                }
            }

            State0 = t0;
            State1 = t1;
        }

        /// <summary>
        /// Advance the internal state by 2^32 steps.
        /// </summary>
        /// <remarks>
        /// This method is equivalent to 2^32 <see cref="Next"/> calls.
        /// It can be executed in the same amount of time as 128 <see cref="Next"/> calls.
        /// </remarks>
        public void Jump32() => Jump(stackalloc ulong[] { 0x40165CBAE9CA6DEB, 0x688E6BFC19485AB1 });

        /// <summary>
        /// Advance the internal state by 2^64 steps.
        /// </summary>
        /// <remarks>
        /// This method is equivalent to 2^64 <see cref="Next"/> calls.
        /// It can be executed in the same amount of time as 128 <see cref="Next"/> calls.
        /// </remarks>
        public void Jump64() => Jump(stackalloc ulong[] { 0xF4DF34E424CA5C56, 0x2FE2DE5C2E12F601 });

        /// <summary>
        /// Advance the internal state by 2^96 steps.
        /// </summary>
        /// <remarks>
        /// This method is equivalent to 2^96 <see cref="Next"/> calls.
        /// It can be executed in the same amount of time as 128 <see cref="Next"/> calls.
        /// </remarks>
        public void Jump96() => Jump(stackalloc ulong[] { 0x185F4DF8B7634607, 0x95A98C7025F908B2 });


        /// <summary>
        /// Rewinds the internal state to the previous state.
        /// </summary>
        public void Prev()
        {
            static ulong rotl(ulong x, int k) => x << k | x >> -k;

            ulong t1 = rotl(State0 ^ State1, 64 - 29);
            t1 ^= t1 << 44 ^ (t1 & ~0xffffful) << 24;
            t1 ^= (t1 & (0x7ffful << 40)) << 4;
            t1 ^= (t1 & (0x7fful << 40)) << 8;
            t1 ^= (t1 & (0x7ul << 40)) << 16;
            t1 ^= (t1 & (0xffffful << 35)) >> 20;
            t1 ^= (t1 & (0x7ffful << 20)) >> 20;

            State0 ^= rotl(t1, 29);
            State1 = t1;
        }
    }
}