Dotnet/FastRng
Dotnet/FastRng
1
0
Fork 0
Code Issues 2 Pull Requests Actions Packages Projects Releases 5 Wiki Activity

Added float implementations

Browse Source
This commit is contained in:
Thorsten Sommer 2020-10-31 23:27:18 +01:00
parent 01ee5900d5
commit c340635c4a
6 changed files with 387 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
FastRng/Float/Distributions/IDistribution.cs Normal file
View File

@ -0,0 +1,12 @@
using System.Threading;
using System.Threading.Tasks;
namespace FastRng.Float.Distributions
{
public interface IDistribution
{
public IRandom Random { get; set; }
public ValueTask<float> GetDistributedValue(CancellationToken token);
}
}

12
FastRng/Float/Distributions/Uniform.cs Normal file
View File

@ -0,0 +1,12 @@
using System.Threading;
using System.Threading.Tasks;
namespace FastRng.Float.Distributions
{
public sealed class Uniform : IDistribution
{
public IRandom Random { get; set; }
public async ValueTask<float> GetDistributedValue(CancellationToken token = default) => this.Random == null ? float.NaN : await this.Random.GetUniform(token);
}
}

21
FastRng/Float/IRandom.cs Normal file
View File

@ -0,0 +1,21 @@
using System.Threading;
using System.Threading.Tasks;
using FastRng.Float.Distributions;
namespace FastRng.Float
{
public interface IRandom
{
public ValueTask<float> GetUniform(CancellationToken cancel = default);
public ValueTask<uint> NextNumber(uint rangeStart, uint rangeEnd, IDistribution distribution, CancellationToken cancel = default);
public ValueTask<ulong> NextNumber(ulong rangeStart, ulong rangeEnd, IDistribution distribution, CancellationToken cancel = default);
public ValueTask<float> NextNumber(float rangeStart, float rangeEnd, IDistribution distribution, CancellationToken cancel = default);
public ValueTask<float> NextNumber(IDistribution distribution, CancellationToken cancel = default);
public void StopProducer();
}
}

61
FastRng/Float/MathTools.cs Normal file
View File

@ -0,0 +1,61 @@
using System;
namespace FastRng.Float
{
public static class MathTools
{
private static readonly float SQRT_2 = MathF.Sqrt(2.0f);
private static readonly float SQRT_PI = MathF.Sqrt(MathF.PI);
public static float Gamma(float z)
{
// Source: http://rosettacode.org/wiki/Gamma_function#Go
const float F1 = 6.5f;
const float A1 = .99999999999980993f;
const float A2 = 676.5203681218851f;
const float A3 = 1259.1392167224028f;
const float A4 = 771.32342877765313f;
const float A5 = 176.61502916214059f;
const float A6 = 12.507343278686905f;
const float A7 = .13857109526572012f;
const float A8 = 9.9843695780195716e-6f;
const float A9 = 1.5056327351493116e-7f;
var t = z + F1;
var x = A1 +
A2 / z -
A3 / (z + 1) +
A4 / (z + 2) -
A5 / (z + 3) +
A6 / (z + 4) -
A7 / (z + 5) +
A8 / (z + 6) +
A9 / (z + 7);
return MathTools.SQRT_2 * MathTools.SQRT_PI * MathF.Pow(t, z - 0.5f) * MathF.Exp(-t) * x;
}
public static float Factorial(float x) => MathTools.Gamma(x + 1.0f);
public static ulong Factorial(uint x)
{
if (x > 20)
throw new ArgumentOutOfRangeException(nameof(x), $"Cannot compute {x}!, since ulong.max is 18_446_744_073_709_551_615.");
ulong accumulator = 1;
for (uint factor = 1; factor <= x; factor++)
accumulator *= factor;
return accumulator;
}
public static ulong Factorial(int x)
{
if(x < 0)
throw new ArgumentOutOfRangeException(nameof(x), "Given value must be greater as zero.");
return MathTools.Factorial((uint) x);
}
}
}

218
FastRng/Float/MultiThreadedRng.cs Normal file
View File

