RoApp/Assets/Src/Plugins/MersenneTwister.cs

using System.Collections;
using System;

namespace NPack
{

    /// <summary>

    /// Generates pseudo-random numbers using the Mersenne Twister algorithm.

    /// </summary>

    /// <remarks>

    /// See <a href="http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html">

    /// http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html%3C/a> for details

    /// on the algorithm.

    /// </remarks>

    public class MersenneTwister : Random
    {

        /// <summary>

        /// Creates a new pseudo-random number generator with a given seed.

        /// </summary>

        /// <param name="seed">A value to use as a seed.</param>

        public MersenneTwister(Int32 seed)
        {

            init((UInt32)seed);

        }

        /// <summary>

        /// Creates a new pseudo-random number generator with a default seed.

        /// </summary>

        /// <remarks>

        /// <c>new <see cref="System.Random"/>().<see cref="Random.Next()"/></c>

        /// is used for the seed.

        /// </remarks>

        public MersenneTwister()

            : this(new Random().Next()) /* a default initial seed is used   */

        { }

        /// <summary>

        /// Creates a pseudo-random number generator initialized with the given array.

        /// </summary>

        /// <param name="initKey">The array for initializing keys.</param>

        public MersenneTwister(Int32[] initKey)
        {

            if (initKey == null)
            {

                throw new ArgumentNullException("initKey");

            }

            UInt32[] initArray = new UInt32[initKey.Length];

            for (int i = 0; i < initKey.Length; ++i)
            {

                initArray[i] = (UInt32)initKey[i];

            }

            init(initArray);

        }



        public MersenneTwister(Random rnd)
        {

            // TODO: Complete member initialization

            this.rnd = rnd;

        }

        /// <summary>

        /// Returns the next pseudo-random <see cref="UInt32"/>.

        /// </summary>

        /// <returns>A pseudo-random <see cref="UInt32"/> value.</returns>

        //[CLSCompliant(false)]

        public virtual UInt32 NextUInt32()
        {

            return GenerateUInt32();

        }

        /// <summary>

        /// Returns the next pseudo-random <see cref="UInt32"/>

        /// up to <paramref name="maxValue"/>.

        /// </summary>

        /// <param name="maxValue">

        /// The maximum value of the pseudo-random number to create.

        /// </param>

        /// <returns>

        /// A pseudo-random <see cref="UInt32"/> value which is at most <paramref name="maxValue"/>.

        /// </returns>

        //[CLSCompliant(false)]

        public virtual UInt32 NextUInt32(UInt32 maxValue)
        {

            return (UInt32)(GenerateUInt32() / ((Double)UInt32.MaxValue / maxValue));

        }

        /// <summary>

        /// Returns the next pseudo-random <see cref="UInt32"/> at least

        /// <paramref name="minValue"/> and up to <paramref name="maxValue"/>.

        /// </summary>

        /// <param name="minValue">The minimum value of the pseudo-random number to create.</param>

        /// <param name="maxValue">The maximum value of the pseudo-random number to create.</param>

        /// <returns>

        /// A pseudo-random <see cref="UInt32"/> value which is at least

        /// <paramref name="minValue"/> and at most <paramref name="maxValue"/>.

        /// </returns>

        /// <exception cref="ArgumentOutOfRangeException">

        /// If <c><paramref name="minValue"/> &gt;= <paramref name="maxValue"/></c>.

        /// </exception>

        //[CLSCompliant(false)]

        public virtual UInt32 NextUInt32(UInt32 minValue, UInt32 maxValue) /* throws ArgumentOutOfRangeException */
        {

            if (minValue >= maxValue)
            {

                throw new ArgumentOutOfRangeException();

            }

            return (UInt32)(GenerateUInt32() / ((Double)UInt32.MaxValue / (maxValue - minValue)) + minValue);

        }

        /// <summary>

        /// Returns the next pseudo-random <see cref="Int32"/>.

        /// </summary>

        /// <returns>A pseudo-random <see cref="Int32"/> value.</returns>

        public override Int32 Next()
        {

            return Next(Int32.MaxValue);

        }

        /// <summary>

