arXiv:1507.00843v2 [math.PR] 27 Sep 2016

Optimal linear Bernoulli factories for small
mean problems∗
Mark Huber
mhuber@cmc.edu
Version: September 29, 2016

Abstract
Suppose a coin with unknown probability p of heads can be flipped
as often as desired. A Bernoulli factory for a function f is an algorithm that uses flips of the coin together with auxiliary randomness to
flip a single coin with probability f (p) of heads. Applications include
perfect sampling from the stationary distribution of certain regenerative processes. When f is analytic, the problem can be reduced to
a Bernoulli factory of the form f (p) = Cp for constant C. Presented
here is a new algorithm that for small values of Cp, requires roughly
only C coin flips. From information theoretic considerations, this is
also conjectured to be (to first order) the minimum number of flips
needed by any such algorithm.
For large values of Cp, the new algorithm can also be used to build
a new Bernoulli factory that uses only 80% of the expected coin flips
of the older method. In addition, the new method also applies to the
more general problem of a linear multivariate Bernoulli factory, where
there are k coins, the kth coin has unknown probability pk of heads,
and the goal is to simulate a coin flip with probability C1 p1 +· · ·+Ck pk
of heads.

Keywords: randomized algorithm, near perfect simulation, regenerative
processes
MSC Classes: 65C50, 68Q17
The
final
publication
is
http://dx.doi.org/10.1007/s11009-016-9518-3.
∗

1

available

at

Springer

via

1

Introduction

The notion of a Bernoulli factory was introduced in Asmussen, Glynn, and Thorisson
(1992) in the context of generating samples exactly from the stationary distribution of a regenerative Markov process. A Bernoulli factory works as
follows. Suppose we have the ability to draw independent identically distributed (iid) Bernoulli random variables, each of which is 1 with probability
p and 0 with probability 1 − p (write X ∼ Bern(p).) Then given a function
f , the goal is to use a random number of draws from X to build a new random variable which is also Bernoulli, but with chance f (p) of being 1 for a
speciﬁed function f . In Asmussen, Glynn, and Thorisson (1992), the needed
function was a linear function, namely a constant times p. This simple case
generalizes: in Nacu and Peres (2005) it was shown that the ability to draw
from f (p) = 2p could be used to build a Bernoulli factory for any analytic f
that was bounded away from 1.
The focus here is on building a nearly optimal linear Bernoulli factory
where Cp is known to be small. This is nearly optimal in the sense that
it uses (to ﬁrst order) only C ﬂips of the coin, and there is strong evidence
to indicate that at least C ﬂips are necessary. As in Huber (to appear), a
Bernoulli factory can be deﬁned as follows.
Definition 1. Given p∗ ∈ (0, 1] and a function f : [0, p∗ ] → [0, 1], let A
be a computable function that takes as input a number u ∈ [0, 1] together
with a sequence of values in {0, 1}, and returns an output in {0, 1}. For any
p ∈ [0, p∗], X1 , X2 , . . . iid Bern(p), and U ∼ Unif([0, 1]), let T be the inﬁmum
of times t such that the value of A(U, X1, X2 , . . .) only depends on the values
of X1 , . . . , Xt . If the following holds, then call A a Bernoulli factory.
1. T is a stopping time with respect to the natural ﬁltration that is ﬁnite
with probability 1.
2. A(U, X1, X2 , . . .) ∼ Bern(f (p)).
Call T the running time of the Bernoulli factory.
Colloquially, a draw X ∼ Bern(p) will be refereed to as a coin ﬂip, or
more speciﬁcally, a p-coin ﬂip. The result X = 1 corresponds to heads on
the coin, while X = 0 indicates tails. So a Bernoulli factory attempts to ﬂip
a coin with f (p) chance of heads, by using a random number of coin ﬂips
from the original coin together with some auxiliary randomness.
2