@ -0,0 +1,218 @@
using System;
using System.Diagnostics.CodeAnalysis;
using System.Threading;
using System.Threading.Channels;
using System.Threading.Tasks;
using FastRng.Float.Distributions;
namespace FastRng.Float
{
/// <summary>
/// This class uses the George Marsaglia's MWC algorithm. The algorithm's implementation based loosely on John D.
/// Cook's (johndcook.com) implementation (https://www.codeproject.com/Articles/25172/Simple-Random-Number-Generation).
/// Thanks John for your work.
/// </summary>
public sealed class MultiThreadedRng : IRandom
{
#if DEBUG
private const int CAPACITY_RANDOM_NUMBERS_4_SOURCE = 10_000;
#else
private const int CAPACITY_RANDOM_NUMBERS_4_SOURCE = 16_000_000;
#endif
private readonly CancellationTokenSource producerTokenSource = new CancellationTokenSource();
private readonly object syncUintGenerators = new object();
private readonly object syncUniformDistributedFloatGenerators = new object();
private readonly Thread[] producerRandomUint = new Thread[2];
private readonly Thread[] producerRandomUniformDistributedFloat = new Thread[2];
private uint mW;
private uint mZ;
private readonly Channel<uint> channelRandomUint = Channel.CreateBounded<uint>(new BoundedChannelOptions(CAPACITY_RANDOM_NUMBERS_4_SOURCE)
{
FullMode = BoundedChannelFullMode.Wait,
SingleReader = false,
SingleWriter = false,
});
private readonly Channel<float> channelRandomUniformDistributedFloat = Channel.CreateBounded<float>(new BoundedChannelOptions(CAPACITY_RANDOM_NUMBERS_4_SOURCE)
{
FullMode = BoundedChannelFullMode.Wait,
SingleReader = false,
SingleWriter = false,
});
#region Constructors
public MultiThreadedRng()
{
//
// Initialize the mW and mZ by using
// the system's time.
//
var now = DateTime.Now;
var ticks = now.Ticks;
this.mW = (uint) (ticks >> 16);
this.mZ = (uint) (ticks % 4294967296);
this.StartProducerThreads(deterministic: false);
}
public MultiThreadedRng(uint seedU)
{
this.mW = seedU;
this.mZ = 362436069;
this.StartProducerThreads(deterministic: true);
}
public MultiThreadedRng(uint seedU, uint seedV)
{
this.mW = seedU;
this.mZ = seedV;
this.StartProducerThreads(deterministic: true);
}
private void StartProducerThreads(bool deterministic = false)
{
this.producerRandomUint[0] = new Thread(() => this.RandomProducerUint(this.channelRandomUint.Writer, this.producerTokenSource.Token)) {IsBackground = true};
this.producerRandomUint[1] = new Thread(() => this.RandomProducerUint(this.channelRandomUint.Writer, this.producerTokenSource.Token)) {IsBackground = true};
this.producerRandomUint[0].Start();
if(!deterministic)
this.producerRandomUint[1].Start();
this.producerRandomUniformDistributedFloat[0] = new Thread(() => this.RandomProducerUniformDistributedFloat(this.channelRandomUint.Reader, channelRandomUniformDistributedFloat.Writer, this.producerTokenSource.Token)) {IsBackground = true};
this.producerRandomUniformDistributedFloat[1] = new Thread(() => this.RandomProducerUniformDistributedFloat(this.channelRandomUint.Reader, channelRandomUniformDistributedFloat.Writer, this.producerTokenSource.Token)) {IsBackground = true};
this.producerRandomUniformDistributedFloat[0].Start();
if(!deterministic)
this.producerRandomUniformDistributedFloat[1].Start();
}
#endregion
#region Producers
[ExcludeFromCodeCoverage]
private async void RandomProducerUint(ChannelWriter<uint> channelWriter, CancellationToken cancellationToken)
{
try
{
var buffer = new uint[CAPACITY_RANDOM_NUMBERS_4_SOURCE];
while (!cancellationToken.IsCancellationRequested)
{
lock (syncUintGenerators)
{
for (var n = 0; n < buffer.Length && !cancellationToken.IsCancellationRequested; n++)
{
this.mZ = 36_969 * (this.mZ & 65_535) + (this.mZ >> 16);
this.mW = 18_000 * (this.mW & 65_535) + (this.mW >> 16);
buffer[n] = (this.mZ << 16) + this.mW;
}
}
for (var n = 0; n < buffer.Length && !cancellationToken.IsCancellationRequested; n++)
await channelWriter.WriteAsync(buffer[n], cancellationToken);
}
}
catch (OperationCanceledException)
{
}
}
[ExcludeFromCodeCoverage]
private async void RandomProducerUniformDistributedFloat(ChannelReader<uint> channelReaderUint, ChannelWriter<float> channelWriter, CancellationToken cancellationToken)
{
try
{
var buffer = new float[CAPACITY_RANDOM_NUMBERS_4_SOURCE];
var randomUint = new uint[CAPACITY_RANDOM_NUMBERS_4_SOURCE];
while (!cancellationToken.IsCancellationRequested)
{
for (var n = 0; n < randomUint.Length; n++)
randomUint[n] = await channelReaderUint.ReadAsync(cancellationToken);
lock (syncUniformDistributedFloatGenerators)
for (var n = 0; n < buffer.Length && !cancellationToken.IsCancellationRequested; n++)
buffer[n] = (randomUint[n] + 1.0f) * 2.328306435454494e-10f; // 2.328 => 1/(2^32 + 2)
for (var n = 0; n < buffer.Length && !cancellationToken.IsCancellationRequested; n++)
await channelWriter.WriteAsync(buffer[n], cancellationToken);
}
}
catch (OperationCanceledException)
{
}
}
#endregion
#region Implementing interface
public async ValueTask<float> GetUniform(CancellationToken cancel = default)
{
try
{
return await this.channelRandomUniformDistributedFloat.Reader.ReadAsync(cancel);
}
catch (OperationCanceledException)
{
return float.NaN;
}
}
public async ValueTask<uint> NextNumber(uint rangeStart, uint rangeEnd, IDistribution distribution, CancellationToken cancel = default)
{
if (rangeStart > rangeEnd)
{
var tmp = rangeStart;
rangeStart = rangeEnd;
rangeEnd = tmp;
}
var range = rangeEnd - rangeStart;
distribution.Random ??= this;
var distributedValue = await distribution.GetDistributedValue(cancel);
return (uint) ((distributedValue * range) + rangeStart);
}
public async ValueTask<ulong> NextNumber(ulong rangeStart, ulong rangeEnd, IDistribution distribution, CancellationToken cancel = default(CancellationToken))
{
if (rangeStart > rangeEnd)
{
var tmp = rangeStart;
rangeStart = rangeEnd;
rangeEnd = tmp;
}
var range = rangeEnd - rangeStart;
distribution.Random = this;
var distributedValue = await distribution.GetDistributedValue(cancel);
return (ulong) ((distributedValue * range) + rangeStart);
}
public async ValueTask<float> NextNumber(float rangeStart, float rangeEnd, IDistribution distribution, CancellationToken cancel = default(CancellationToken))
{
if (rangeStart > rangeEnd)
{
var tmp = rangeStart;
rangeStart = rangeEnd;
rangeEnd = tmp;
}
var range = rangeEnd - rangeStart;
distribution.Random = this;
var distributedValue = await distribution.GetDistributedValue(cancel);
return (distributedValue * range) + rangeStart;
}
public async ValueTask<float> NextNumber(IDistribution distribution, CancellationToken cancel = default) => await this.NextNumber(0.0f, 1.0f, distribution, cancel);
public void StopProducer() => this.producerTokenSource.Cancel();
#endregion
}
}

