We're given a coin that shows heads with an unknown probability, $\lambda$. The goal is to use that coin (and possibly also a fair coin) to build a "new" coin that shows heads with a probability that depends on $\lambda$, call it $f(\lambda)$. This is the Bernoulli factory problem, and it can be solved only for certain functions $f$. (For example, flipping the coin twice and taking heads only if exactly one coin shows heads, we can simulate the probability $2\lambda(1-\lambda)$.)

Specifically, the only functions that can be simulated this way are continuous and polynomially bounded, and map $[0, 1]$ to $[0, 1]$, as well as $f=0$ and $f=1$. These functions are called factory functions in this question. (A function $f(x)$ is polynomially bounded if both $f$ and $1-f$ are bounded below by min($x^n$, $(1-x)^n$) for some integer $n$ (Keane and O'Brien 1994). This implies that $f$ admits no roots on (0, 1) and can't take on the value 0 or 1 except possibly at 0 and/or 1.)

Polynomials that approach a factory function

There are two algorithms (one by Thomas and Blanchet 2012, another by Łatuszyński et al.) to simulate a factory function via two sequences of polynomials that converge from above and below to that function. Roughly speaking, the algorithms work as follows:

  1. Generate U, a uniform random number in $[0, 1]$.
  2. Flip the input coin (with a probability of heads of $\lambda$), then build an upper and lower bound for $f(\lambda)$, based on the outcomes of the flips so far. In this case, these bounds come from two degree-$n$ polynomials that approach $f$ as $n$ gets large, where $n$ is the number of coin flips so far in the algorithm.
  3. If U is less than or equal to the lower bound, return 1. If U is greater than the upper bound, return 0. Otherwise, go to step 2.

The result of the algorithm is 1 with probability exactly equal to $f(\lambda)$, or 0 otherwise.

However, both algorithms require the polynomial sequences to meet certain requirements: Roughly speaking, the sequences must be of Bernstein-form polynomials that converge from above and below to a factory function as follows:

  • Each sequence's polynomials must have coefficients lying in [0, 1], and be of increasing degree.
  • The degree-n polynomials' coefficients must lie at or "inside" those of the previous upper polynomial and the previous lower one (once the polynomials are elevated to degree n).

See below for a formal statement of these polynomials.

The answer below and my supplemental notes include formulas for computing these polynomials for large classes of factory functions, but none of them ensure a finite expected number of coin flips in general, and it is suspected that a finite number of flips isn't possible unless the factory function is C2 continuous. (Indeed, see related results that only $C^2$ continuous functions can be approximated by polynomials at a rate of $O(1/n^{2+\epsilon})$ or better, for positive $\epsilon$, by Holtz et al. (2011).)


Thus the questions are:

  1. Given a $C^2$ continuous function, are there practical formulas to compute polynomials that—

    • meet the formal statement below, and
    • can be used to simulate that function with a finite expected running time?
  2. Are there other practical formulas to approximate specific factory functions with polynomials that meet the formal statement below?

One example to question 1 is the function $\min(2\lambda,1-\epsilon)$ on the domain $(0, 1/2-\epsilon)$ given in Nacu and Peres 2005. On the other hand, the method of Holtz et al. 2011 is hard to apply since, among other things, the paper's results are only asymptotic.


  • A related question seeks a practical way to apply the Holtz method.
  • A related question seeks ways to approximate concave functions.
  • This question is part of a broad request for algorithms to solve the Bernoulli factory problem for specific factory functions and irrational constants (see question 1 of the linked section). Especially algorithms that rely solely on rational arithmetic and don't calculate square roots or transcendental functions directly. However, this question is not the place to answer this request directly; see my Cross Validated question instead.

Formal Statement

More formally, for the approximation schemes I am looking for, there exist two sequences of polynomials, namely—

  • $g_{n}(\lambda): =\sum_{k=0}^{n}a(n, k){n \choose k}\lambda^{k}(1-\lambda)^{n-k}$, and
  • $h_{n}(\lambda): =\sum_{k=0}^{n}b(n, k){n \choose k}\lambda^{k}(1-\lambda)^{n-k}$,

for every integer $n\ge1$, such that—

  1. $0\le a(n, k)\le b(n, k)\le1$,
  2. $\lim_{n}g_{n}(\lambda)=\lim_{n}h_{n}(\lambda)=f(\lambda)$ for every $\lambda\in[0,1]$, and
  3. for every $m<n$, both $(g_{n} - g_{m})$ and $(h_{m} - h_{n})$ have non-negative coefficients once $g_{n}$, $g_{m}$, $h_{n}$, and $h_{m}$ are rewritten as degree-$n$ polynomials in Bernstein form,

where $f(\lambda)$ is continuous on $[0, 1]$ (Nacu and Peres 2005; Holtz et al. 2011), and the goal is to find the appropriate values for $a(n, k)$ and $b(n, k)$.

It is allowed for $a(n, k)\lt0$ for a given $n$ and some $k$, in which case all $a(n, k)$ for that $n$ are taken to be 0 instead. It is allowed for $b(n, k)\gt1$ for a given $n$ and some $k$, in which case all $b(n, k)$ for that $n$ are taken to be 1 instead.

Equivalently, for the approximation schemes I am looking for, $f(\lambda)$ can be rewritten as— $$\sum_{n\ge 0} F_n(\lambda),$$ where—

  • $\sum_{n>N} F_n(\lambda) \le \psi_N(\lambda)$ for every $\lambda \in [0,1]$,
  • $(\psi_n)$ is a "nonincreasing sequence of positive functions ... that converges uniformly to zero", such as a sequence of constants, and
  • each $F_n$ is a degree-$n$ polynomial with Bernstein coefficients 0 or greater,

and the goal is to find the appropriate functions $F_n$ and $\psi_n$ for each $n$ (Holtz et al. 2011, Lemma 4).


  • Thomas, A.C., Blanchet, J., "A Practical Implementation of the Bernoulli Factory", arXiv:1106.2508v3 [stat.AP], 2012.
  • Łatuszyński, K., Kosmidis, I., Papaspiliopoulos, O., Roberts, G.O., "Simulating events of unknown probabilities via reverse time martingales", arXiv:0907.4018v2 [stat.CO], 2009/2011.
  • Keane, M. S., and O'Brien, G. L., "A Bernoulli factory", ACM Transactions on Modeling and Computer Simulation 4(2), 1994.
  • Holtz, O., Nazarov, F., Peres, Y., "New Coins from Old, Smoothly", Constructive Approximation 33 (2011).
  • Nacu, Şerban, and Yuval Peres. "Fast simulation of new coins from old", The Annals of Applied Probability 15, no. 1A (2005): 93-115.
  • Farouki, R.T., and Rajan, V.T., "Algorithms for polynomials in Bernstein form", Computer Aided Geometric Design 5(1), 1988.
  • $\begingroup$ Also on the subject of bounded polynomials: mathoverflow.net/questions/79921/… $\endgroup$
    – Matt F.
    Jan 6 at 19:50
  • $\begingroup$ Can we assume that $f$ is non-decreasing, and if so do you require that each approximating polynomial is non-decreasing? That would seem natural for a Bernoulli factory problem with approximations for $\min(\lambda,x)$. $\endgroup$
    – Matt F.
    Jan 6 at 20:04
  • $\begingroup$ @MattF: No, f is not necessarily non-decreasing; it can be any continuous function that is bounded away from 0 and 1 except at the points 0 and 1. Also, each approximating upper polynomial is no greater than the previous and each approximating lower polynomial is no less than the previous. However, if you have an answer that works for any non-decreasing f, you are welcome to post it. $\endgroup$
    – Peter O.
    Jan 7 at 1:18
  • $\begingroup$ @MattF: Also, there is an intriguing suggestion from Thomas and Blanchet (2012): use multiple pairs of polynomial sequences that converge to $f$ in the manner given in question 1 of this post, where each pair is optimized for particular ranges of $\lambda$. Answers to this question that suggest multiple approximation schemes of this kind for the same function are also appreciated. $\endgroup$
    – Peter O.
    Jan 7 at 1:31

While the following does not fully answer my question, I make the following notes.

The following notes relate to finding polynomials for certain functions, but none of them have a finite expected running time in general. Other users are encouraged to add other answers to my question.

Polynomial Schemes

The following inequality from Nacu and Peres 2005 is useful: $$|\mathbb{E}[f(X/n)]-f(k/(2n))| \le \kappa_f(n), \tag{1}$$ where—

  • $\kappa_f(n)$ is a function that depends only on $f$ and $n$ and makes the inequality hold true for all $f$ belonging to a given class of functions (such as Lipschitz continuous or twice-differentiable functions), and
  • $X$ is a hypergeometric($2n$, $k$, $n$) random variable.

Notably, if a function $f$ is such that (1) holds true, then in general, for all $n$ that are powers of 2—

  • The $k$th Bernstein coefficient for the upper polynomial of $n$th degree is $f(k/n) + \delta(n)$, and
  • the $k$th Bernstein coefficient for the lower polynomial of $n$th degree is $f(k/n) - \delta(n)$,

where $n$ is the polynomial's degree and $\delta(n)$ is a solution to the following functional equation: $$\delta(n) = \delta(n*2) + \kappa_f(n), $$or equivalently, the linear recurrence $\delta(n) = \delta(n+1) + \kappa_f(2^n)$.

For example—

  1. If $f$ is Lipschitz continuous with Lipschitz constant $C$
    • $\kappa_f(n) = C/\sqrt{2*n}$ , so
    • $\delta(n) = (1+\sqrt{2})C/\sqrt{n}$ , and
  2. If $f$ is twice differentiable and $M$ is not less than the highest value of $|f′′(x)|$ for any $x$ in [0, 1]—
    • $\kappa_f(n) = M/(4*n)$, so
    • $\delta(n) = M/(2*n)$

(Nacu and Peres 2005, Proposition 10). As experiments show, these bounds are far from tight, but they can generally be improved only by an order of magnitude.

A new result of mine exploits the results of Nacu and Peres to derive a new approximation scheme for Hölder continuous functions. For example, if $f$ is $\alpha$-Hölder continuous with Hölder constant $M$

  • $\kappa_f(n) = M(1/(7n))^{\alpha/2}$, so
  • $\delta(n) = \frac{M(2/7)^{\alpha/2}}{(2^{\alpha/2}-1)n^{\alpha/2}}$.

(For proofs and additional notes on approximation schemes I'm looking for, see my supplemental notes.)

Moreover, due to Jensen's inequality:

  • If $f$ is concave in $[0, 1]$, the Bernstein coefficients for the lower polynomials are simply $f(k/n)$.
  • If $f$ is convex in $[0, 1]$, the Bernstein coefficients for the upper polynomials are simply $f(k/n)$.

If $f$ equals 0 or 1 at the points 0 and/or 1, then $f$ can be transformed to bound it away from 0 and/or 1. For example, if $f(0)=0$, then it can be divided by a function that bounds $f$ from above. This case is too complicated to detail in this answer; see my supplemental notes for details.

Of course, the methods above are not the only way to build polynomials that converge to a function in the manner asked for by my question, and this answer doesn't solve all the issues I mention in my question.

Also, the procedure shifts the whole approximation by a constant, whereas it should be possible to provide different approximations to the same function that are optimized for particular probabilities of heads, as suggested in Thomas and Blanchet 2012.

Also, I have found the following result, which shows that any factory function admits an approximation scheme with polynomials (in a manner that solves the Bernoulli factory problem). This includes continuous functions that are not Hölder continuous. The method of proving the result goes back to Nacu and Peres (2005). This is one of several new results relating to this problem; see the appendix on my page on approximation schemes. However, again, this scheme doesn't have a finite expected running time in general.

Result: Let $f(\lambda)$ be a continuous function that maps [0, 1] to [0, 1], and let $X$ be a hypergeometric($2n$, $k$, $n$) random variable. Let $\omega(x)$ be a modulus of continuity of $f$ (a non-negative and nondecreasing function in the interval [0, 1], for which $\omega(x) = 0$, and for which $|f(x)-f(y)|\le\omega(|x-y|)$ for all $x$ in [0, 1] and all $y$ in [0, 1]). If $\omega$ is concave on [0, 1], then the expression—$$|\mathbb{E}[f(X/n)]-f(k/(2n))|, $$is bounded from above by—

  • $\omega(\sqrt{\frac{1}{8n-4}})$, for all n≥1 that are integer powers of 2,
  • $\omega(\sqrt{\frac{1}{7n}})$, for all n≥4 that are integer powers of 2, and
  • $\omega(\sqrt{\frac{1}{2n}})$, for all n≥1 that are integer powers of 2.

Proof: ω is assumed to be non-negative because absolute values are non-negative. To prove the first and second bounds: abs(E[f(X/n)] − f(k/(2 * n))) ≤ E[abs(f(X/n) − f(k/(2 * n))] ≤ E[ω(abs(X/nk/(2 * n))] ≤ ω(E[abs(X/nk/(2 * n))]) (by Jensen's inequality and because ω is concave) ≤ ω(sqrt(E[abs(X/nk/(2 * n))]2)) = ω(sqrt(Var[X/n])) = ω(sqrt((k*(2 * nk)/(4*(2 * n −1)*n2)))) ≤ ω(sqrt((n2/(4*(2 * n −1)*n2)))) = ω(sqrt((1/(8*n −4)))) = ρ, and for all n ≥4 that are integer powers of 2, ρω(sqrt(1/(7*n))). To prove the third bound: abs(E[f(X/n)] − f(k/(2 * n))) ≤ ω(sqrt(Var[X/n])) ≤ ω(sqrt(1/(2*n))). □

The only technical barrier, though, is to solve the functional equation—$$\delta(n) = \delta(2n) + \kappa(n), $$ where $\kappa(n)$ is one of the bounds proved above, such as $\omega\left(\sqrt{\frac{1}{8n-4}}\right)$.

In general, the solution to this equation is—$$\delta(2^m) = \sum_{k=m}^\infty \kappa(2^k), $$ where $n = 2^m$ and $m\ge0$ are integers, provided the sum converges.

In this case, the third bound has a trivial solution when $\omega(x)$ is of the form $cx^\alpha$, but not in general.

Now, we approximate $f$ with polynomials in Bernstein form of power-of-two degrees. These are two sequences of polynomials that converge to $f$ from above and below, such that their coefficients "increase" and "decrease" just as the polynomials themselves do. In general, for all $n\ge1$ that are integer powers of 2:

  • The $k$th Bernstein coefficient for the upper polynomial of $n$th degree is $f(k/n) + \delta(n)$.
  • The $k$th Bernstein coefficient for the lower polynomial of $n$th degree is $f(k/n) - \delta(n)$.

Thus, for the polynomials of degree $n$, $\delta(n)$ is the amount by which to shift the approximating polynomials upward and downward.


Theorem 26 of Nacu and Peres (2005) and the proof of Keane and O'Brien (1994) give general ways to simulate continuous factory functions $f(\lambda)$ on the interval $[0, 1]$. The former is limited to functions that are bounded away from 0 and 1, while the latter is not. However, both methods don't provide formulas of the form $f(k/n) \pm \delta(n)$ that work for a whole class of factory functions (such as formulas of the kind given earlier in this answer). For this and other reasons, given below, both methods are impractical:

  • Before a given function $f$ can be simulated, the methods require computing the necessary degree of approximation (finding $k_a$ or $s_i$ for each polynomial $a$ or $i$, respectively). This work has to be repeated for each function $f$ to be simulated.
  • Computing the degree of approximation involves, among other things, checking whether the approximating polynomial is "close enough" to $f$, which can require a symbolic maximization.
  • This computation gets more and more time-intensive with increasing degree.
  • For a given $f$, it's not guaranteed whether the $k_a$'s (or $s_i$'s) will show a pattern or keep that pattern "forever" — especially since only a finite number of approximation degrees can be computed with these methods.



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.