Asmussen, Glynn, and Thorisson (1992) introduced Bernoulli factories
for an application in perfect simulation, but did not show that they exist.
Keane and O’Brien (1994) constructed the ﬁrst general Bernoulli factories,
showing that such a factory with ﬁnite running time existed if and only if
f (p) was continuous over [0, p∗ ] for some p∗ ∈ (0, 1], and either it holds that
f (p) is identically 0 or 1, or that both f (p) and 1 − f (p) are polynomially
bounded away from 0 and 1 over the allowable range of p.
Their strategy for building a Bernoulli factory was to construct Bernstein
polynomials that approximated the function f (p) as closely as possible. Bernstein polynomials are linear combinations of functions of the form pk (1−p)n−k
where both n and k ≤ n are nonnegative integers. Such polynomials can be
created from the coin by ﬂipping it n times and seeing if exactly k heads
and n − k tails result. Keane and O’Brien could show that the running time
T was ﬁnite with probability 1 for their algorithm, but not much more. In
particular, they could not show any bounds on the average running time, or
even that it was ﬁnite.
Nacu and Peres (2005) developed this approach further, and showed that
Bernstein polynomials could be constructed tightly enough that the running
time would have a ﬁnite expectation. In addition, they showed that the tail
of the distribution of their running time declined exponentially. Moreover,
their work contained a proof that f (p) = 2p is in a sense the most important
function, since it can be used to construct a Bernoulli factory for any function
that is both real analytic over [0, 1] and bounded away from 1.
However, their approach was not a practical algorithm. While the number of coin ﬂips had ﬁnite expectation, the amount of memory and time
needed to compute the function A grew exponentially with the number of
ﬂips. Work of Latuszyński, Kosmidis, Papspiliopoulos, and Roberts (2011)
solved this issue, and gave the ﬁrst practical implementation of the Nacu
and Peres approach. Their approach created a pair of reverse time processes,
one a supermartingale, the other a submartingale, that converged on the target f (p). The values could be computed without the exponential overhead
associated with the Nacu-Peres algorithm.
To bound f (p) = Cp away from 1, they considered the function f (p) =
min{Cp, 1 − ǫ} so that the function was deﬁned over the entirety of [0, 1].
However, this was not strictly necessary, as in the original application of
Asmussen, Glynn, and Thorisson (1992), it was possible to easily insure that
f (p) ≤ 1−ǫ. By not trying to sample from the function f (p) = min{Cp, 1−ǫ}
for all values of p, but only for those with Cp ≤ 1−ǫ, a new approach became
3

possible.
This new approach was developed in Huber (to appear), and was the ﬁrst
that did not begin with Bernstein polynomial approximations. Instead, this
approach used ﬂips of the coin to alter the problem in ways that insured that
the ﬁnal output had the correct distribution. For instance, suppose the goal
was to generate a coin with probability of heads 2p. Then ﬂip the original
coin once. If the coin is heads, the the output is heads. Otherwise, it is
necessary to ﬂip a p/(1 − p)-coin.
That way, the chance the ﬁnal output is heads is p(1)+(1−p)·p/(1−p) =
2p. By advancing carefully in this manner, it was shown how to build a
Bernoulli factory such that for Cp ≤ 1 − ǫ where ǫ is a known constant,
E[T ] ≤ 9.5Cǫ−1 .
Moreover, the same work showed that this running time is the best possible up to a constant. Speciﬁcally, in Huber (to appear) it was shown that that
any Bernoulli factory (to ﬁrst order) must use on average at least 0.04Cǫ−1
coin ﬂips. It remains an open question what the best constants for the
lower and upper bounds are, although there is strong reason to believe (see
Section 4) that on average at least C ﬂips (to ﬁrst order) are necessary to
generate a Cp coin.
The primary application of the Bernoulli factory, starting with Asmussen,
Glynn, and Thorisson (1992), is to generate perfect samples from the stationary distribution of regenerative Markov chains. In Lee, Doucet and Latusyzński
(2014), it was shown how to use the Bernoulli factory in Huber (to appear)
to generate coins where Cp ≤ 2/3. Under this condition, the older algorithm
gave a bound of 38 coin ﬂips on the average number needed.
Without going into the details of their algorithm, by tripling the expected
running time of there algorithm, it is possible to ensure that Cp ≤ 2/9.
Under these conditions, the expected number of ﬂips for the new algorithm
is bounded above by 7.8 (see Theorem 1) giving an algorithm that only
requires 62% as many work on average as the old one.
This work gives the following results.
1. For Cp small, an algorithm will be given that uses only C coin ﬂips on
average.
2. For Cp at most 1 − ǫ for known ǫ, an algorithm will be given that uses
only 7.57Cǫ−1 coin ﬂips on average.
4

3. The new algorithm can be extended from the function Cp for single
variate coins to the multivariate coin problem where there are k coins
with unknown means p1 , . . . , pk . Suppose the goal is to generate Bern(r)
where
r(p1 , . . . , pk ) = C1 p1 + · · · + Ck pk .
(1)
Then setting C = C1 + · · · + Ck , the algorithm in the multivariate case
has running time equal to the single coin case.
More precisely, the running time of the new algorithm is given as follows.
Theorem 1. Suppose it is known that Cp ≤ M for a constant M < 1/2.
Then there exists an algorithm for producing a Cp-coin that uses on average
at most


C
15.2
+ Cp · C
(1 − 2M)(1 + Cp)
1 − 2M + Cp
coin flips.
This theorem is shown in Section 3.5. The multivariate version is similar,
and is shown in Section 5.
Theorem 2. Let C = C1 + · · · + Ck , and r = C1 p1 + · · · + Ck pk , Suppose it is
known that r ≤ M for a constant M < 1/2. Then there exists an algorithm
for producing an r-coin that uses on average at most


15.2
C
+ Cp · C
(1 − 2M)(1 + r)
1 − 2M + r
flips from among the k coins.
The remainder of this paper is organized as follows. Section 2 presents the
algorithm for small r, and shows correctness and the bound on the running
time. Section 3 gives the extension to the multivariate problem and for larger
values of r, and also includes the proofs of correctness and the bound on the
running time. Finally, Section 4 considers why C ﬂips is likely the best
possible.

