# Generating uniformly distributed random numbers using a coin

You have one coin. You may flip it as many times as you want.

You want to generate a random number $r$ such that $a \leq r < b$ where $r,a,b\in \mathbb{Z}^+$.

Distribution of the numbers should be uniform.

It is easy if $b -a = 2^n$:

r = a + binary2dec(flip n times write 0 for heads and 1 for tails)


What if $b-a \neq 2^n$?

-

What you're looking for is based on Rejection sampling or the accept-reject method (note that the Wiki page is a bit technical).

This method is useful in these kinds of situations: you want to pick some random object from a set (a random integer in the set $[a,b[$ in your case), but you don't know how to do that, but you can pick some random object from a larger set containing the first set (in your case, $[a, 2^k + a[$ for some $k$ such that $2^k + a \ge b$; this corresponds to $k$ coin flips).

In such a scenario, you just keep picking elements from the bigger set until you randomly pick an element in the smaller set. If your smaller set is big enough compared to your larger set (in your case, $[a, 2^k + a[$ contains at most twice as many integers as $[a,b[$, which is good enough), this is efficient.

An alternative example: suppose you want to pick a random point inside a circle with radius 1. Now, it isn't really easy to come up with a direct method for this. We turn to the accept-reject method: we sample points in a 1x1 square encompassing the circle, and test if the number we draw lies inside the circle.

-
Note that if we reject samples from $A$ in order to get a distribution on $B$, the expected number of iterations is $\frac{|A|}{|B|}$ (as we perform an experiment with geometric distribution). – Raphael Mar 21 '12 at 15:23
I recall seeing somewhere that this can't be done exactly unless the range is a power of 2 (as stands to reason, e.g. 1 / 3 has no terminating binary expansion). – vonbrand Jan 25 at 14:22

(technicalities: the answer fits selection of number $a \le x < b$)

Since you are allowed to flip your coin as many times as you wish, you can get your probability as-close-as-you-wish to uniform by picking a fraction $r\in [0,1]$ (using binary radix: you flip the coin for each digit after the point) and multiply $r$ by $b-a$ to get a number between 0 and [b-a-1] (rounding down to an integer). Add this number to $a$ and you're done.

Example: say $b-a=3$. 1/3 in binary is 0.0101010101.... Then, if your flip is between 0 and 0.010101... your pick is $b$. if it is beween 0.010101.. and 0.10101010... your pick will be $a+1$, and otherwise it will be $a+2$.

If you flip your coin $t$ times then each number between $a$ and $b$ will be chosen with probability $\frac{1}{b-a}\pm 2^{-(t+1)}$.

-
This doesn't give a uniform distribution. For some applications (e.g. crypto, sometimes), this can be very bad. – Gilles Mar 21 '12 at 9:58
@Gilles: It can be fixed to give a perfectly uniform distribution by flipping until it's no longer possible for the result to change. This is the most efficient answer. – Neil G Mar 21 '12 at 11:24
@NeilG I know it can be fixed, but fixing it would be a crucial part of the answer. – Gilles Mar 21 '12 at 12:33
@Gilles: You're right. He could modify the answer to say that you can produce a perfectly uniform distribution if you flip while $\lfloor (b-a)(f + 2^{-t-1})\rfloor \ne \lfloor (b-a)(f - 2^{-t-1})\rfloor$. +1 from me for the best average case and worst case time. – Neil G Mar 21 '12 at 12:45
@NeilG, it can't be "fixed", as there is quite a sizeable set of integers which don't have a terminating binary fraction. – vonbrand Mar 20 at 20:02
show 1 more comment

Choose a number in the next larger power of 2 range, and discard answers greater than $b-a$.

n = b-a;
N = round_to_next_larger_power_of_2(n)
while (1) {
x = random(0 included to N excluded);
if (x < n) break;
}
r = a + x;

-
And why does this work? – Raphael Mar 21 '12 at 13:20
@Raphael are you skeptical, or do you merely want the poster to explain in more detail ? – Suresh Mar 21 '12 at 15:30
@Suresh: The latter. The pseudo code could be polished a bit, but it implements what other answerers explain. Without justification, this answer is not worth much on its own. – Raphael Mar 21 '12 at 15:39

If b-a is not a power of 2 then you may have to flip many coins in order to get a result. You may even never get a result, but that is unlikely in the extreme.

Methods

The simplest method is to generate a number in [a, a+2^n), where 2^n >= b-a, until one happens to land in [a, b). You throw away a lot of entropy with this method.

A more expensive method allows you to keep all of the entropy, but becomes very expensive computationally as the number of coin flips / dice rolls increases. Intuitively it is like treating the coin flips as the digits of a binary number to the right of the decimal point, converting that number from base 2 to base a-b after, and returning the digits of that number as they become 'stuck'.

Example

The following code converts rolls of a fair n-sided die into rolls of a fair m-sided die (in your case n=2, m=a-b), with increasing marginal cost as the number of rolls increases. Note the need for a Rational number type with arbitrary precision. One nice property is that converting from n-sided to m-sided and back to n-sided will return the original stream, though perhaps delayed by a couple rolls due to digits having to become stuck.

public static IEnumerable<BigInteger> DigitConversion(this IEnumerable<BigInteger> inputStream, BigInteger modIn, BigInteger modOut) {
//note: values are implicitly scaled so the first unfixed digit of the output ranges from 0 to 1
Rational b = 0; //offset of the chosen range
Rational d = 1; //size of the chosen range
foreach (var r in inputStream) {
//narrow the chosen range towards the real value represented by the input
d /= modIn;
b += d * r;
//check for output digits that have become fixed
while (true) {
var i1 = (b * modOut).Floor();
var i2 = ((b + d) * modOut).Floor(); //note: ideally b+d-epsilon, but another iteration makes that correction unnecessary
if (i1 != i2) break; //digit became fixed?
//fix the next output digit (rescale the range to make next digit range from 0 to 1)
d *= modOut;
b *= modOut;
b -= i1;
yield return i1;
}
}
}

-
 "but that is unlikely in the extreme" -- the probability for this event is $0$; we say it "almost surely" does not happen. – Raphael♦ Mar 21 '12 at 15:21

Nobody mentioned this, so let me formally prove that if $D$ is a domain whose size is not a power of two, then finitely many fair coin tosses aren't enough to generate a uniformly random member of $D$. Suppose you use $k$ coins to generate a member of $D$. For each $d \in D$, the probability that your algorithm generated $d$ is of the form $A/2^k$ for some integer $A$. The fundamental theorem of arithmetic shows that $A/2^k \neq 1/{|D|}$.

If you want to generate $n$ independent uniform samples of $D$, then the expected number of coin tosses you need (under the optimal algorithm) is roughly $n \log_2 |D|$. More generally, if you want to generate from a distribution of entropy $H$, you need roughly $nH$ random bits in expectation. An algorithm achieving this is arithmetic decoding, applied to an (ostensibly) infinite sequence of random bits.

-

If $b-a$ is not a power of 2, then after generating $r$, check whether it is in-range, and if not, generate it again.

The most likely time you will have to re-generate $r$ is when $b - a = 2^n + 1$ for some integer $n$, but even then it will have a greater than 50% chance of falling within range on each generation.

-
 This doesn't seem complete. – Dave Clarke Mar 21 '12 at 9:47

This is a proposed solution for the case when b - a does not equal 2^k. It is supposed to work in fixed number of steps (no need to throw away candidates which are outside your expected range).

However, I am not sure this correct. Please critique and help describe the exact non-uniformity in this random number generator (if any), and how to measure/quantify it.

Firstly convert to equivalent problem of generating uniformly distributed random numbers in range [0, z-1] where z = b - a.

Also, let m = 2^k be the smallest power of 2 >= z.

As per the solution above, we already have a uniformly distributed random number generator R(m) in range [0,m-1] (can be done by tossing k coins, one for each bit).

    Keep a random seed s and initialize with s = R(m).

function random [0, z-1] :
x = R(m) + s
while x >= z:
x -= z
s = x
return x


The while loop runs at most 3 times, giving next random number in fixed number of steps (best case = worst case).

See a test program for numbers [0,2] here: http://pastebin.com/zuDD2V6H

-
 This is not uniform. Take $z=3$. You get $m=4$ and the probabilities are $1/2,1/4,1/4$. – Yuval Filmus Jan 24 at 21:49 Please have a look at the pseudo code as well as the linked code more closely. It does emit 0, 1, and 2 almost with equal frequency... – vpathak Jan 25 at 2:10 You're right, when you look at the outputs individually, they are uniform (the stationary probability is uniform). But they're not independent. If you have just output $0$, say, then the next output is zero with probability $1/2$ and one or two with probability $1/4$ each. – Yuval Filmus Jan 25 at 4:08 You can replace the entire function with a single line: return s = (s + R(m)) % z; – Yuval Filmus Jan 25 at 4:10