// <copyright file="MersenneTwister.cs" company="Math.NET">
|
// Math.NET Numerics, part of the Math.NET Project
|
// http://numerics.mathdotnet.com
|
// http://github.com/mathnet/mathnet-numerics
|
//
|
// Copyright (c) 2009-2017 Math.NET
|
//
|
// Permission is hereby granted, free of charge, to any person
|
// obtaining a copy of this software and associated documentation
|
// files (the "Software"), to deal in the Software without
|
// restriction, including without limitation the rights to use,
|
// copy, modify, merge, publish, distribute, sublicense, and/or sell
|
// copies of the Software, and to permit persons to whom the
|
// Software is furnished to do so, subject to the following
|
// conditions:
|
//
|
// The above copyright notice and this permission notice shall be
|
// included in all copies or substantial portions of the Software.
|
//
|
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
|
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
|
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
|
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
|
// OTHER DEALINGS IN THE SOFTWARE.
|
// </copyright>
|
|
/*
|
Original code's copyright and license:
|
Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
|
All rights reserved.
|
|
Redistribution and use in source and binary forms, with or without
|
modification, are permitted provided that the following conditions
|
are met:
|
|
1. Redistributions of source code must retain the above copyright
|
notice, this list of conditions and the following disclaimer.
|
|
2. Redistributions in binary form must reproduce the above copyright
|
notice, this list of conditions and the following disclaimer in the
|
documentation and/or other materials provided with the distribution.
|
|
3. The names of its contributors may not be used to endorse or promote
|
products derived from this software without specific prior written
|
permission.
|
|
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
|
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
|
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
|
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
|
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
|
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
|
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
|
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
|
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
|
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
|
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
|
|
Any feedback is very welcome.
|
http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
|
email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
|
*/
|
|
using System.Collections.Generic;
|
using System.Runtime.Serialization;
|
using System.Threading;
|
|
#if !NETSTANDARD1_3
|
using System;
|
using System.Runtime;
|
#endif
|
|
namespace IStation.Numerics.Random
|
{
|
/// <summary>
|
/// Random number generator using Mersenne Twister 19937 algorithm.
|
/// </summary>
|
[Serializable]
|
[DataContract(Namespace = "urn:IStation/Numerics/Random")]
|
public class MersenneTwister : RandomSource
|
{
|
/// <summary>
|
/// Mersenne twister constant.
|
/// </summary>
|
const uint LowerMask = 0x7fffffff;
|
|
/// <summary>
|
/// Mersenne twister constant.
|
/// </summary>
|
const int M = 397;
|
|
/// <summary>
|
/// Mersenne twister constant.
|
/// </summary>
|
const uint MatrixA = 0x9908b0df;
|
|
/// <summary>
|
/// Mersenne twister constant.
|
/// </summary>
|
const int N = 624;
|
|
/// <summary>
|
/// Mersenne twister constant.
|
/// </summary>
|
const double Reciprocal = 1.0/4294967296.0; // 1.0/(uint.MaxValue + 1.0)
|
|
/// <summary>
|
/// Mersenne twister constant.
|
/// </summary>
|
const uint UpperMask = 0x80000000;
|
|
/// <summary>
|
/// Mersenne twister constant.
|
/// </summary>
|
static readonly uint[] Mag01 = { 0x0U, MatrixA };
|
|
/// <summary>
|
/// Mersenne twister constant.
|
/// </summary>
|
[DataMember(Order = 1)]
|
readonly uint[] _mt = new uint[N];
|
|
/// <summary>
|
/// Mersenne twister constant.
|
/// </summary>
|
[DataMember(Order = 2)]
|
int _mti = N + 1;
|
|
/// <summary>
|
/// Initializes a new instance of the <see cref="MersenneTwister"/> class using
|
/// a seed based on time and unique GUIDs.
|
/// </summary>
|
/// <remarks>If the seed value is zero, it is set to one. Uses the
|
/// value of <see cref="Control.ThreadSafeRandomNumberGenerators"/> to
|
/// set whether the instance is thread safe.</remarks>
|
public MersenneTwister() : this(RandomSeed.Robust())
|
{
|
}
|
|
/// <summary>
|
/// Initializes a new instance of the <see cref="MersenneTwister"/> class using
|
/// a seed based on time and unique GUIDs.
|
/// </summary>
|
/// <param name="threadSafe">if set to <c>true</c> , the class is thread safe.</param>
|
public MersenneTwister(bool threadSafe) : this(RandomSeed.Robust(), threadSafe)
|
{
|
}
|
|
/// <summary>
|
/// Initializes a new instance of the <see cref="MersenneTwister"/> class.
|
/// </summary>
|
/// <param name="seed">The seed value.</param>
|
/// <remarks>Uses the value of <see cref="Control.ThreadSafeRandomNumberGenerators"/> to
|
/// set whether the instance is thread safe.</remarks>
|
public MersenneTwister(int seed)
|
{
|
init_genrand((uint)seed);
|
}
|
|
/// <summary>
|
/// Initializes a new instance of the <see cref="MersenneTwister"/> class.
|
/// </summary>
|
/// <param name="seed">The seed value.</param>
|
/// <param name="threadSafe">if set to <c>true</c>, the class is thread safe.</param>
|
public MersenneTwister(int seed, bool threadSafe) : base(threadSafe)
|
{
|
init_genrand((uint)seed);
|
}
|
|
static readonly ThreadLocal<MersenneTwister> DefaultInstance = new ThreadLocal<MersenneTwister>(() => new MersenneTwister(RandomSeed.Robust(), true));
|
|
/// <summary>
|
/// Default instance, thread-safe.
|
/// </summary>
|
public static MersenneTwister Default => DefaultInstance.Value;
|
|
/*/// <summary>
|
/// Initializes a new instance of the <see cref="MersenneTwister"/> class.
|
/// </summary>
|
/// <param name="init_key">The initialization key.</param>
|
public MersenneTwister(int[] init_key)
|
{
|
if (init_key == null)
|
{
|
throw new ArgumentNullException("init_key");
|
}
|
uint[] array = new uint[init_key.Length];
|
for (int i = 0; i < array.Length; i++)
|
{
|
array[i] = (uint) init_key[i];
|
}
|
init_by_array(array);
|
}
|
*/
|
/* initializes _mt[_n] with a seed */
|
|
void init_genrand(uint s)
|
{
|
_mt[0] = s & 0xffffffff;
|
for (_mti = 1; _mti < N; _mti++)
|
{
|
_mt[_mti] = 1812433253*(_mt[_mti - 1] ^ (_mt[_mti - 1] >> 30)) + (uint)_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] &= 0xffffffff;
|
/* for >32 bit machines */
|
}
|
}
|
|
/* initialize by an array with array-length */
|
/* init_key is the array for initializing keys */
|
/* slight change for C++, 2004/2/26 */
|
|
/* private void init_by_array(uint[] init_key)
|
{
|
uint key_length = (uint) init_key.Length;
|
init_genrand(19650218);
|
uint i = 1;
|
uint j = 0;
|
uint k = (_n > key_length ? _n : key_length);
|
for (; k > 0; k--)
|
{
|
_mt[i] = (_mt[i] ^ ((_mt[i - 1] ^ (_mt[i - 1] >> 30))*1664525)) + init_key[j] + j; //non linear
|
_mt[i] &= 0xffffffff; // for WORDSIZE > 32 machines
|
i++;
|
j++;
|
if (i >= _n)
|
{
|
_mt[0] = _mt[_n - 1];
|
i = 1;
|
}
|
if (j >= key_length) j = 0;
|
}
|
for (k = _n - 1; k > 0; k--)
|
{
|
_mt[i] = (_mt[i] ^ ((_mt[i - 1] ^ (_mt[i - 1] >> 30))*1566083941)) - i; // non linear
|
_mt[i] &= 0xffffffff; // for WORDSIZE > 32 machines
|
i++;
|
if (i >= _n)
|
{
|
_mt[0] = _mt[_n - 1];
|
i = 1;
|
}
|
}
|
|
_mt[0] = 0x80000000; // MSB is 1; assuring non-zero initial array
|
}*/
|
|
/* generates a random number on [0,0xffffffff]-interval */
|
|
uint genrand_int32()
|
{
|
uint y;
|
|
/* mag01[x] = x * MATRIX_A for x=0,1 */
|
|
if (_mti >= N)
|
{
|
/* generate _n words at one time */
|
int kk;
|
|
if (_mti == N + 1) /* if init_genrand() has not been called, */
|
{
|
init_genrand(5489); /* a default initial seed is used */
|
}
|
|
for (kk = 0; 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++];
|
|
/* Tempering */
|
y ^= y >> 11;
|
y ^= (y << 7) & 0x9d2c5680;
|
y ^= (y << 15) & 0xefc60000;
|
y ^= y >> 18;
|
|
return y;
|
}
|
|
/// <summary>
|
/// Returns a random double-precision floating point number greater than or equal to 0.0, and less than 1.0.
|
/// </summary>
|
protected sealed override double DoSample()
|
{
|
return genrand_int32()*Reciprocal;
|
}
|
|
/// <summary>
|
/// Returns a random 32-bit signed integer greater than or equal to zero and less than <see cref="F:System.Int32.MaxValue"/>
|
/// </summary>
|
protected sealed override int DoSampleInteger()
|
{
|
uint uint32 = genrand_int32();
|
int int31 = (int)(uint32 >> 1);
|
if (int31 == int.MaxValue)
|
{
|
return DoSampleInteger();
|
}
|
|
return int31;
|
}
|
|
/// <summary>
|
/// Fills the elements of a specified array of bytes with random numbers in full range, including zero and 255 (<see cref="F:System.Byte.MaxValue"/>).
|
/// </summary>
|
protected sealed override void DoSampleBytes(byte[] buffer)
|
{
|
for (var i = 0; i < buffer.Length; i++)
|
{
|
buffer[i] = (byte)(genrand_int32() % 256);
|
}
|
}
|
|
/* /// <summary>
|
/// Generates a random number on [0,1) with 53-bit resolution.
|
/// </summary>
|
/// <returns>A random number on [0,1) with 53-bit resolution.</returns>
|
public double NextDoubleResolution53()
|
{
|
ulong a = genrand_int32() >> 5, b = genrand_int32() >> 6;
|
return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0);
|
}*/
|
|
/// <summary>
|
/// Fills an array with random numbers greater than or equal to 0.0 and less than 1.0.
|
/// </summary>
|
/// <remarks>Supports being called in parallel from multiple threads.</remarks>
|
public static void Doubles(double[] values, int seed)
|
{
|
uint[] t = new uint[624];
|
int k;
|
uint s = (uint)seed;
|
|
t[0] = s & 0xffffffff;
|
for (k = 1; k < N; k++)
|
{
|
t[k] = 1812433253*(t[k - 1] ^ (t[k - 1] >> 30)) + (uint)k;
|
t[k] &= 0xffffffff;
|
}
|
|
for (int i = 0; i < values.Length; i++)
|
{
|
uint y;
|
|
if (k >= N)
|
{
|
int kk;
|
for (kk = 0; kk < N - M; kk++)
|
{
|
y = (t[kk] & UpperMask) | (t[kk + 1] & LowerMask);
|
t[kk] = t[kk + M] ^ (y >> 1) ^ Mag01[y & 0x1];
|
}
|
|
for (; kk < N - 1; kk++)
|
{
|
y = (t[kk] & UpperMask) | (t[kk + 1] & LowerMask);
|
t[kk] = t[kk + (M - N)] ^ (y >> 1) ^ Mag01[y & 0x1];
|
}
|
|
y = (t[N - 1] & UpperMask) | (t[0] & LowerMask);
|
t[N - 1] = t[M - 1] ^ (y >> 1) ^ Mag01[y & 0x1];
|
|
k = 0;
|
}
|
|
y = t[k++];
|
|
/* Tempering */
|
y ^= y >> 11;
|
y ^= (y << 7) & 0x9d2c5680;
|
y ^= (y << 15) & 0xefc60000;
|
y ^= y >> 18;
|
|
values[i] = y*Reciprocal;
|
}
|
}
|
|
/// <summary>
|
/// Returns an array of random numbers greater than or equal to 0.0 and less than 1.0.
|
/// </summary>
|
/// <remarks>Supports being called in parallel from multiple threads.</remarks>
|
[TargetedPatchingOptOut("Performance critical to inline this type of method across NGen image boundaries")]
|
public static double[] Doubles(int length, int seed)
|
{
|
var data = new double[length];
|
Doubles(data, seed);
|
return data;
|
}
|
|
/// <summary>
|
/// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0.
|
/// </summary>
|
/// <remarks>Supports being called in parallel from multiple threads, but the result must be enumerated from a single thread each.</remarks>
|
public static IEnumerable<double> DoubleSequence(int seed)
|
{
|
uint[] t = new uint[624];
|
int k;
|
uint s = (uint)seed;
|
|
t[0] = s & 0xffffffff;
|
for (k = 1; k < N; k++)
|
{
|
t[k] = 1812433253*(t[k - 1] ^ (t[k - 1] >> 30)) + (uint)k;
|
t[k] &= 0xffffffff;
|
}
|
|
while (true)
|
{
|
uint y;
|
|
if (k >= N)
|
{
|
int kk;
|
for (kk = 0; kk < N - M; kk++)
|
{
|
y = (t[kk] & UpperMask) | (t[kk + 1] & LowerMask);
|
t[kk] = t[kk + M] ^ (y >> 1) ^ Mag01[y & 0x1];
|
}
|
|
for (; kk < N - 1; kk++)
|
{
|
y = (t[kk] & UpperMask) | (t[kk + 1] & LowerMask);
|
t[kk] = t[kk + (M - N)] ^ (y >> 1) ^ Mag01[y & 0x1];
|
}
|
|
y = (t[N - 1] & UpperMask) | (t[0] & LowerMask);
|
t[N - 1] = t[M - 1] ^ (y >> 1) ^ Mag01[y & 0x1];
|
|
k = 0;
|
}
|
|
y = t[k++];
|
|
/* Tempering */
|
y ^= y >> 11;
|
y ^= (y << 7) & 0x9d2c5680;
|
y ^= (y << 15) & 0xefc60000;
|
y ^= y >> 18;
|
|
yield return y*Reciprocal;
|
}
|
}
|
}
|
}
|