5

2

The algorithm for small Cp

Let r = Cp. The ﬁrst piece of the algorithm is a method for drawing from
the logistic Bernoulli factory
f (p) =

r
1+r

that uses T coins, where E[T ] = C/(1 + r).
As usual, say that X is exponential with rate λ (write X ∼ Exp(λ)) if X
has density fX (s) = exp(−λs)1(s ≥ 0). Here 1(·) is the indicator function
that evaluates to 1 when the argument is true and 0 when the argument is
false. The following basic facts about exponentials will prove useful.
Fact 1. Let X ∼ Exp(λ1 ) and Y ∼ Exp(λ2 ) be independent. Then P(X ≤
Y ) = λ1 /(λ1 + λ2 ).
Fact 2 (Memoryless). If X ∼ Exp(λ), then for s > 0, the conditional distribution of X − s given X > s is exponential with rate λ as well. That is,
[X − s|X > s] ∼ Exp(λ).
Exponentials can be employed to deﬁne a one dimensional Poisson point
process.
Definition 2. Let A1 , A2 , . . . be independent and identically distributed (iid)
exponential random variables with rate λ. Then
P = {A1 , A1 + A2 , A1 + A2 + A3 , . . .}
forms a Poisson point process on [0, ∞) of rate λ. For [a, b] ⊂ [0, ∞), P ∩[a, b]
is a Poisson point process on [a, b] of rate λ.
Several well known facts about Poisson point processes are useful.
Fact 3. The converse of the definition holds: any Poisson point process
P ⊂ [0, ∞) of rate λ with points 0 < P1 < P2 < · · · has P1 ∼ Exp(λ) and
Pi − Pi−1 ∼ Exp(λ), and all these exponentials are independent.
Fact 4. Let P = {P1 , P2 , . . .} be a Poisson point process. Let B1 , B2 , . . .
be a sequence of iid Bern(p) random variables. Then P ′ = {Pi : Bi = 1}
is a Poisson point process of rate λp. [The process P ′ is called the thinned
process.]
6

Fact 5. Let P1 and P2 be independent Poisson point processes of rate λ1 and
λ2 over [0, ∞). Then P1 ∪ P2 is a Poisson point process of rate λ1 + λ2 over
[0, ∞).
Fact 6. The expected number of points in a Poisson point process of rate λ
over [a, b] is Poisson distributed with mean λ(b − a).
These ideas can be used to build the logistic Bernoulli factory for r/(1+r).
Logistic Bernoulli Factory
1)
2)
3)
4)
5)
6)

Input: C

X ← 0, draw A ← Exp(1)
Draw T ← Exp(C)
While X = 0 and T < A
Draw B ← Bern(p)
If B = 1 then X = 1, else T ← T + Exp(C)
Return X

Note that line 4 can be accomplished in constant time (with Θ(k) preprocessing time) using the Alias method of Walker (1974).
Lemma 1. The output of Logistic Bernoulli Factory is a Bernoulli with
mean r/(1 + r).
Proof. Let T1 , T2 , . . . be the successive values of T taken on in the algorithm,
and B1 , B2 , . . . the successive values of B. Since Ti+1 − Ti is Exp(C), the {Ti }
form a Poisson point process P of rate C. Let P ′ be the points Ti ∈ P with
Bi = 1. Then P ′ is a point process with rate Cp = r.
Let T1′ = min{P ′}. The while loop examines the P ′ process, and returns
1 if T1′ < A, and 0 otherwise. By Fact 3, T1′ ∼ Exp(r), and A ∼ Exp(1), so
P(T1′ < A) = r/(1 + r) by Fact 1.
Lemma 2. In one call to Logistic Bernoulli Factory, the expected number of coin flips needed is C/(1 + r).
Proof. The Poisson process of rate C combined with the process of rate 1
forms a Poisson point process of rate C + 1. The chance that any point of
this process is from the thinned rate Cp process combined with the rate 1
process is (Cp + 1)/(C + 1). Therefore, the number of points generated in the
rate C + 1 process has a geometric distribution with mean (C + 1)/(Cp + 1).
7

Each of these points has a C/(C + 1) chance of coming from the rate C
process initially, and so requires a coin ﬂip. Therefore, combining these eﬀects
gives an expected number of coin ﬂips of [C/(C + 1)][(C + 1)/(Cp + 1)] =
C/(r + 1).
Now suppose that there is a known M < 1/2 such that r ≤ M. Let
BF(C) denote the Bernoulli Factory from Huber (to appear) that ﬂips a Cp
coin using on average 9.5C(1 − M)−1 ﬂips of the original coin. Consider the
following algorithm.
Small r 1D Bernoulli Factory
1)
2)
3)
4)
5)
6)

Input: C, M

