// <copyright file="Xoshiro256StarStar.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-2014 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:
|
Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org)
|
|
To the extent possible under law, the author has dedicated all copyright
|
and related and neighboring rights to this software to the public domain
|
worldwide. This software is distributed without any warranty.
|
|
See <http://creativecommons.org/publicdomain/zero/1.0/>.
|
*/
|
|
using System.Collections.Generic;
|
using System.Runtime.Serialization;
|
|
#if !NETSTANDARD1_3
|
using System;
|
using System.Runtime;
|
#endif
|
|
namespace IStation.Numerics.Random
|
{
|
/// <summary>
|
/// Xoshiro256** pseudo random number generator.
|
/// A random number generator based on the <see cref="System.Random"/> class in the .NET library.
|
/// </summary>
|
/// <remarks>
|
/// This is xoshiro256** 1.0, our all-purpose, rock-solid generator. It has
|
/// excellent(sub-ns) speed, a state space(256 bits) that is large enough
|
/// for any parallel application, and it passes all tests we are aware of.
|
///
|
/// For generating just floating-point numbers, xoshiro256+ is even faster.
|
///
|
/// The state must be seeded so that it is not everywhere zero.If you have
|
/// a 64-bit seed, we suggest to seed a splitmix64 generator and use its
|
/// output to fill s.
|
///
|
/// For further details see:
|
/// David Blackman & Sebastiano Vigna (2018), "Scrambled Linear Pseudorandom Number Generators".
|
/// https://arxiv.org/abs/1805.01407
|
/// </remarks>
|
[Serializable]
|
[DataContract(Namespace = "urn:IStation/Numerics/Random")]
|
public class Xoshiro256StarStar : RandomSource
|
{
|
// Constants.
|
const double REAL_UNIT_UINT = 1.0 / (1UL << 53);
|
|
// RNG state.
|
[DataMember(Order = 1)]
|
ulong _s0;
|
[DataMember(Order = 2)]
|
ulong _s1;
|
[DataMember(Order = 3)]
|
ulong _s2;
|
[DataMember(Order = 4)]
|
ulong _s3;
|
|
/// <summary>
|
/// Construct a new random number generator with a random seed.
|
/// </summary>
|
public Xoshiro256StarStar() : this(RandomSeed.Robust())
|
{
|
}
|
|
/// <summary>
|
/// Construct a new random number generator with random seed.
|
/// </summary>
|
/// <param name="threadSafe">if set to <c>true</c> , the class is thread safe.</param>
|
public Xoshiro256StarStar(bool threadSafe) : this(RandomSeed.Robust(), threadSafe)
|
{
|
}
|
|
/// <summary>
|
/// Construct a new random number generator with random seed.
|
/// </summary>
|
/// <param name="seed">The seed value.</param>
|
public Xoshiro256StarStar(int seed)
|
{
|
Initialise(seed);
|
}
|
|
/// <summary>
|
/// Construct a new random number generator with random seed.
|
/// </summary>
|
/// <param name="seed">The seed value.</param>
|
/// <param name="threadSafe">if set to <c>true</c> , the class is thread safe.</param>
|
public Xoshiro256StarStar(int seed, bool threadSafe) : base(threadSafe)
|
{
|
Initialise(seed);
|
}
|
|
/// <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()
|
{
|
// Note. Here we generate a random integer between 0 and 2^53-1 (i.e. 53 binary 1s) and multiply
|
// by the fractional unit value 1.0 / 2^53, thus the result has a max value of
|
// 1.0 - (1.0 / 2^53), or 0.99999999999999989 in decimal.
|
return (NextInnerULong() >> 11) * REAL_UNIT_UINT;
|
}
|
|
/// <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 override int DoSampleInteger()
|
{
|
retry:
|
// Handle the special case where the value int.MaxValue is generated; this is outside
|
// the range of permitted return values for this method.
|
ulong rtn = NextInnerULong() & 0x7fff_ffffUL;
|
if (rtn == 0x7fff_ffffUL)
|
{
|
goto retry;
|
}
|
return (int)rtn;
|
}
|
|
/// <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 override void DoSampleBytes(byte[] buffer)
|
{
|
// For improved performance the below loop operates on these stack allocated copies of the heap variables.
|
// Notes. doing this means that these heavily used variables are located near to other local/stack variables,
|
// thus they will very likely be cached in the same CPU cache line.
|
ulong s0 = _s0;
|
ulong s1 = _s1;
|
ulong s2 = _s2;
|
ulong s3 = _s3;
|
|
int i = 0;
|
|
// Fill up the bulk of the buffer in chunks of 8 bytes at a time.
|
for (int bound = buffer.Length - 3; i < bound;)
|
{
|
// Generate 64 random bits.
|
ulong x = RotateLeft(s1 * 5, 7) * 9;
|
|
// Update PRNG state.
|
ulong t = s1 << 17;
|
s2 ^= s0;
|
s3 ^= s1;
|
s1 ^= s2;
|
s0 ^= s3;
|
s2 ^= t;
|
s3 = RotateLeft(s3, 45);
|
|
// Assign bits to a segment of 8 bytes.
|
buffer[i++] = (byte)x;
|
buffer[i++] = (byte)(x >> 8);
|
buffer[i++] = (byte)(x >> 16);
|
buffer[i++] = (byte)(x >> 24);
|
buffer[i++] = (byte)(x >> 32);
|
buffer[i++] = (byte)(x >> 40);
|
buffer[i++] = (byte)(x >> 48);
|
buffer[i++] = (byte)(x >> 56);
|
}
|
|
// Fill up any remaining bytes in the buffer.
|
if (i < buffer.Length)
|
{
|
// Generate 64 random bits.
|
ulong x = RotateLeft(s1 * 5, 7) * 9;
|
|
// Update PRNG state.
|
ulong t = s1 << 17;
|
s2 ^= s0;
|
s3 ^= s1;
|
s1 ^= s2;
|
s0 ^= s3;
|
s2 ^= t;
|
s3 = RotateLeft(s3, 45);
|
|
// Allocate one byte at a time until we reach the end of the buffer.
|
while (i < buffer.Length)
|
{
|
buffer[i++] = (byte)x;
|
x >>= 8;
|
}
|
}
|
|
// Update the state variables on the heap.
|
_s0 = s0;
|
_s1 = s1;
|
_s2 = s2;
|
_s3 = s3;
|
}
|
|
/// <summary>
|
/// Returns a random N-bit signed integer greater than or equal to zero and less than 2^N.
|
/// N (bit count) is expected to be greater than zero and less than 32 (not verified).
|
/// </summary>
|
protected override int DoSampleInt32WithNBits(int bitCount)
|
{
|
return (int)(NextInnerULong() >> (64 - bitCount));
|
}
|
|
/// <summary>
|
/// Returns a random N-bit signed long integer greater than or equal to zero and less than 2^N.
|
/// N (bit count) is expected to be greater than zero and less than 64 (not verified).
|
/// </summary>
|
protected override long DoSampleInt64WithNBits(int bitCount)
|
{
|
return (long)(NextInnerULong() >> (64 - bitCount));
|
}
|
|
private void Initialise(int seed)
|
{
|
// Notes.
|
// xoroshiro256** requires that at least one of the state variable be non-zero, use of splitmix64
|
// satisfies that requirement because its outputs are equidistributed, i.e. if a zero is output
|
// then the next zero will be after a further 2^64 outputs.
|
// Splitmix will also accept a zero input, thus all possible seeds can be accepted here and will
|
// all generate good initial state for xoshiro256**.
|
ulong longSeed = (ulong)seed;
|
|
// Use the splitmix64 RNG to hash the seed.
|
_s0 = Splitmix64(ref longSeed);
|
_s1 = Splitmix64(ref longSeed);
|
_s2 = Splitmix64(ref longSeed);
|
_s3 = Splitmix64(ref longSeed);
|
}
|
|
private ulong NextInnerULong()
|
{
|
ulong s0 = _s0;
|
ulong s1 = _s1;
|
ulong s2 = _s2;
|
ulong s3 = _s3;
|
|
ulong result = RotateLeft(s1 * 5, 7) * 9;
|
|
ulong t = s1 << 17;
|
s2 ^= s0;
|
s3 ^= s1;
|
s1 ^= s2;
|
s0 ^= s3;
|
s2 ^= t;
|
s3 = RotateLeft(s3, 45);
|
|
_s0 = s0;
|
_s1 = s1;
|
_s2 = s2;
|
_s3 = s3;
|
|
return result;
|
}
|
|
/// <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)
|
{
|
// Init state.
|
ulong longSeed = (ulong)seed;
|
ulong s0 = Splitmix64(ref longSeed);
|
ulong s1 = Splitmix64(ref longSeed);
|
ulong s2 = Splitmix64(ref longSeed);
|
ulong s3 = Splitmix64(ref longSeed);
|
|
for (int i = 0; i < values.Length; i++)
|
{
|
// Generate sample.
|
values[i] = ((RotateLeft(s1 * 5, 7) * 9) >> 11) * REAL_UNIT_UINT;
|
|
// Update PRNG state.
|
ulong t = s1 << 17;
|
s2 ^= s0;
|
s3 ^= s1;
|
s1 ^= s2;
|
s0 ^= s3;
|
s2 ^= t;
|
s3 = RotateLeft(s3, 45);
|
}
|
}
|
|
/// <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)
|
{
|
// Init state.
|
ulong longSeed = (ulong)seed;
|
ulong s0 = Splitmix64(ref longSeed);
|
ulong s1 = Splitmix64(ref longSeed);
|
ulong s2 = Splitmix64(ref longSeed);
|
ulong s3 = Splitmix64(ref longSeed);
|
|
while (true)
|
{
|
// Generate sample.
|
double x = ((RotateLeft(s1 * 5, 7) * 9) >> 11) * REAL_UNIT_UINT;
|
|
// Update PRNG state.
|
ulong t = s1 << 17;
|
s2 ^= s0;
|
s3 ^= s1;
|
s1 ^= s2;
|
s0 ^= s3;
|
s2 ^= t;
|
s3 = RotateLeft(s3, 45);
|
|
// Yield sample.
|
yield return x;
|
}
|
}
|
|
/// <summary>
|
/// Splitmix64 RNG.
|
/// </summary>
|
/// <param name="x">RNG state. This can take any value, including zero.</param>
|
/// <returns>A new random UInt64.</returns>
|
/// <remarks>
|
/// Splitmix64 produces equidistributed outputs, thus if a zero is generated then the
|
/// next zero will be after a further 2^64 outputs.
|
/// </remarks>
|
private static ulong Splitmix64(ref ulong x)
|
{
|
ulong z = (x += 0x9E3779B97F4A7C15UL);
|
z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9UL;
|
z = (z ^ (z >> 27)) * 0x94D049BB133111EBUL;
|
return z ^ (z >> 31);
|
}
|
|
private static ulong RotateLeft(ulong x, int k)
|
{
|
// Note. RyuJIT will compile this to a single rotate CPU instruction (as of about .NET 4.6.1 and dotnet core 2.0).
|
return (x << k) | (x >> (64 - k));
|
}
|
}
|
}
|
