Sampling a biased coin with a fair coin

August 15, 2015revised September 20, 2026

How would you use a fair coin to choose something with probability $1/3$?

Flip until the first heads. If it arrives on an even-numbered flip, choose the thing. Otherwise, don’t. The successful sequences are TH, TTTH, TTTTTH, and so on, so their probabilities add up to

$$\frac14+\frac1{16}+\frac1{64}+\cdots =\frac{1/4}{1-1/4}=\frac13.$$

There is no rejection and restart. Every terminating run supplies an answer, including a heads on the first flip. The chance of still flipping after $k$ tosses is $2^{-k}$.

# Let the binary digits choose the outcome

Write the desired probability as

$$p=(0.d_1d_2d_3\ldots)_{2} =\sum_{k=1}^{\infty}\frac{d_k}{2^k}, \qquad d_k\in\{0,1\}.$$

Flip a fair coin until the first heads. If that happens on flip $K$, return $d_K$. A returned 1 means success; a returned 0 means failure. This is a Bernoulli trial with parameter $p$.

The first heads occurs at position $k$ with probability $2^{-k}$, so

$$\Pr(\text{return }1) =\sum_{k=1}^{\infty}\Pr(K=k)d_k =\sum_{k=1}^{\infty}2^{-k}d_k=p.$$

The construction assigns the probability mass $1/2$ to the first digit, $1/4$ to the second, $1/8$ to the third, and so on. Those are exactly the place values in the binary expansion. For $1/3=(0.010101\ldots)_{2}$, the even-numbered digits are 1.

Reproduce this figure

Download source bundleextract, then runtar -xf coin-sampling-gen_figures.tar && cd coin-sampling-gen_figures/coin-sampling && uv run --locked gen_figures.py

Requires uv and Python 3.12; the locked packages download on first run.

Individual files:gen_figures.pylockfilefigure_quality.py

Lumbroso gives this geometric-index construction and an exact rational implementation in Appendix B of his paper on random-bit generation. It is an established sampling method, rather than a special property of one third.1

# Two flips on average, but no fixed limit

For any nonnegative integer-valued running time, its expectation is the sum of its tail probabilities. Each run of length $K$ contributes one to the first $K$ terms of that sum. Here this gives

$$\mathbb E[K] =\sum_{k=0}^{\infty}\Pr(K>k) =\sum_{k=0}^{\infty}2^{-k}=2.$$

This counts fair flips. Computing a requested digit of $p$ is separate work. For rational probabilities the digit calculation is a short integer operation; an arbitrary real number is not automatically an executable input.

The algorithm terminates with probability one. That does not give a fixed maximum runtime: after 20 flips the chance of still waiting is $2^{-20}$, about one in a million.

That distinction is unavoidable for $p=1/3$. Any algorithm that always finishes within $m$ fair flips can be described using $2^m$ equally likely bit strings. Pad shorter runs with ignored flips. Its success probability must then be an integer multiple of $2^{-m}$, which cannot equal $1/3$.

# When fewer flips suffice

A probability with a terminating binary expansion is called dyadic. For example,

$$\frac38=(0.011000\ldots)_{2}.$$

The first-heads procedure returns 0 on H, and 1 on TH or TTH. After TTT, every remaining digit is zero. There is no reason to keep flipping: return 0 immediately. The four leaves have masses $1/2$, $1/4$, $1/8$, and $1/8$; the two success leaves sum to $3/8$.

If $p=a/2^m$ is in lowest terms and $0<p<1$, this modification has cost

$$\mathbb E[T] =\sum_{k=0}^{m-1}2^{-k} =2-2^{1-m}.$$

Use the terminating expansion for a dyadic number. Although $(0.1000\ldots)_{2}=(0.0111\ldots)_{2}$, recognizing a constant tail is what allows the algorithm to finish without waiting for another heads.

# Why this cost is optimal

Imagine any exact algorithm as a binary tree: internal nodes request a fair flip, and leaves return 0 or 1. After $k$ flips, a path still in progress carries probability $2^{-k}$.

If $p$ is not dyadic, there must be an unfinished path at every finite depth. Otherwise, the algorithm would have the finite flip limit ruled out above. Consequently every exact sampler has

$$\Pr(T>k)\ge 2^{-k}, \qquad \mathbb E[T]\ge2.$$

The digit construction attains equality at every depth. For a reduced dyadic $a/2^m$, the same argument applies up to depth $m-1$, and the modified construction again attains every bound.

This is optimality for one exact sample starting with no saved randomness, using independent fair bits as the sole random source. It says nothing by itself about processor time or the cost of a batch of samples. Knuth and Yao’s discrete distribution generating trees extend this reasoning to any finite collection of outcomes; Saad and colleagues state their binary-digit characterization in Theorem 2.1. 2

