Edit: Dunno what PCG is and don't want to read the paper? Maybe this video of Melissa O'Neill (the author) explaining things will be palatable instead.
Original: I attempted to answer this question last week and was thrown off when someone mentioned in a comment that my method would show bias. Eager to prove the naysayer wrong I discovered that things were even worse than he suggested as certain ranges will cause things to break down entirely. So I decided to start from scratch, do some more research, and try again.
Performance was also much more important to me this time around so I decided to try and find a way to avoid calling RandomNumberGenerator.GetBytes() to generate every single value. Eventually I settled on PCG as an appropriate modern algorithm. Floating point math has also been entirely avoided while working with integers in order to avoid the sorts of issues I encountered in my previous attempts.
Looking for potential ways to speed things up or simplify the code; the branch statements bother me quite a bit but can't envision a better way to include the full range of values without them (the upper bound cannot be generated by theGetUInt32(uint ExclusiveHigh) method). Haven't been able to find any bugs so far but that doesn't mean they're not lurking either...
public class PcgRandom
{
private const double INVERSE_32_BIT = 2.32830643653869628906e-010d;
private const double INVERSE_52_BIT = 2.22044604925031308085e-016d;
private static RandomNumberGenerator Rng = new RNGCryptoServiceProvider();
private ulong m_state;
private ulong m_stream;
[CLSCompliant(false)]
public PcgRandom(ulong state, ulong stream) {
m_state = state;
m_stream = (stream | 1UL);
}
[CLSCompliant(false)]
public PcgRandom(ulong state) : this(state, GetSeed()) { }
public PcgRandom() : this(GetSeed(), GetSeed()) { }
/// <summary>
/// Generates a uniformly distributed double between the range (0, 1).
/// </summary>
public double GetDouble() {
return CreateDouble(GetInt32(), GetInt32());
}
/// <summary>
/// Generates a uniformly distributed 32-bit signed integer between the range of int.MaxValue and int.MinValue.
/// </summary>
public int GetInt32() {
return ((int)GetUInt32());
}
/// <summary>
/// Generates a uniformly distributed 32-bit signed integer between the range [min, max].
/// </summary>
public int GetInt32(int x, int y) {
var min = Math.Min(x, y);
var max = Math.Max(x, y);
var range = (max + 1L - min);
if (uint.MaxValue > range) {
return ((int)(GetUInt32((uint)range) + min));
}
else {
return GetInt32();
}
}
/// <summary>
/// Generates a uniformly distributed 32-bit unsigned integer between the range of uint.MaxValue and uint.MinValue.
/// </summary>
[CLSCompliant(false)]
public uint GetUInt32() {
return Pcg32(ref m_state, m_stream);
}
/// <summary>
/// Generates a uniformly distributed 32-bit unsigned integer between the range [min, max].
/// </summary>
[CLSCompliant(false)]
public uint GetUInt32(uint x, uint y) {
var min = Math.Min(x, y);
var max = Math.Max(x, y);
var range = (max + 1UL - min);
if (uint.MaxValue > range) {
return (GetUInt32((uint)range) + min);
}
else {
return GetUInt32();
}
}
private uint GetUInt32(uint exclusiveHigh) {
var threshold = ((uint)((0x100000000UL - exclusiveHigh) % exclusiveHigh));
var sample = GetUInt32();
while (sample < threshold) {
sample = GetUInt32();
}
return (sample % exclusiveHigh);
}
private static double CreateDouble(int x, int y) {
// reference: https://www.doornik.com/research/randomdouble.pdf
return (0.5d + (INVERSE_52_BIT / 2) + (x * INVERSE_32_BIT) + ((y & 0x000FFFFF) * INVERSE_52_BIT));
}
private static ulong GetSeed() {
var buffer = new byte[sizeof(ulong)];
Rng.GetBytes(buffer);
return BitConverter.ToUInt64(buffer, 0);
}
private static uint Pcg32(ref ulong state, ulong stream) {
// reference: http://www.pcg-random.org/paper.html
state = unchecked(state * 6364136223846793005UL + stream);
return RotateRight((uint)(((state >> 18) ^ state) >> 27), (int)(state >> 59));
}
private static uint RotateRight(uint value, int count) {
return ((value >> count) | (value << ((-count) & 31)));
}
}
Distribution Results:
Distribution Code:
class Program
{
static void Main(string[] args) {
var results = Sample(2000000000);
for (var i = 0; i < results.Length; i++) {
Console.WriteLine($" N: {i} | Count: {results[i]}");
}
Console.ReadKey();
}
static int[] Sample(int count) {
var pcg = new PcgRandom();
var results = new int[10];
Parallel.For(0, count, (index) => {
Interlocked.Increment(ref results[pcg.GetInt32(0, 9)]);
});
return results;
}
}
Benchmark.NET Results:
Benchmark.NET Code:
class Program
{
static void Main(string[] args) {
var summary = BenchmarkRunner.Run<RandomComparison>();
Console.ReadKey();
}
}
public class RandomComparison
{
private static Random m_random = new Random();
private static PcgRandom m_pcg = new PcgRandom();
[Benchmark()]
public double DotNetRandom_Double() {
return m_random.NextDouble();
}
[Benchmark()]
public int DotNetRandom_Int32_RangeMax() {
return m_random.Next(int.MinValue, int.MaxValue);
}
[Benchmark()]
public int DotNetRandom_Int32_RangeHalf() {
return m_random.Next(0, int.MaxValue);
}
[Benchmark()]
public double PcgRandom_Double() {
return m_pcg.GetDouble();
}
[Benchmark()]
public int PcgRandom_Int32_RangeMax() {
return m_pcg.GetInt32();
}
[Benchmark()]
public int PcgRandom_Int32_RangeHalf() {
return m_pcg.GetInt32(0, int.MaxValue);
}
}
