8
\$\begingroup\$

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);
    }
}
\$\endgroup\$
2
\$\begingroup\$

If I'm not wrong m_stream may be readonly. There are few chances that compiler will optimize code because of this but it's - as usual - more about communication than performance.

You can add one empty line after each constructors and also move { } to their own line(s), it increase readability and hard-disk storage is cheap:

[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()) {
}

I see you're using CLSCompliant attribute. If you really care about other .NET languages without unsigned integers support then you may want to add a specific constructor just for them:

public PcgRandom(int state, int stream) {
    checked {
        m_state = (ulong)state;
        m_stream = stream ¦ 1;
    }
}

Performance shouldn't be an issue while creating the object then I'm explicitly asking for a checked conversion.

In GetUInt32(int, int) you need to cast range to uint just because you're using + 1UL, just change + 1UL to + 1.) Also else after a return is usually redundant. For unsigned numbers overflow should do the job. Let's also see the similar GetInt32():

public int GetInt32(int x, int y) {
    var min = Math.Min(x, y);
    var max = Math.Max(x, y);
    var range = max + 1 - min;

    if (uint.MaxValue != range) 
        return (int)(GetUInt32((uint)range) + min);

    return GetInt32();
}

If you're really looking for performance then you should avoid those calls to Math.Min() and Math.Max(). You may be surprised about branchless code JIT compiler may generate (when CMOV & friends are supported):

int min = x, max = y;
if (min > max) {
    min = y;
    max = x;
}

Or maybe even:

int min = x < y ? x : y;
int max = x < y ? y : x;

Unfortunately benchmarking this can be tricky because generated code might be different for different CPUs. Just try all of them and pick the easier one to read...What about Math.Abs()? Given that:

$$ |a-b| = \max(a,b) - \min(a,b) $$

And that compiler is pretty good to avoid branches we may also want to try performance of this (pseudo-code without casts and only for full range):

public int GetInt32(int x, int y) {
    return GetUInt32(Math.Abs(x - y + 1)) + Math.Min(x, y);
}

You may also try to replace Math.Min() with the usual (x < y ? x : y).

\$\endgroup\$
  • \$\begingroup\$ (m_stream | 1) isn't a micro-optimization, it's there to guarantee that the value is odd; a requirement of the PCG construction. Though, this probably means that I should make that particular constructor private since a user might be surprised that all even streams give the same result as the next odd stream. The only reason I didn't mark m_stream as readonly is because I eventually want to expose the stream functionality to end users but haven't figured out how to do so in a way that hides the implementation details. \$\endgroup\$ – Kittoes0124 Apr 25 '17 at 16:31
  • \$\begingroup\$ Correct me if I'm wrong but the GetUInt32(uint x, uint y)requires the cast to ulong in order to prevent integer wraparound in the case where uint.MinValue and uint.MaxValue are passed as arguments. Though, I suppose there's probably a way to allow the wraparound and still get the same results if I change the if... \$\endgroup\$ – Kittoes0124 Apr 25 '17 at 16:37
  • \$\begingroup\$ True about odd number, I didn't do my homework! Casting should not matter for an unchecked operation (it will just overflow) but I already turned off my brain then... \$\endgroup\$ – Adriano Repetti Apr 25 '17 at 16:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.