Lumbroso’s Appendix B concludes that the optimal cost is two for every Bernoulli parameter. The dyadic exception matters: the derivation uses $\{x\}+\{-x\}=1$, which fails when $x$ is an integer. Here the braces mean fractional part, and $x=2^kp$ becomes an integer once a dyadic expansion terminates. The general tree result includes this exception.

# Exact arithmetic for rational probabilities

For $p=a/b$, generate the binary digits by doubling the numerator. The quotient is the next digit; the remainder is the numerator for the next step. For $1/3$, the remainders alternate between 2 and 1, producing digits 0 and 1 forever. Only the digits reached by the coin flips need to be computed.

def bernoulli(a, b, flip):
    """Return 1 with probability a/b; flip supplies independent fair bits."""
    if b <= 0 or not 0 <= a <= b:
        raise ValueError("require b > 0 and 0 <= a <= b")
    if a == b:
        return 1
    while a:
        digit, a = divmod(2 * a, b)
        if flip() == 1:
            return digit
    return 0

The inputs a and b are integers; flip() returns 0 or 1. The loop keeps the remainder below b, so its integer storage is proportional to $\log b$, even when the expansion repeats forever. Python’s integers avoid overflow in 2 * a; a fixed-width implementation must leave room for that doubling.

The terminating condition is important. Returning the current digit before its coin flip when the remainder becomes zero would be wrong: for $p=1/2$, it would always return 1. We must still distinguish the heads leaf from the tails branch.

For a simulation, a call might look like this:

import random

rng = random.Random(20150815)
result = bernoulli(1, 3, lambda: rng.getrandbits(1))

This is a reproducible pseudorandom simulation. The probability proof assumes independent fair bits. Also, requesting one bit from a generator does not necessarily consume only one bit of its internal output: efficient implementations can buffer a word and consume its bits across calls. Neither the proof nor this Python example establishes a speed advantage.

# What a floating-point comparison changes

The familiar random.random() < p uses a finite grid. Python’s random() supplies 53 bits of precision, with possible values $j/2^{53}$ for $0\le j<2^{53}$. For a threshold $0<p<1$, comparison against that grid accepts with probability $\lceil 2^{53}p\rceil/2^{53}$, assuming equiprobable grid points. An intended decimal or fraction may also have been rounded when converted to the floating-point threshold. 3

The discrepancy is tiny for many simulations. It is still a different contract from exact $1/3$, and the distinction becomes visible for very small probabilities: any positive threshold below $2^{-53}$ accepts the grid point zero, yielding probability $2^{-53}$.

Passing (1, 3) to the integer algorithm preserves the fraction. Passing the exact integer ratio of a stored float instead preserves that float’s dyadic value. Neither representation reconstructs an intended real number that was already rounded away.

# Weighted choices

Suppose three outcomes have probabilities $1/2$, $1/3$, and $1/6$. First accept A with probability $1/2$. If that fails, choose B with conditional probability

$$\frac{1/3}{1-1/2}=\frac23.$$

If B also fails, choose C. The unconditional probabilities are

$$\Pr(B)=\frac12\frac23=\frac13, \qquad \Pr(C)=\frac12\frac13=\frac16.$$

For nonnegative integer weights, no floating-point normalization is needed:

remaining = sum(weights)
for i, weight in enumerate(weights):
    if bernoulli(weight, remaining, flip):
        return i
    remaining -= weight

Require a nonempty list with at least one positive weight. Zero weights are skipped without using a flip; the last positive weight is selected with certainty if reached. In probability notation the $k$th trial uses $p_k/(1-\sum_{j<k}p_j)$, the probability remaining after all previous failures.

This extension is exact, but sequential choices are not generally optimal in bits. Ordering matters. Four equal categories cost $3/2+(3/4)2+(1/2)1=7/2$ flips on average under this procedure. Assigning the four outcomes to 00, 01, 10, and 11 uses exactly two.

# An irrational probability without its digits

There is another way to build exact probabilities: compose simpler coins. For example, we can make a coin with success probability $1/e$ without calculating $e$, taking a logarithm, or reading its binary expansion.

Imagine an initial success that costs no flip. Then try independent coins with success probabilities $1/2$, $1/3$, $1/4$, and so on. Stop at the first failure. Return 1 if the total number of successes, including the initial one, is even.

def inverse_e(flip):
    n = 2
    while bernoulli(1, n, flip):
        n += 1
    return int((n - 1) % 2 == 0)

If $N$ is that success count, reaching at least $n$ successes has probability

$$\Pr(N\ge n)=\frac12\frac13\cdots\frac1n=\frac1{n!}.$$

