//
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
//
// Copyright (c) 2009-2013 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.
//
using System.Collections.Generic;
using System.Runtime.Serialization;
#if !NETSTANDARD1_3
using System;
using System.Runtime;
#endif
namespace IStation.Numerics.Random
{
///
/// Multiplicative congruential generator using a modulus of 2^59 and a multiplier of 13^13.
///
[Serializable]
[DataContract(Namespace = "urn:IStation/Numerics/Random")]
public class Mcg59 : RandomSource
{
const ulong Modulus = 576460752303423488;
const ulong Multiplier = 302875106592253;
const double Reciprocal = 1.0/Modulus;
[DataMember(Order = 1)]
ulong _xn;
///
/// Initializes a new instance of the class using
/// a seed based on time and unique GUIDs.
///
public Mcg59() : this(RandomSeed.Robust())
{
}
///
/// Initializes a new instance of the class using
/// a seed based on time and unique GUIDs.
///
/// if set to true , the class is thread safe.
public Mcg59(bool threadSafe) : this(RandomSeed.Robust(), threadSafe)
{
}
///
/// Initializes a new instance of the class.
///
/// The seed value.
/// If the seed value is zero, it is set to one. Uses the
/// value of to
/// set whether the instance is thread safe.
public Mcg59(int seed)
{
if (seed == 0)
{
seed = 1;
}
_xn = (uint)seed%Modulus;
}
///
/// Initializes a new instance of the class.
///
/// The seed value.
/// The seed is set to 1, if the zero is used as the seed.
/// if set to true , the class is thread safe.
public Mcg59(int seed, bool threadSafe) : base(threadSafe)
{
if (seed == 0)
{
seed = 1;
}
_xn = (uint)seed%Modulus;
}
///
/// Returns a random double-precision floating point number greater than or equal to 0.0, and less than 1.0.
///
protected sealed override double DoSample()
{
double ret = _xn*Reciprocal;
_xn = (_xn*Multiplier)%Modulus;
return ret;
}
///
/// Fills an array with random numbers greater than or equal to 0.0 and less than 1.0.
///
/// Supports being called in parallel from multiple threads.
public static void Doubles(double[] values, int seed)
{
if (seed == 0)
{
seed = 1;
}
ulong xn = (uint)seed%Modulus;
for (int i = 0; i < values.Length; i++)
{
values[i] = xn*Reciprocal;
xn = (xn*Multiplier)%Modulus;
}
}
///
/// Returns an array of random numbers greater than or equal to 0.0 and less than 1.0.
///
/// Supports being called in parallel from multiple threads.
[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;
}
///
/// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0.
///
/// Supports being called in parallel from multiple threads, but the result must be enumerated from a single thread each.
public static IEnumerable DoubleSequence(int seed)
{
if (seed == 0)
{
seed = 1;
}
ulong xn = (uint)seed%Modulus;
while (true)
{
yield return xn*Reciprocal;
xn = (xn*Multiplier)%Modulus;
}
}
}
}