β ← 1/(1 − 2M)
Draw Y ← Logistic Bernoulli Factory(βC)
Draw B ← Bern(1/β)
If Y = 0, then X ← 0
Elseif Y = 1 and B = 1, then X ← 1
Else X ← BF(βC/(β − 1))

Lemma 3. Algorithm Small r 1D Bernoulli Factory produces a Bernoulli
distributed output with mean Cp ≤ M < 1/2, and requires at most (on
average)


C
1
+ Cp · 19C
(1 − 2M)(1 + Cp)
1 − 2M + Cp
coin flips to do so.
Note that for small p and M, this running time is to ﬁrst order just C.
Proof. First show correctness. Let A1 be the event that Y = 1 and B = 1 in
the algorithm (in which case line 5 sets X to be 1), and A2 be the event that
Y = 1, B = 0, and a call to BF(Cβ/(β − 1)) returns a 1 (in which case line 7
sets X to be 1). These are disjoint events, and the output of the algorithm
is X = 1(A1 ) + 1(A2 ). Therefore,
P(X = 1) = P(A1 ) + P(A2 ).
The value of Y is the call to Logistic Bernoulli Factory(βC), and
so P(Y = 1) = βCp/(1 + βCp). For the Bernoulli B, P(B = 1) = 1/β.
Therefore P(A) = P(Y = 1)P(B = 1) = Cp/(1 + βCp).
8

The output of BF(Cβ/(β − 1)) is 1 with probability equal to βCp, so
P(A2 ) =

βCp
βCp
βCp
(1 − 1/β)
= Cp
1 + βCp
β−1
1 + βCp

Therefore
P(X = 1) = Cp

1
βCp
+ Cp
= Cp
1 + βCp
1 + βCp

as desired.
Now for the running time. Lemma 2 gives a running time of βC/(1 + Cp)
for line 1. Line 5 is executed with probability (β − 1)Cp/(1 + βCp). By the
way β was chosen, Cpβ/(β − 1) ≤ 1/2. Therefore, Theorem 1.1 of Huber
(to appear) gives that the call to BF(Cβ/(β −1)) requires at most 19[Cβ/(β −
1)] ﬂips. Therefore, the total number of ﬂips is on average at most


19C
C
.
+ Cp ·
(1 − 2M)(1 + Cp)
1 − 2M + Cp

This is close to Theorem 1, but the constant of 19 in the second term is
larger. To improve this algorithm, and eliminate the need for the call to the
old BF algorithm, it is necessary to consider what happens for larger r.

3

Large r algorithm

In this section, the algorithm of the previous section is improved to allow for
all r ∈ [0, 1 − ǫ], where ǫ is arbitrarily close to 0. Along the way, the older
9.5Cǫ−1 algorithm of Huber (to appear) is improved to a 7.5Cǫ−1 algorithm.
The ﬁrst step is to build a random coin ﬂip whose mean is slightly larger
than r. If this coin is tails, return tails for r. If the coin returns heads, heads
will be returned with probability close to 1. Otherwise, a new coin will need
to be ﬂipped.
3.1

A coin flip with mean slightly larger than r

Consider an asymmetric random walk on the integers Ω = {0, 1, . . . , m},
where given the current state Xt , the next state is either max{0, Xt − 1}, or
min{Xt + 1, m}. The transition probabilities are
P(Xt+1 = min{i + 1, m}|Xt = i) = pr , P(Xt+1 = max{i − 1, 0}|Xt = i) = qr ,
9

where pr + qr = 1. This is also called the Gambler’s Ruin walk.
The following facts about this well known process will be helpful.
Fact 7. Suppose pr 6= qr and T = inf{t : Xt ∈ {0, m}}. Then
1 − (qr /pr )X0
1 − (qr /pr )m
X0
m
E[T ] =
−
· P(XT = m).
qr − pr qr − pr

P(XT = m) =

(2)
(3)

Fact 8. Suppose X0 = m, pr < qr , and T = inf{t : Xt = 0}. Then
E[T ] ≤ m/(qr − pr ).
Consider the following Bernoulli factory that begins a Gambler’s Ruin
walk starting at state 1, and returns heads if the state reaches 0 before it
reaches m.
A

Input: m, C

1) s ← 1
2) While s ∈ {1, 2, . . . , m − 1}
3)
B ← Logistic Bernoulli Factory(C)
4)
s ← s − 2B + 1
5) Return 1(s = 0)
Lemma 4. The output of A is a Bernoulli with mean r(1 − r m−1 )/(1 − r m ).
Proof. Logistic Bernoulli Factory outputs a Bernoulli that has mean
r/(1 + r). Hence pr = 1/(1 + r), qr = r/(1 + r), and qr /pr = r. From Fact 7,
in line 5 that makes P(I = 0) = 1 − (1 − r)/(1 − r m ) = (r − r m )/(1 − r m ) =
r(1 − r m−1 )/(1 − r m ).
Lemma 5. The expected number of coin flips used by A is at most C(m − 1).
Proof. As in the last proof pr = 1/(1 + r) and qr = r/(1 + r). So
qr − pr =