        /// Returns the next pseudo-random <see cref="Int32"/> up to <paramref name="maxValue"/>.

        /// </summary>

        /// <param name="maxValue">The maximum value of the pseudo-random number to create.</param>

        /// <returns>

        /// A pseudo-random <see cref="Int32"/> value which is at most <paramref name="maxValue"/>.

        /// </returns>

        /// <exception cref="ArgumentOutOfRangeException">

        /// When <paramref name="maxValue"/> &lt; 0.

        /// </exception>

        public override Int32 Next(Int32 maxValue)
        {

            if (maxValue <= 1)
            {

                if (maxValue < 0)
                {

                    throw new ArgumentOutOfRangeException();

                }

                return 0;

            }

            return (Int32)(NextDouble() * maxValue);

        }

        /// <summary>

        /// Returns the next pseudo-random <see cref="Int32"/>

        /// at least <paramref name="minValue"/>

        /// and up to <paramref name="maxValue"/>.

        /// </summary>

        /// <param name="minValue">The minimum value of the pseudo-random number to create.</param>

        /// <param name="maxValue">The maximum value of the pseudo-random number to create.</param>

        /// <returns>A pseudo-random Int32 value which is at least <paramref name="minValue"/> and at

        /// most <paramref name="maxValue"/>.</returns>

        /// <exception cref="ArgumentOutOfRangeException">

        /// If <c><paramref name="minValue"/> &gt;= <paramref name="maxValue"/></c>.

        /// </exception>

        public override Int32 Next(Int32 minValue, Int32 maxValue)
        {

            if (maxValue <= minValue)
            {

                throw new ArgumentOutOfRangeException();

            }

            if (maxValue == minValue)
            {

                return minValue;

            }

            return Next(maxValue - minValue) + minValue;

        }

        /// <summary>

        /// Fills a buffer with pseudo-random bytes.

        /// </summary>

        /// <param name="buffer">The buffer to fill.</param>

        /// <exception cref="ArgumentNullException">

        /// If <c><paramref name="buffer"/> == <see langword="null"/></c>.

        /// </exception>

        public override void NextBytes(Byte[] buffer)
        {

            // [codekaizen: corrected this to check null before checking length.]

            if (buffer == null)
            {

                throw new ArgumentNullException();

            }

            Int32 bufLen = buffer.Length;

            for (Int32 idx = 0; idx < bufLen; ++idx)
            {

                buffer[idx] = (Byte)Next(256);

            }

        }

        /// <summary>

        /// Returns the next pseudo-random <see cref="Double"/> value.

        /// </summary>

        /// <returns>A pseudo-random double floating point value.</returns>

        /// <remarks>

        /// <para>

        /// There are two common ways to create a double floating point using MT19937:

        /// using <see cref="GenerateUInt32"/> and dividing by 0xFFFFFFFF + 1,

        /// or else generating two double words and shifting the first by 26 bits and

        /// adding the second.

        /// </para>

        /// <para>

        /// In a newer measurement of the randomness of MT19937 published in the

        /// journal "Monte Carlo Methods and Applications, Vol. 12, No. 5-6, pp. 385 –§C 393 (2006)"

        /// entitled "A Repetition Test for Pseudo-Random Number Generators",

        /// it was found that the 32-bit version of generating a double fails at the 95%

        /// confidence level when measuring for expected repetitions of a particular

        /// number in a sequence of numbers generated by the algorithm.

        /// </para>

        /// <para>

        /// Due to this, the 53-bit method is implemented here and the 32-bit method

        /// of generating a double is not. If, for some reason,

        /// the 32-bit method is needed, it can be generated by the following:

        /// <code>

        /// (Double)NextUInt32() / ((UInt64)UInt32.MaxValue + 1);

        /// </code>

        /// </para>

        /// </remarks>

        public override Double NextDouble()
        {

            return compute53BitRandom(0, InverseOnePlus53BitsOf1s);

        }

        /// <summary>

        /// Returns a pseudo-random number greater than or equal to zero, and

        /// either strictly less than one, or less than or equal to one,

        /// depending on the value of the given parameter.

