Generating uniform random integers with a certain maximum

Question

I want to generate uniform integers that satisfy 0 <= result <= maxValue.

I already have a generator that returns uniform values in the full range of the built in unsigned integer types. Let's call the methods for this byte Byte(), ushort UInt16(), uint UInt32() and ulong UInt64(). Assume that the result of these methods is perfectly uniform.

The signature of the methods I want are uint UniformUInt(uint maxValue) and ulong UniformUInt(ulong maxValue).

What I'm looking for:

1. Correctness
   I'd prefer the return values to be distributed in the given interval.
   But a very small bias is acceptable if it increases performance significantly. By that I mean a bias of an order that allows distinguisher with probability 2/3 given 2^64 values.
   It must work correctly for any maxValue.
2. Performance
   The method should be fast.
3. Efficiency
   The method does consume little raw randomness, since depending on the underlying generator, generating the raw bytes might be costly. Wasting a few bits is fine, but consuming say 128 bits to generate a single number is probably excessive.

It's also possible to cache some left over randomness from the previous call in some member variables.

Be careful with int overflows, and wrapping behavior.

I already have a solution(I'll post it as an answer), but it's a bit ugly for my tastes. So I'd like to get ideas for better solutions.

Suggestions on how to unit test with large maxValues would be nice too, since I can't generate a histogram with 2^64 buckets and 2^74 random values. Another complication is that with certain bugs, only some maxValue distributions are biased a lot, and others only very slightly.

Roscoeroscommon

What is the rationale for not using floor(maxValue/maxfordatatype*existingrandom()) ? – Ergocalciferol

@EugenRieck It violates the correctness property, unless maxfordatatype is significantly larger than maxValue, in which case it'd have efficiency issues. – Roscoeroscommon

I understood from your OQ, that a small bias is acceptable, so I should understand, that an 1/maxValue bias is NOT acceptably small ? – Ergocalciferol

@EugenRieck since each return value only has a probability of 1/maxValue, a bias of the same magnitude means that some fields might not be reached at all, and others twice as often as they should. System.Random uses a similar algorithm in Next(maxValue), and I can distinguish it from true random numbers using only a few samples. – Roscoeroscommon

Note that I quantified small bias, so that even a perfect distinguisher would need 2^64 samples answer correctly in 2/3 of the cases when asked: Is this sample set truly random, or was it generated by the generator with that know bias. So my requirement is very strict, and I prefer not having a bias at all. – Roscoeroscommon

One more: Can 64bit operations be assumed to be efficient? – Ergocalciferol

Since this is intended to be a general purpose library, I like being fast on at least x86 and x64. – Roscoeroscommon

Answer 1

How about something like this as a general-purpose solution? The algorithm is based on that used by Java's nextInt method, rejecting any values that would cause a non-uniform distribution. So long as the output of your UInt32 method is perfectly uniform then this should be too.

uint UniformUInt(uint inclusiveMaxValue)
{
    unchecked
    {
        uint exclusiveMaxValue = inclusiveMaxValue + 1;

        // if exclusiveMaxValue is a power of two then we can just use a mask
        // also handles the edge case where inclusiveMaxValue is uint.MaxValue
        if ((exclusiveMaxValue & (~exclusiveMaxValue + 1)) == exclusiveMaxValue)
            return UInt32() & inclusiveMaxValue;

        uint bits, val;
        do
        {
            bits = UInt32();
            val = bits % exclusiveMaxValue;

            // if (bits - val + inclusiveMaxValue) overflows then val has been
            // taken from an incomplete chunk at the end of the range of bits
            // in that case we reject it and loop again
        } while (bits - val + inclusiveMaxValue < inclusiveMaxValue);

        return val;
    }
}

The rejection process could, theoretically, keep looping forever; in practice the performance should be pretty good. It's difficult to suggest any generally applicable optimisations without knowing (a) the expected usage patterns, and (b) the performance characteristics of your underlying RNG.

For example, if most callers will be specifying a max value <= 255 then it might not make sense to ask for four bytes of randomness every time. On the other hand, the performance benefit of requesting fewer bytes might be outweighed by the additional cost of always checking how many you actually need. (And, of course, once you do have specific information then you can keep optimising and testing until your results are good enough.)

Contumelious

(I should also point out that I haven't tested this particular piece of code for correctness or bias. I just adapted it from some older code that I had knocking around. The other code had been tested once-upon-a-time, so this should hopefully be ok too.) – Contumelious

My solution uses a similar rejection based algorithm, but I calculate the rejection bound explicitly. So your rejection condition is probably slightly faster. – Roscoeroscommon

I'll probably add the same cases like in my own code, do reduce entropy consumption. My profiling indicates that the special casing is worthwhile, especially when using slower randomness providers. – Roscoeroscommon

Answer 2

I am not sure, that his is an answer. It definitly needs more space than a comment, so I have to write it here, but I am willing to delete if others think this is stupid.

From the OQ I get, that

1. Entropy bits are very expensive
2. Everything else should be considered expensive, but less so than entropy.

My idea is to use binary digits to half, quater ... the maxValue space, until it is reduced to a number. Somthing like

I'l use maxValue=333 (decimal) as an example and assume a function getBit(), that randomly returns 0 or 1

offset:=0
space:=maxValue

while (space>0)

  //Right-shift the value, keeping the rightmost bit this should be 
  //efficient on x86 and x64, if coded in real code, not pseudocode
  remains:=space & 1
  part:=floor(space/2)
  space:=part

  //In the 333 example, part is now 166, but 2*166=332 If we were to simply chose one
  //half of the space, we would be heavily biased towards the upper half, so in case
  //we have a remains, we consume a bit of entropy to decide which half is bigger

  if (remains)
    if(getBit())
      part++;

  //Now we decide which half to chose, consuming a bit of entropy
  if (getBit())
    offset+=part;

  //Exit condition: The remeinind number space=0 is guaranteed to be met
  //In the 333 example, offset will be 0, 166 or 167, remaining space will be 166
}

randomResult:=offset

getBit() can either come from your entropy source, if it is bit-based, or by consuming n bits of entropy at once on first call (obviously with n being the optimum for your entropy source), and shifting this until empty.

Ergocalciferol

very expensive is a bit of an overstatement. It's up to 40 clocks per byte depending on the chosen provider. But with other providers it's much faster. – Roscoeroscommon

If it is much faster, you can use the simple division method for 8 and 16 bit starting from a 32 resp. 64 bit source, but for examples like 333333 I stand to my concept. – Ergocalciferol

Certainly an interesting implementation. Not sure how fast it is. I'll need to benchmark. – Roscoeroscommon

As I said in the preface, I thought it more important to avoid using excess entropy, but I suspect a shift with carry, 2-3 if, 0-2 additions in place per ceil(ld(maxValue)) to be reasonably fast. Much will depend on the implementation of getBit() – Ergocalciferol

Answer 3

My current solution. A bit ugly for my tastes. It also has two divisions per generated number, which might negatively impact performance (I haven't profiled this part yet).

uint UniformUInt(uint maxResult)
{
    uint rand;
    uint count = maxResult + 1;

    if (maxResult < 0x100)
    {
        uint usefulCount = (0x100 / count) * count;
        do
        {
            rand = Byte();
        } while (rand >= usefulCount);
        return rand % count;
    }
    else if (maxResult < 0x10000)
    {
        uint usefulCount = (0x10000 / count) * count;
        do
        {
            rand = UInt16();
        } while (rand >= usefulCount);
        return rand % count;
    }
    else if (maxResult != uint.MaxValue)
    {
        uint usefulCount = (uint.MaxValue / count) * count;//reduces upper bound by 1, to avoid long division
        do
        {
            rand = UInt32();
        } while (rand >= usefulCount);
        return rand % count;
    }
    else
    {
        return UInt32();
    }
}

ulong UniformUInt(ulong maxResult)
{
    if (maxResult < 0x100000000)
        return InternalUniformUInt((uint)maxResult);
    else if (maxResult < ulong.MaxValue)
    {
        ulong rand;
        ulong count = maxResult + 1;
        ulong usefulCount = (ulong.MaxValue / count) * count;//reduces upper bound by 1, since ulong can't represent any more
        do
        {
            rand = UInt64();
        } while (rand >= usefulCount);
        return rand % count;
    }
    else
        return UInt64();
}

Roscoeroscommon