r
1
1−r
−
=−
,
1+r 1+r
1+r

and (qr /pr ) = r. Using (3),




1−r
m
1
1+r
.
· m
− 1 = (1 + r)
−
E[T ] =
1−r
1 − rm
1 − rm 1 − r
10

Let f (r) = m/(1 − r m ) − 1/(1 − r). Then it holds that f (r) < m − 1 for
all r ∈ (0, 1) and m ≥ 1. Note that for all r ∈ (0, 1):
(1 − r m )
< (m − 1)(1 − r m )
1−r
⇔ m − mr − 1 − r − · · · − r m−1 < m − 1 − r m (m − 1)
⇔ (m − 1)r m < r + r 2 + · · · + r m−1 + mr
⇔ m − 1 < r 1−m + r 2−m + · · · + r −1 + mr 1−m .

f (r) < m − 1 ⇔ m(1 − r) −

Since r ∈ (0, 1), r i−m ≥ 1 for all i ∈ {1, . . . , m − 1}, so the right hand side is
strictly greater than the left hand side. Note f (0) = m−1, so for r ∈ [0, 1−ǫ],
the function is at most m − 1.
Each call to line 3 requires on average C/(1 + r) time by Lemma 2.
Together, the overall number of steps (on average) is at most
C
(1 + r)(m − 1)
= C(m − 1).
1+r

3.2

After the flip

Here is how A can be useful. Using A, it is possible to generate a Bernoulli
random variable that is 1 with probability
pβ = βr

1 − (βr)m−1
1 − (βr)m

for any constant β > 1. By choosing β large enough, pβ ≥ r. Note that
pβ /β ≤ r. So r ∈ [pβ /β, pβ ].
So the algorithm works as follows. First ﬂip a pβ -coin. If it is tails, then
return tails for the r-coin as well. If it is heads, then ﬂip a (1/β)-coin. If
that is heads as well, return heads for the r-coin. Otherwise, ﬂip a p′ -coin,
and return the same value for the r-coin.
For this algorithm to work, p′ must satisfy:
pβ
r=
+ p′ pβ (1 − (1/β)).
β
Solving for p′ gives


1
(βr)m−1
′
p =
,
β − 1 1 + (βr)1 + · · · + (βr)m−2
so the next step of the algorithm is ﬁguring out how to generate a p′ -coin.
11

3.3

Generating a p′ -coin

Fortunately, we do not have to actually generate a p′ -coin for all possible
values of β, as we are allowed to choose the value of β to use, as long as
pβ ≥ r for our choice of β. Let
β =1+

1
,
m−1

so (β − 1)−1 = m − 1. Then
p′ =

(m − 1)(βr)m−1
.
1 + (βr) + · · · + (βr)m−2

(4)

Note that p′ ≤ 1 which gives that pβ ≥ r for this choice of β.
The algorithm for generating a p′ -coin for p′ as in (4) will be called B here,
and is shown graphically in Figure 1. Notice that if the ﬁrst ﬂip is heads and
the second ﬂip is tails, then our problem has changed to the same problem,
but with m reduced to m − 1.
Flip a
βr-coin

Heads

Flip a
(βr)m−2
1 + · · · + (βr)m−2
coin

Tails
Return
tails

Heads

Return
heads

Tails

Flip an
(m − 2)(βr)m−2
1 + · · · + (βr)m−3
coin

Figure 1: A graphical illustration of Algorithm B.
To utilize this procedure, it is necessary to be able to generate a (βr)m−2 /(1+
· · ·+ (βr)m−2 coin. Fortunately, this can be accomplished fairly quickly using
the Gambler’s ruin chain from earlier.

12

High Power Logistic BF
Input: m, β, C Output: X ∼ Bern((βr)m /(1 + · · · + (βr)m ))
1) s ← 1
2) While s ∈ {1, . . . , m}
3)
Draw B ← Logistic Bernoulli Factory(βC)
4)
s ← s + 2B − 1
5) X ← 1(s = m + 1)
Lemma 6. The output of High Power Logistic BF has distribution Bern((βr)m/(1+
· · ·+ (βr)m)). The expected number of coin flips used is at most βC/(1 − βr).
Proof. This is a Gambler’s ruin where p = βr/(1 + βr) and q = 1/(1 + βr),
so q/p = 1/(βr). Hence from Fact 7,
P(s = m + 1) =

(βr)m(1 − βr)
(βr)m
1 − (1/(βr))1
=
=
.
1 − (1/(βr))m+1
1 − (βr)m+1
1 + · · · + (βr)m

Also from Fact 7, note q − p = (1 + βr)/(1 − βr), so if T is the number
of times line 3 is called,


1 + βr
(βr)m (1 − βr)
1 + βr
E[T ] =
1 − (m + 1)
≤
.
m+1
1 − βr
1 − (βr)
1 − βr
Each call to Logistic Bernoulli Factory takes time βC/(1 + βr), so
the overall number of coin ﬂips (on average) is at most βC/(1 − βr).
In pseudocode, algorithm B looks like this.
B