        /// </summary>

        /// <param name="includeOne">

        /// If <see langword="true"/>, the pseudo-random number returned will be

        /// less than or equal to one; otherwise, the pseudo-random number returned will

        /// be strictly less than one.

        /// </param>

        /// <returns>

        /// If <paramref name="includeOne"/> is <see langword="true"/>,

        /// this method returns a double-precision pseudo-random number greater than

        /// or equal to zero, and less than or equal to one.

        /// If <paramref name="includeOne"/> is <see langword="false"/>, this method

        /// returns a double-precision pseudo-random number greater than or equal to zero and

        /// strictly less than one.

        /// </returns>

        public Double NextDouble(Boolean includeOne)
        {

            return includeOne ? compute53BitRandom(0, Inverse53BitsOf1s) : NextDouble();

        }

        /// <summary>

        /// Returns a pseudo-random number greater than 0.0 and less than 1.0.

        /// </summary>

        /// <returns>A pseudo-random number greater than 0.0 and less than 1.0.</returns>

        public Double NextDoublePositive()
        {

            return compute53BitRandom(0.5, Inverse53BitsOf1s);

        }

        /// <summary>

        /// Returns a pseudo-random number between 0.0 and 1.0.

        /// </summary>

        /// <returns>

        /// A single-precision floating point number greater than or equal to 0.0,

        /// and less than 1.0.

        /// </returns>

        public Single NextSingle()
        {

            return (Single)NextDouble();

        }

        /// <summary>

        /// Returns a pseudo-random number greater than or equal to zero, and either strictly

        /// less than one, or less than or equal to one, depending on the value of the

        /// given boolean parameter.

        /// </summary>

        /// <param name="includeOne">

        /// If <see langword="true"/>, the pseudo-random number returned will be

        /// less than or equal to one; otherwise, the pseudo-random number returned will

        /// be strictly less than one.

        /// </param>

        /// <returns>

        /// If <paramref name="includeOne"/> is <see langword="true"/>, this method returns a

        /// single-precision pseudo-random number greater than or equal to zero, and less

        /// than or equal to one. If <paramref name="includeOne"/> is <see langword="false"/>,

        /// this method returns a single-precision pseudo-random number greater than or equal to zero and

        /// strictly less than one.

        /// </returns>

        public Single NextSingle(Boolean includeOne)
        {

            return (Single)NextDouble(includeOne);

        }

        /// <summary>

        /// Returns a pseudo-random number greater than 0.0 and less than 1.0.

        /// </summary>

        /// <returns>A pseudo-random number greater than 0.0 and less than 1.0.</returns>

        public Single NextSinglePositive()
        {

            return (Single)NextDoublePositive();

        }

        /// <summary>

        /// Generates a new pseudo-random <see cref="UInt32"/>.

        /// </summary>

        /// <returns>A pseudo-random <see cref="UInt32"/>.</returns>

        //[CLSCompliant(false)]

        protected UInt32 GenerateUInt32()
        {

            UInt32 y;

            /* _mag01[x] = x * MatrixA  for x=0,1 */

            if (_mti >= N) /* generate N words at one time */
            {

                Int16 kk = 0;

                for (; kk < N - M; ++kk)
                {

                    y = (_mt[kk] & UpperMask) | (_mt[kk + 1] & LowerMask);

                    _mt[kk] = _mt[kk + M] ^ (y >> 1) ^ _mag01[y & 0x1];

                }

                for (; kk < N - 1; ++kk)
                {

                    y = (_mt[kk] & UpperMask) | (_mt[kk + 1] & LowerMask);

                    _mt[kk] = _mt[kk + (M - N)] ^ (y >> 1) ^ _mag01[y & 0x1];

                }

                y = (_mt[N - 1] & UpperMask) | (_mt[0] & LowerMask);

                _mt[N - 1] = _mt[M - 1] ^ (y >> 1) ^ _mag01[y & 0x1];

                _mti = 0;

            }

            y = _mt[_mti++];

            y ^= temperingShiftU(y);

            y ^= temperingShiftS(y) & TemperingMaskB;

            y ^= temperingShiftT(y) & TemperingMaskC;

            y ^= temperingShiftL(y);

            return y;

        }