For $n=1$, the empty product is one. Therefore $\Pr(N=n)=1/n!-1/(n+1)!$, and the probability of an even count is

$$\Pr(N\text{ even}) =\frac1{2!}-\frac1{3!}+\frac1{4!}-\frac1{5!}+\cdots =e^{-1}.$$

The series also explains why the process terminates: its continuation probability decreases factorially. With our rational sampler, the expected number of fair flips is

$$\sum_{n=2}^{\infty}\frac{c_n}{(n-1)!} \approx2.35318,$$

where $c_n$ is the cost of the $1/n$ coin: two unless $n$ is a power of two, in which case the finite-tail saving applies. This is more than the two-flip digit sampler would use for $1/e$, but the only probabilities we need to represent are integer fractions.

Flajolet, Pelletier, and Soria study this broader idea under the name Buffon machines: finite procedures driven by fair coins that produce specified probabilities, including irrational ones. Their paper develops exponential-probability generators in Section 2.2. The rational-coin parity construction above is a composition whose proof is given here, rather than a transcription of their algorithm.5

# Sampling in batches and the compression limit

A coin with success probability $p$ has binary entropy

$$H(p)=-p\log_2p-(1-p)\log_2(1-p).$$

This is the average information in one output, measured in bits. It is at most one, and approaches zero for an almost deterministic coin. Why, then, does a single exact nondyadic sample still need two fair bits on average?

The decision tree must allocate whole leaves, each with probability $2^{-k}$. For one output, that indivisibility has a cost. Knuth and Yao’s bound, as stated in 2, places the optimum for a finite distribution between its entropy and its entropy plus two bits.

Reproduce this figure

Download source bundleextract, then runtar -xf coin-sampling-gen_figures.tar && cd coin-sampling-gen_figures/coin-sampling && uv run --locked gen_figures.py

Requires uv and Python 3.12; the locked packages download on first run.

Individual files:gen_figures.pylockfilefigure_quality.py

Now treat $r$ independent Bernoulli outputs as one outcome from a distribution on $2^r$ strings. Their total entropy is $rH(p)$. Applying the same theorem to that product distribution gives

$$rH(p)\le\mathbb E[L_r]<rH(p)+2.$$

The expected number of fair bits per output therefore lies between $H(p)$ and $H(p)+2/r$. For $p=1/3$, entropy is about $0.9183$ bits; with blocks of 100, the theorem guarantees a rate below $0.9383$ bits per output. Calling the one-coin procedure 100 times still costs 200 bits on average. Batching changes the sampler.

This is a theoretical bit-cost statement. Constructing the product distribution’s decision tree may be expensive, and the theorem is not a recommendation to materialize $2^{100}$ possibilities.

The connection to compression is more direct than a shared entropy formula. A compressor gives likely strings short descriptions. A sampler allocates more random-bit paths to likely strings. Arithmetic coding expresses both through intervals: encoding narrows an interval according to observed symbols; decoding random bits selects symbols according to the interval lengths. That exact interpretation assumes exact interval boundaries; finite-precision coding implementations may round them.4

The first-heads trick is a small instance of that correspondence. Each binary digit says whether its share of random-bit space belongs to success. Grouping outputs lets those shares be allocated more efficiently.

# Reproducing the calculations

The reproduction bundle contains the integer samplers and exact decision-tree checks. It explores every rational input $a/b$ with $1\le b\le64$ and $0\le a\le b$ through depth 12, comparing leaf masses and unfinished mass with the target probabilities and the optimal tail bound. These are exact rational calculations, so the checks do not depend on simulation noise. It also checks recycled dice and bounds the $1/e$ sampler’s success probability by accounting for every still-unfinished bit prefix. The figure source downloads include the dependency lockfile and shared plotting helper.

# References

[1] Jérémie Lumbroso (2013). Optimal Discrete Uniform Generation from Coin Flips, and Applications. Section 1 and Appendix B. ↩

[2] Feras A. Saad, Cameron E. Freer, Martin C. Rinard, and Vikash K. Mansinghka (2020). The Fast Loaded Dice Roller. Section 2, Theorem 2.1, states the Knuth–Yao result from The Complexity of Nonuniform Random Number Generation (1976). ↩

[3] Python documentation. random: Generate pseudo-random numbers. Precision and distribution of random(). ↩

[4] Paul G. Howard and Jeffrey Scott Vitter (1992). Practical Implementations of Arithmetic Coding. Arithmetic coding as a bridge between random bits and prescribed distributions. ↩

[5] Philippe Flajolet, Maryse Pelletier, and Michèle Soria (2011). On Buffon Machines and Numbers. SODA 2011, pp. 172–183, especially Section 2.2. ↩