Input: ǫ, m, β, C

Output: X

1)
2)
3)
4)
5)
6)
7)
8)

X ← 0.5
While X ∈
/ {0, 1}
Draw B1 ← Linear Bernoulli Factory(1 − (1 − ǫ)β, β · C)
If B1 = 0 then X ← 0
Else
B2 ← High Power Logistic BF(m − 2, β, C)
If B2 = 1 then X ← 1
Else m ← m − 1

13

The most important thing to note here is that like many perfect simulation algorithms, this method employs recursion. We do not yet have an
algorithm for completing line 3! However, this algorithm B can be used as
a subroutine to create such an algorithm, and then this subroutine will call
the ﬁnished algorithm.
Lemma 7. The output of B has distribution
Bern((m − 1)(βr)m−1 /(1 + · · · + (βr)m−2 )).
Proof. The proof is by induction. When m = 2, if B1 = 1 then B2 ∼ Bern(1),
so X = 1 with probability βr as desired.
Now suppose that the result holds for m, consider m + 1. Then


(βr)m−2
1 + · · · + (βr)m−3 (m − 2)(βr)m−2
P(X = 1) = βr
+
·
1 + · · · + (βr)m−2 1 + · · · + (βr)m−2 1 + · · · + (βr)m−3
(m − 1)(βr)m−1
,
=
1 + · · · + (βr)m−2
completing the induction.
3.4

The new linear Bernoulli Factory

With these preliminaries in place, the overall algorithm is as follows.
Linear Bernoulli Factory
1)
2)
3)
4)
5)
6)
7)
8)

Input: ǫ, C

Output: B

m ← ⌈4.5ǫ−1 ⌉ + 1, β ← 1 + 1/(m − 1)
B1 ← A(m, β · C)
If B1 = 1
Draw B2 ← Bern(1/β)
If B2 = 1 then B ← 1
Else
Draw B ← B(m, β, C)
Else B ← 0

Now consider the expected number of coin ﬂips used by the algorithm.
As will become clear in the proofs of Lemmas 8 and 10 below, making m =
Θ(ǫ−1 ) is the correct choice. That leaves the choice of constant up to us, and
14

the constant of 4.5 from Line 1 was chosen to make the overall running time
as small as possible.
This algorithm calls A and B. Line 2 of B needs to draw a Bern(βr) random
variable. The best way to draw these random variables is to recursively call
Linear Bernoulli Factory. In order to ensure that this back in forth calling
eventually comes to a halt with probability 1, it is easiest to bound the total
expected number of calls to Linear Bernoulli Factory.
Lemma 8. The expected number of calls to Linear Bernoulli Factory is
at most 1.4.
Proof. Let m1 be the value of m in the ﬁrst call to Linear Bernoulli Factory,
and β1 = 1 + 1/(m1 − 1). From this ﬁrst call there is a chance of calling B,
which in turn calls Linear Bernoulli Factory with m2 and β2 . Each of
those second generation calls might call a third generation, and so on. To
bound the expected number of calls to Linear Bernoulli Factory sum over
all possible calls of the probability that that call is executed. Let Ni denote
the number of ith generation calls.
The expected number of calls in the ﬁrst generation is 1. Consider a call
in the second generation. In order for that call to be made, there must have
been a call to B from the ﬁrst generation, and all prior second generation calls
from line 3 of B must have had B1 = 0. The number of times the while loop
in B is executed is stochastically dominated by a geometric random variable
with mean 1/(1 − β1 r). Since
1 − β1 r ≥ 1 − (1 + 1/⌈4.5ǫ−1 ⌉)(1 − ǫ)
= (7/9)ǫ + (2/9)ǫ2 ,

(5)
(6)

the number of calls made is bounded (in expectation) by (9/7)ǫ−1.
But before B is even called, ﬁrst it must have held that B1 = 1 and B2 = 0
in lines 2 and 4 of the ﬁrst generation call to Linear Bernoulli Factory.
The probability that a call to B is made is at most
(1 − 1/β1 )βr(1 − (β1 r)m−1 )/(1 − (β1 r)m ) ≤

1
m−1

So the expected number of calls to Linear Bernoulli Factory in the
second generation is bounded by
E[N2 |N1 ] ≤ N1 [1/(m1 − 1)](9/7)ǫ−1 ≤ (2/7)N1 .
15

i
This step forms
Pthe basis of an induction that gives E[Ni ] ≤ (2/7) N1 =
i
(2/7) . Therefore
E[Ni ] ≤ 1/(1 − 2/7) = 1.4.

Lemma 9. The output B of Linear Bernoulli Factory has B ∼ Bern(r).