63
FastRng/Float/ShapeFitter.cs Normal file
View File

@ -0,0 +1,63 @@
using System;
using System.Threading;
using System.Threading.Tasks;
using FastRng.Float.Distributions;
namespace FastRng.Float
{
/// <summary>
/// ShapeFitter is a rejection sampler, cf. https://en.wikipedia.org/wiki/Rejection_sampling
/// </summary>
public sealed class ShapeFitter
{
private readonly float[] probabilities;
private readonly IRandom rng;
private readonly float max;
private readonly float sampleSize;
private readonly IDistribution uniform = new Uniform();
public ShapeFitter(Func<float, float> shapeFunction, IRandom rng, ushort sampleSize = 50)
{
this.rng = rng;
this.sampleSize = sampleSize;
this.probabilities = new float[sampleSize];
var sampleStepSize = 1.0f / sampleSize;
var nextStep = 0.0f + sampleStepSize;
var maxValue = 0.0f;
for (var n = 0; n < sampleSize; n++)
{
this.probabilities[n] = shapeFunction(nextStep);
if (this.probabilities[n] > maxValue)
maxValue = this.probabilities[n];
nextStep += sampleStepSize;
}
this.max = maxValue;
}
public async ValueTask<float> NextNumber(CancellationToken token = default)
{
while (!token.IsCancellationRequested)
{
var x = await this.rng.GetUniform(token);
if (float.IsNaN(x))
return x;
var nextBucket = (int)MathF.Floor(x * this.sampleSize);
var threshold = this.probabilities[nextBucket];
var y = await this.rng.NextNumber(0.0f, this.max, this.uniform, token);
if (float.IsNaN(y))
return y;
if(y > threshold)
continue;
return x;
}
return float.NaN;
}
}
}