        /* Period parameters */

        private const Int32 N = 624;

        private const Int32 M = 397;

        private const UInt32 MatrixA = 0x9908b0df; /* constant vector a */

        private const UInt32 UpperMask = 0x80000000; /* most significant w-r bits */

        private const UInt32 LowerMask = 0x7fffffff; /* least significant r bits */

        /* Tempering parameters */

        private const UInt32 TemperingMaskB = 0x9d2c5680;

        private const UInt32 TemperingMaskC = 0xefc60000;

        private static UInt32 temperingShiftU(UInt32 y)
        {

            return (y >> 11);

        }

        private static UInt32 temperingShiftS(UInt32 y)
        {

            return (y << 7);

        }

        private static UInt32 temperingShiftT(UInt32 y)
        {

            return (y << 15);

        }

        private static UInt32 temperingShiftL(UInt32 y)
        {

            return (y >> 18);

        }

        private readonly UInt32[] _mt = new UInt32[N]; /* the array for the state vector  */

        private Int16 _mti;

        private static readonly UInt32[] _mag01 = { 0x0, MatrixA };

        private void init(UInt32 seed)
        {

            _mt[0] = seed & 0xffffffffU;

            for (_mti = 1; _mti < N; _mti++)
            {

                _mt[_mti] = (uint)(1812433253U * (_mt[_mti - 1] ^ (_mt[_mti - 1] >> 30)) + _mti);

                // See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.

                // In the previous versions, MSBs of the seed affect  

                // only MSBs of the array _mt[].                       

                // 2002/01/09 modified by Makoto Matsumoto            

                _mt[_mti] &= 0xffffffffU;

                // for >32 bit machines

            }

        }

        private void init(UInt32[] key)
        {

            Int32 i, j, k;

            init(19650218U);

            Int32 keyLength = key.Length;

            i = 1; j = 0;

            k = (N > keyLength ? N : keyLength);

            for (; k > 0; k--)
            {

                _mt[i] = (uint)((_mt[i] ^ ((_mt[i - 1] ^ (_mt[i - 1] >> 30)) * 1664525U)) + key[j] + j); /* non linear */

                _mt[i] &= 0xffffffffU; // for WORDSIZE > 32 machines

                i++; j++;

                if (i >= N) { _mt[0] = _mt[N - 1]; i = 1; }

                if (j >= keyLength) j = 0;

            }

            for (k = N - 1; k > 0; k--)
            {

                _mt[i] = (uint)((_mt[i] ^ ((_mt[i - 1] ^ (_mt[i - 1] >> 30)) * 1566083941U)) - i); /* non linear */

                _mt[i] &= 0xffffffffU; // for WORDSIZE > 32 machines

                i++;

                if (i < N)
                {

                    continue;

                }

                _mt[0] = _mt[N - 1]; i = 1;

            }

            _mt[0] = 0x80000000U; // MSB is 1; assuring non-zero initial array

        }

        // 9007199254740991.0 is the maximum double value which the 53 significand

        // can hold when the exponent is 0.

        private const Double FiftyThreeBitsOf1s = 9007199254740991.0;

        // Multiply by inverse to (vainly?) try to avoid a division.

        private const Double Inverse53BitsOf1s = 1.0 / FiftyThreeBitsOf1s;

        private const Double OnePlus53BitsOf1s = FiftyThreeBitsOf1s + 1;

        private const Double InverseOnePlus53BitsOf1s = 1.0 / OnePlus53BitsOf1s;

        private Random rnd;

        private Double compute53BitRandom(Double translate, Double scale)
        {

            // get 27 pseudo-random bits

            UInt64 a = (UInt64)GenerateUInt32() >> 5;

            // get 26 pseudo-random bits

            UInt64 b = (UInt64)GenerateUInt32() >> 6;

            // shift the 27 pseudo-random bits (a) over by 26 bits (* 67108864.0) and

            // add another pseudo-random 26 bits (+ b).

            return ((a * 67108864.0 + b) + translate) * scale;
        }
    }
}