Proof. Line 2 of B requires a draw B1 ← Bern(βr). Suppose that for the
ﬁrst L times this line is called, the Linear Bernoulli Factory is called to
generate this random variable. Then, from the L + 1st time onwards, an
oracle generates the random variable.
We show by strong induction that Linear Bernoulli Factory generates
from Bern(r) for any ﬁnite M. The base case when M = 0 operates as follows.
Lemma 7 immediately gives in this case that a call to B returns a random
variable with distribution Bern((m−1)(βr)m−1 /(1+· · ·+(βr)m−2)). Lemma 4
gives that B1 from line 2 has distribution Bern((βr)m/(1 + · · · + (βr)m )).
Putting this together gives




1 − (βr)m−1 1
1
(m − 1)(βr)m−1
P(B = 1) = (βr)
+ 1−
1 − (βr)m β
β 1 + · · · + (βr)m−2


1
(βr)m−1
1 − (βr)m−1 1
+
= (βr)
1 − (βr)m β β 1 + · · · + (βr)m−2
(1 − (βr)m−1 ) 1 + · · · + (βr)m−1
·
=r·
1 − (βr)m
1 + · · · + (βr)m−2
= r.
This is the rare induction proof where the base case is just as hard as
the induction step. Suppose it holds for L, and consider what happens for
call limit L + 1. Then the ﬁrst call to Linear Bernoulli Factory might call
B, which might call Linear Bernoulli Factory. But the ﬁrst such call has
used up one call, so only has L + 1 − 1 calls remaining, so by strong induction
each returns the correct distribution. Hence Lemma 7 holds, and the ﬁrst
call returns the correct distribution by the same argument as the base case.
Let N be the random number of calls to Linear Bernoulli Factory
needed by the algorithm. Then let B be the output when N is unbounded,
and BM be the output when a limit on calls equal to L is in place. Then
P(B = 1) = P(B = 1, N ≤ L) + P(B = 1, N > L)
= P(BL = 1, N ≤ M) + P(B = 1, N > L)
= P(BL = 1) − P(BL = 1, N > L) + P(B = 1, N > L).
16

Both P(BL , N > L) and P(B = 1, N > L) are bounded above by P(N > L).
Since by the last lemma E[N] ≤ 1.4, limL→∞ P(N > L) = 0. The only way
this can hold for all L is if P(B = 1) = P(BL = 1) for all L, so B has the
correct distribution.
Lemma 10. Linear Bernoulli Factory uses on average at most 7.67Cǫ−1
coin flips to generate B ∼ Bern(r).
Proof. From the proof of Lemma 8, the expected number of calls to the ith
generation of Linear Bernoulli Factory is bounded above by (2/7)i .
From (6), at each successive generation of calls, ǫ is being multiplied by a
factor of at least 7/9. So an ith generation call to Linear Bernoulli Factory
has an m value of at most ⌈4.5(9/7)i ǫ−1 ⌉ + 1, where ǫ was the input for the
0th generation.
Coin ﬂips occur during the call to A, and Lemma 4 bounds the expected number of coin ﬂips by C⌈4.5(9/7)i ǫ−1 ⌉. So the expected total ﬂips
coming from the ith generation of Linear Bernoulli Factory is at most
C[4.5(18/49)iǫ−1 + (2/7)i ].
Now look at the ﬂips coming from an ith generation call to B. This generation is only called from an ith generation call to Linear Bernoulli Factory,
of which there the expected number is at most (2/7)i . The call to B occurs
with probability at most (β − 1)r, so at most r/m. The while loop inside
B is run (on average) at most ǫ−1 times, each of which could make a call
to High Power Logistic BF. By Lemma 6 this requires at most βCǫ−1 coin
ﬂips. So the total number of coin ﬂips from an ith generation call to B is at
most
(2/7)i [r/m][βCǫ−1 (9/7)i+1 ]ǫ−1 ≤ (9/7)4.5−1(18/49)i Cǫ−1 .
Summing over these ﬂips and the ones from Linear Bernoulli Factory
gives a total sum of
(469/62)Cǫ−1 ≤ 7.57Cǫ−1
coin ﬂips on average.
3.5

Small r

Now that a recursive algorithm for large r has been built, a recursive analogue
for Small r 1D Bernoulli Factory works as follows.
17

Small r Bernoulli Factory
1)
2)
3)
4)
5)
6)

Input: C, M

β ← 1/(1 − 2M)
Draw Y ← Logistic Bernoulli Factory(βC)
Draw B ← Bern(1/β)
If Y = 0, then X ← 0
Elseif Y = 1 and B = 1, then X ← 1
Else X ← Linear Bernoulli Factory(Cβ/(β − 1))

Lemma 11. Algorithm Small r Bernoulli Factory produces a Bernoulli
distributed output with mean Cp ≤ M < 1/2, and requires at most (on
average)


C
15.2
+ Cp · C
(1 − 2M)(1 + Cp)
1 − 2M + Cp
coin flips to do so.
Proof. The proof is essentially the same as that of Lemma 3.

4

Lower bound

