From c340635c4a2d036d5dd91b1ad0c708fa07df1aa9 Mon Sep 17 00:00:00 2001 From: Thorsten Sommer Date: Sat, 31 Oct 2020 23:27:18 +0100 Subject: [PATCH] Added float implementations --- FastRng/Float/Distributions/IDistribution.cs | 12 + FastRng/Float/Distributions/Uniform.cs | 12 + FastRng/Float/IRandom.cs | 21 ++ FastRng/Float/MathTools.cs | 61 ++++++ FastRng/Float/MultiThreadedRng.cs | 218 +++++++++++++++++++ FastRng/Float/ShapeFitter.cs | 63 ++++++ 6 files changed, 387 insertions(+) create mode 100644 FastRng/Float/Distributions/IDistribution.cs create mode 100644 FastRng/Float/Distributions/Uniform.cs create mode 100644 FastRng/Float/IRandom.cs create mode 100644 FastRng/Float/MathTools.cs create mode 100644 FastRng/Float/MultiThreadedRng.cs create mode 100644 FastRng/Float/ShapeFitter.cs diff --git a/FastRng/Float/Distributions/IDistribution.cs b/FastRng/Float/Distributions/IDistribution.cs new file mode 100644 index 0000000..9b80c08 --- /dev/null +++ b/FastRng/Float/Distributions/IDistribution.cs @@ -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 GetDistributedValue(CancellationToken token); + } +} \ No newline at end of file diff --git a/FastRng/Float/Distributions/Uniform.cs b/FastRng/Float/Distributions/Uniform.cs new file mode 100644 index 0000000..fc593a2 --- /dev/null +++ b/FastRng/Float/Distributions/Uniform.cs @@ -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 GetDistributedValue(CancellationToken token = default) => this.Random == null ? float.NaN : await this.Random.GetUniform(token); + } +} \ No newline at end of file diff --git a/FastRng/Float/IRandom.cs b/FastRng/Float/IRandom.cs new file mode 100644 index 0000000..af8d319 --- /dev/null +++ b/FastRng/Float/IRandom.cs @@ -0,0 +1,21 @@ +using System.Threading; +using System.Threading.Tasks; +using FastRng.Float.Distributions; + +namespace FastRng.Float +{ + public interface IRandom + { + public ValueTask GetUniform(CancellationToken cancel = default); + + public ValueTask NextNumber(uint rangeStart, uint rangeEnd, IDistribution distribution, CancellationToken cancel = default); + + public ValueTask NextNumber(ulong rangeStart, ulong rangeEnd, IDistribution distribution, CancellationToken cancel = default); + + public ValueTask NextNumber(float rangeStart, float rangeEnd, IDistribution distribution, CancellationToken cancel = default); + + public ValueTask NextNumber(IDistribution distribution, CancellationToken cancel = default); + + public void StopProducer(); + } +} \ No newline at end of file diff --git a/FastRng/Float/MathTools.cs b/FastRng/Float/MathTools.cs new file mode 100644 index 0000000..63115a4 --- /dev/null +++ b/FastRng/Float/MathTools.cs @@ -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); + } + } +} \ No newline at end of file diff --git a/FastRng/Float/MultiThreadedRng.cs b/FastRng/Float/MultiThreadedRng.cs new file mode 100644 index 0000000..8b1d178 --- /dev/null +++ b/FastRng/Float/MultiThreadedRng.cs @@ -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 +{ + /// + /// 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. + /// + 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 channelRandomUint = Channel.CreateBounded(new BoundedChannelOptions(CAPACITY_RANDOM_NUMBERS_4_SOURCE) + { + FullMode = BoundedChannelFullMode.Wait, + SingleReader = false, + SingleWriter = false, + }); + + private readonly Channel channelRandomUniformDistributedFloat = Channel.CreateBounded(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 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 channelReaderUint, ChannelWriter 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 GetUniform(CancellationToken cancel = default) + { + try + { + return await this.channelRandomUniformDistributedFloat.Reader.ReadAsync(cancel); + } + catch (OperationCanceledException) + { + return float.NaN; + } + } + + public async ValueTask 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 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 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 NextNumber(IDistribution distribution, CancellationToken cancel = default) => await this.NextNumber(0.0f, 1.0f, distribution, cancel); + + public void StopProducer() => this.producerTokenSource.Cancel(); + + #endregion + } +} \ No newline at end of file diff --git a/FastRng/Float/ShapeFitter.cs b/FastRng/Float/ShapeFitter.cs new file mode 100644 index 0000000..025f583 --- /dev/null +++ b/FastRng/Float/ShapeFitter.cs @@ -0,0 +1,63 @@ +using System; +using System.Threading; +using System.Threading.Tasks; +using FastRng.Float.Distributions; + +namespace FastRng.Float +{ + /// + /// ShapeFitter is a rejection sampler, cf. https://en.wikipedia.org/wiki/Rejection_sampling + /// + 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 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 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; + } + } +} \ No newline at end of file