## Background

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)$.)

This question is a continuation of [another question of mine](https://mathoverflow.net/questions/379858/what-are-ways-to-compute-polynomials-that-converge-from-above-and-below-to-a-con), which seeks ways to compute polynomials that converge from above and below to $f$ in a manner that solves the Bernoulli factory problem for $f$. These polynomials form an _approximation scheme_ for $f$. See that question for a formal statement of such approximation schemes.

## New coins from old, smoothly

This time, though, we focus on a **specific approximation scheme**, the one **presented by [Holtz et al. 2011](https://link.springer.com/content/pdf/10.1007/s00365-010-9108-5.pdf), in the paper "New coins from old, smoothly"**.

The scheme involves building polynomials that are shifted upward and downward to approximate $f$ from above and below.

The scheme achieves a convergence rate that generally depends on the smoothness of $f$; in fact, it can achieve the highest convergence rate possible for functions with that smoothness. Roughly speaking, the rate is $O((1/n)^{2(r+\alpha)})$ when $\lambda$ is close to 0 or 1, and $O((1/n)^{r+\alpha})$ elsewhere.

The scheme is as follows:

Let $f$ be a continuous and $r$-times differentiable function—

- that maps [0, 1] to the open interval (0, 1), and
- whose $r$th derivative is $\alpha$-Hölder continuous.

Let $b = 2^s$ and let $s\ge0$ be an integer. Let $Q_{n, r}f$ be a degree $n+r$ approximating polynomial called a _Lorentz operator_ (see the paper for details on the Lorentz operator). Let $n_0$ be the smallest $n$ such that $Q_{n_0, r}f$ has coefficients within [0, 1]. Define the following for every integer $n \ge n_0$ divisible by $n_{0}b$:

- $f_{n_0} = Q_{n_0, r}f$.
- $f_{n} = f_{n/b} + Q_{n, r}(f-f_{n/b})$ for each integer $n > n_0$.

- $\phi(n, \alpha, \lambda) = \frac{\theta_{\alpha}}{n^{\alpha}}+(\frac{\lambda(1-\lambda)}{n})^{\alpha/2}$.

Then the coefficients for the degree $n+r$ polynomial that approximates $f$ are—

- $g(n, r, k) = f_{n}(k/(n+r)) - D * B_{n, \phi(n, r+\alpha, \lambda)}(k/(n+r))$, and
- $h(n, r, k) = f_{n}(k/(n+r)) + D * B_{n, \phi(n, r+\alpha, \lambda)}(k/(n+r))$,

where $B_{n, F}$ is the degree-$n$ Bernstein polynomial of $F$, which I know how to compute.

The problem here is:

- This method uses three constants, namely $s$, $\theta_{\alpha}$, and $D$, that are vaguely defined in the paper. For example, the paper says only that $D$ should be chosen "large enough".
- The paper has no examples of how the scheme works for a selection of functions $f$, making it hard to understand and apply the Holtz method (especially in an automated way).

## Questions

1. What are practical lower bounds for $s$, $\theta_{\alpha}$, and $D$ for the method above, given a function $f$?
2. What are other practical ways to apply the Holtz method to certain functions to find polynomials that converge to those functions?
3. To help me and others (especially programmers) understand and apply the Holtz method, can you provide examples of the method (in terms of Bernstein polynomial approximation schemes) for various functions or classes of functions (such as $C^2$ continuous, twice differentiable, or real analytic functions)?

## My Attempt

The following Python code shows my attempt to implement the Holtz approximation scheme using the SymPy computer algebra library. It includes an implementation of the Lorentz operator as well as a method to find Bernstein polynomials for a function.

```
TAURJ_CACHE = {}
TNJ_CACHE = {}
TAURJ_SYMBOL = symbols("T")

def taurj_x(r, j, n):
    if j == 0:
        return S(1)
    if j == 1:
        return S(0)
    if (r, j, n) in TAURJ_CACHE:
        return TAURJ_CACHE[r, j, n]
    ret = 0
    for l in range(2, j + 1):
        tnj = None
        if (n, l) in TNJ_CACHE:
            tnj = TNJ_CACHE[(n, l)]
        else:
            tnj = sum(
                (k - n * TAURJ_SYMBOL) ** l
                * binomial(n, k)
                * TAURJ_SYMBOL ** k
                * (1 - TAURJ_SYMBOL) ** (n - k)
                for k in range(0, n + 1)
            )
        ret -= tnj * taurj_x(r - l, j - l, n) / factorial(l)
    TAURJ_CACHE[(r, j, n)] = ret.simplify()
    return ret


def taurj(r, j, x, n):
    if j == 0:
        return 1
    if j == 1:
        return 0
    return taurj_x(r, j, n).subs(TAURJ_SYMBOL, x)
    
def lorentz_n_k(func, x, n, r, k, pt=None):
    # kth Bernstein coefficient of the Lorentz operator
    # of degree n and smoothness r at the point
    # pt (or at x if not given).
    if pt == None:
        pt = x
    ret = 0
    for j in range(0, r + 1):
        ret += diff(func, (x, j)).subs(x, S(k) / n) * taurj(r, j, pt, n) / n ** j
    return ret  


def lorentz_n(func, x, n, r, pt=None):
    # Lorentz operator.
    # Create a polynomial that approximates func, which in turn uses
    # the symbol x.  The polynomial's degree is n, is assumed to
    # be at least r times differentiable, and is evaluated
    # at the point pt (or at x if not given).
    # NOTE: If r is 0, this is simply the Bernstein operator.
    if pt == None:
        pt = x
    ret = 0
    diffs = [diff(func, (x, j)) for j in range(0, r + 1)]
    for k in range(0, n + 1):
        rm = 0
        for j in range(0, r + 1):
            rm += diffs[j].subs(x, S(k) / n) * taurj(r, j, pt, n) / n ** j
        ret += rm * bernstein_n_k(n, k, pt)
    return ret



def ncfo_lower(func, x, n, pt, alpha, theta, d):
    # NOTE: Assumes func is bounded away from 0 and 1,
    # but it's trivial to bound any function this way, such as
    # by doing a linear transformation such as
    # $f(x) * \frac{8}{10} + \frac{1}{10}$, then undoing
    # the transformation once the polynomials are found
    r = floor(alpha)
    return lorentz_n(func, x, n, r, pt) - d * bernstein_n(
        phi_n(n, theta, alpha, x), x, n
    )


def ncfo_upper(func, x, n, pt, alpha, theta, d):
    r = floor(alpha)
    return lorentz_n(func, x, n, r, pt) + d * bernstein_n(
        phi_n(n, theta, alpha, x), x, n
    )
      
```

REFERENCES:

- 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).