To see why it is unlikely that a method that uses fewer than Cǫ−1 coin ﬂips
can be constructed, consider building an unbiased estimate of p.
The standard estimate is to generate X1 , . . . , Xn iid Bern(p), and then
use the sample average p̂n = (X1 + · · · + Xn )/n as an unbiased estimate of
p. This estimate is unbiased, and has variance p(1 − p)/n.
Now consider the estimate Y /C, where Y ∼ Bern(Cp). Then E[Y /C] =
Cp/C = p so this estimate is also unbiased, and the variance is Cp(1 −
Cp)/C 2 = p(1 − Cp)/C. Therefore, this estimate that used one draw from
Bern(Cp) has the variance of the estimate that used n = C(1 − p)/(1 − Cp)
draws from the p-coin.
The Cramér-Rao lower bound (see, for instance Bickel and Doksum (1977))
on the variance of an unbiased estimate of p is
p(1 − p)
.
n
That is, any unbiased estimate that uses up to n ﬂips of the p-coin must have
variance at least p(1 − p)/n. That immediately gives that any algorithm for
generating a Cp-coin that uses a deterministic number n of coin ﬂips must
18

have n ≥ C(1 − p)/(1 − Cp). Of course, this does not quite apply to a
Bernoulli Factory, because here a random number of coin ﬂips is used.
However, it is strong evidence that C(1 − p)/(1 − Cp) is a lower bound
on the expected number of coin ﬂips needed by an algorithm.

5

Multivariate Bernoulli Factory

This new algorithm was designed for the single coin problem, in this section
consider generating a coin ﬂip whose probability of heads is the sum of the
probability of heads on two diﬀerent coins each of which has an unknown
probability of heads. Unlike the single coin ﬂip, there is no immediate application, however, it is useful to know that the single coin algorithm can be
easily generalized to solve this problem should the need arise.
More generally, the goal is now to generate a coin ﬂip with probability
r = C 1 p1 + · · · + C k pk
of heads, where r is bounded away from 1, using as few ﬂips of the coins
as possible. When k = 1, this is the linear Bernoulli factory studied in the
previous sections. Formally, a multivariate Bernoulli factory is deﬁned as
follows.
Definition 3. Given a computable function f : [0, 1]k → [0, 1], a multivariate
Bernoulli factory is a computable function
A : [0, 1] × {0, 1}{1,2,...}

k

→ {0, 1},

such that if U ∼ Unif([0, 1]) and the Xi,j are independent random variables
with (∀i ∈ {1, . . . , k})(∀j ∈ {1, 2, . . .})(Xi,j ∼ Xi ), then the following properties hold.
1. There exist random variables (T1 , . . . , Tk ) ∈ {1, 2, . . .}k such that the
∞
value of A(U, {X1,i }∞
i=1 , . . . , {Xk,i }i=1 ) only depends on the values of
1
k
{X1,i }Ti=1
, . . . , {Xk,i }Ti=1
, and for all (t1 , . . . , tk ), the event (T1 , . . . , Tk ) =
1
k
(t1 , . . . , tk ) is measurable with respect to {X1,i }ti=1
, . . . , {Xk,i }ti=1
.
∞
2. A(U, {X1,i }∞
i=1 , . . . , {Xk,i }i=1 ) ∼ Bern(f (p1 , p2 , . . . , pk )).

Call T1 + · · · + Tk the running time of the algorithm.
19

The key to the single coin algorithm was generating a random variable
that was exponential with rate parameter r. In the single coin case, this was
done by generating a Poisson process of rate C, then thinning.
For the multivariate coin case, consider generating k independent Poisson
point processes P1 , . . . , Pk , where Pi has rate Ci . Then thin each process Pi
with coin i to obtain a Poisson point process of rate Ci pi . The union of these is
a new Poisson point process of rate r, and the rest of the algorithm operates
as before. The proofs of all the lemmas and Theorem 2 then proceeds in
exactly the same way as given earlier.

6

Acknowledgments

This research supported by NSF DMS-1418495.

References
Asmussen S, Glynn PW, and Thorisson H (1992) Stationarity detection in the
initial transient problem. ACM Trans Modeling and Computer Simulation
2(2):130–157.
Bickel PJ and Doksum KA (1977) Mathematical Statistics. Prentice Hall,
New Jersey.
Huber M (to appear). Nearly optimal Bernoulli factories for linear functions.
Combin Probab Comput, arXiv:1308.1562.
Keane MS and O’Brien GL (1994) A Bernoulli factory. ACM Trans Modeling
and Computer Simulation 4:213–219.
Latuszyński K, Kosmidis I, Papspiliopoulos O, and Roberts G (2011) Simulating events of unknown probabilities via reverse time martingales. Random Structures Algorithms 38(4):441–452.
Lee A, Doucet A, and Latuszyński K (2014) Perfect simulation using atomic regeneration with application to Sequential Monte Carlo.
arXiv:1407.5770v1.
Nacu S and Peres Y (2005) Fast simulation of new coins from old. Ann Appl
Probab 15(1A):93–115.
20

Walker AJ (1974) New fast method for generating discrete random numbers
with arbitrary frequency distributions. Electronics Letters 10(8):127.

21

