Computer Science Stack Exchange is a question and answer site for students, researchers and practitioners of computer science. It only takes a minute to sign up.
Sign up to join this community
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.
(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)}$.
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;
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 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;
}
}
}
Generate a binary decimal. Instead of storing it explicitly, just keep track of the min and max possible values. Once those values lie within the same integer, return that integer. Sketch of the code is below.
(Edit) Fuller explanation: Say you want to generate a random integer from 1 to 3 inclusive with 1/3 probability each. We do this by generating a random binary decimal real, x in the range (0, 1). If x < 1/3, return 1, else if x < 2/3 return 2, else return 3. Instead of generating the digits for x explicitly, we just keep track of the minimum and maximum possible values of x. Initially, the min value of x is 0 and the max is 1. If you flip heads first then the first digit of x behind the decimal point (in binary) is 1. The minimum possible value of x (in binary) then becomes 0.100000 = 1/2 and the max is 0.111111111 = 1. Now if your next flip is tails x starts with 0.10. The min possible value is 0.1000000 = 1/2 and the max is 0.1011111 = 3/4. The min possible value of x is 1/2 so you know there's no chance of returning 1 since that requires x < 1/3. You can still return 2 if x ends up being 1/2 < x < 2/3 or 3 if 2/3 < x < 3/4. Now suppose the third flip is tails. Then x must start with 0.100. Min = 0.10000000 = 1/2 and max = 0.100111111 = 5/8. Now since 1/3 < 1/2 < 5/8 < 2/3 we know that x must fall in the interval (1/3, 2/3), so we can stop generating digits of x and just return 2.
The code does essentially this except instead of generating x between 0 and 1 it generates x between a and b, but the principle is the same.
def gen(a, b):
min_possible = a
max_possible = b
while True:
floor_min_possible = floor(min_possible)
floor_max_possible = floor(max_possible)
if max_possible.is_integer():
floor_max_possible -= 1
if floor_max_possible == floor_min_possible:
return floor_max_possible
mid = (min_possible + max_possible)/2
if coin_flip():
min_possible = mid
else:
max_possible = mid
Remark: I tested this code against the accept/reject method and both yield uniform distributions. This code requires less coin flips than accept reject except when b - a is close to the next power of 2. E.x. if you want to generate a = 0, b = 62 then accept/reject does better. I was able to prove that this code can use on average no more than 2 more coin flips than accept/reject. From my reading, it looks like Knuth and Yao (1976) gave a method to solve this problem and proved that their method is optimal in the expected number of coin flips. They further proved that expected number of flips must be greater than the Shannon entropy of the distribution. I could not find a copy of the text of the paper however and would be curious to see what their method is. (Update: just found an exposition of Knuth Yao 1976 here: http://www.nrbook.com/devroye/Devroye_files/chapter_fifteen_1.pdf but I have not yet read it). Someone also mentioned Han Hoshi in this thread which seems to be more general and solves it using a biased coin. See also http://paper.ijcsns.org/07_book/200909/20090930.pdf by Pae (2009) for a good discussion of the literature.
Theoretically optimal algorithm
Here is an improvement of the other answer I posted. The other answer has the advantage that it is easier to extend to the more general case of generating one discrete distribution from another. In fact, the other answer is a special case of the algorithm due to Han and Hoshi.
The algorithm I will describe here is based on Knuth and Yao (1976). In their paper, they also proved that this algorithm achieves the minimum possible expected number of coin flips.
To illustrate it, consider the Rejection sampling method described by other answers. As an example, suppose you want to generate one of 5 numbers uniformly [0, 4]. The next power of 2 is 8 so you flip the coin 3 times and generate a random number up to 8. If the number is 0 to 4 then you return it. Otherwise, you throw it out and generate another number up to 8 and try again until you succeed. But when you throw out the number, you just wasted some entropy. You can instead condition on the number you threw out to reduce the number of future coin flips you'll need in expectation. Concretely, once you generate the number [0, 7], if it is [0, 4], return. Otherwise, it's 5, 6, or 7, and you do something different in each case. If it's 5, flip the coin again and return either 0 or 1 based on the flip. If it's 6, flip the coin and return either 2 or 3. If it's 7, flip the coin; if it's heads, return 4, if it's tails start over.
The leftover entropy from our initial failed attempt gave us 3 cases (5, 6, or 7). If we just throw this out, we throw away log2(3) coin flips. We instead keep it, and combine it with the outcome of another flip to generate 6 possible cases (5H, 5T, 6H, 6T, 7H, 7T) which let's us immediately try again to generate a final answer with probability of success 5/6.
Here is the code:
# returns an int from [0, b)
def __gen(b):
rand_num = 0
num_choices = 1
while True:
num_choices *= 2
rand_num *= 2
if coin.flip():
rand_num += 1
if num_choices >= b:
if rand_num < b:
return rand_num
num_choices -= b
rand_num -= b
# returns an int from [a, b)
def gen(a, b):
return a + __gen(b - a)
The simple answer?
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 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
If $2^{n-1} < b - a ≤ 2^n$ then pick a number x, $0 ≤ x < 2^n$, using n coin throws to determine the highest, second highest etc. bit of x. If x < b-a then return a + x otherwise try again. As an optimisation, try again as soon as the throws so far tell you that $x < b - a$ is not possible. For example if b - a = 33, and the first coin is 1 = 32, and the second coin is 1 = 16, then x ≥ 48 so restart already.
