Numerically stable Cohen-Villegas-Zagier algorithm for alternating series

This post is a continuation of the work previously presented in Faster numerical inverse Laplace transform in mpmath. The Cohen-Villegas-Zagier algorithm has gained popularity for accelerating the convergence of alternating series, due to its the ease of implementation and effectiveness. However, the forward recursion algorithm proposed in (Cohen et al, 2000) implemented in computer algebra systems such as PARI or mpmath, suffers from cancellation issues when the number of terms $n$ exceeds the threshold $1.31 D$, for $D$ for decimal digits.

In this post, we derive a numerically stable backward algorithm suitable for fixed and arbitrary precision computation, accompanied with the derivation of rigorous error bounds for the evaluation of the last term.

Cohen-Villegas-Zagier algorithm

[Proposition 1 in (Cohen et al, 2000)] For integers $n\ge k\ge 0$ set $$ d_n = \frac{(3+\sqrt{8})^n + (3-\sqrt{8})^n}{2} $$ and $$ c_{n,k} = (-1)^k \left(d_n - \sum_{m=0}^k \frac{n}{n+m} \binom{n+m}{2m} 2^{2m}\right) = (-1)^k \sum_{m=k+1}^n \frac{n}{n+m} \binom{n+m}{2m} 2^{2m} $$

Forward recurrence

Initialize: $$ d = (3 + \sqrt{8})^n; \quad d = (d + 1/d)/2; \quad b = -1; \quad c = -d; \quad s = 0; $$ For $k=0,\ldots, n-1$, repeat: $$ c = b - c, \quad b = \frac{2 (k + n) (k-n)}{(2k + 1) (k+1)}b, \quad s = s + c \cdot a_k; $$

Output: $s / d$.

To compute the coefficient $w_k = c_k / d_n$, so that $S = \sum_{k=0}^{\infty}(-1)^k a_k \approx \sum_{k=0}^n w_k a_k$, change the initialization to: $b = -1/d$ and $c = -1$. For large $n$ and fixed precision, the evaluation of $(d + 1/d)/2$ can overflow, limiting the maximum number of terms for a prescribed precision.

The Python mpmath code to generate the coefficients $w_{n,k}$ using the forward recurrence is

def forward_algorithm(n, digits=100):
    prec = mp.dps
    mp.dps = digits

    d = (mp.mpf('3') + mp.sqrt(mp.mpf('8'))) ** n
    d = (d + mp.one / d) / mp.mpf('2')
    g = -mp.one
    c = -d

    wk = []
    for k in range(n):
        c = g - c
        g = 2 * (k + n) * (k - n) * g / ((2 * k + 1) * (k + mp.one))

        wk.append(c / d)

    mp.dps = prec

    return wk

A numerically stable recursive algorithm

First, observe that the last term $k=n-1$ is simply

$$ c_{n,n-1} = \frac{n}{2n} 2^{2n} = 2^{2n-1}. $$

Also, consider the trigonometric form of the term $d_n$

$$ d_n = \cosh \left(2n \text{ arcsinh}(1)\right) = \cosh\left(2n \log(1+ \sqrt{2})\right). $$

Using the previously defined normalized coefficients $w_{n,k} = c_{n,k}/d_n$, the last normalized coefficient $w_{n, n-1}$ is given by

$$ w_{n, n-1} = \frac{2^{2n-1}}{d_n} = 2^{2n-1} \text{sech}\left(2n \log(1 +\sqrt{2})\right). $$

As mentioned, for large values $n$, a direct computation of the last term $w_{n, n-1}$ in double-precision floating-point arithmetic poses numerical issues (overflow on $2^{2n-1}$ and underflow when computing $\text{sech}$). To circumvent this problem, consider the alternative representation of $w_{n, n-1}$ in terms of the hypergeometric function $_1F_0$

$$ 2^{2n-1}\text{sech}\left(2n \log(1 +\sqrt{2})\right) = 2^{2n} e^{-2n \log(1 +\sqrt{2})} \, _1F_0\left(1;;-e^{-2n \log(1 +\sqrt{2})}\right). $$ Then $$ w_{n, n-1} = \left(\frac{4}{(1 + \sqrt{2})^2}\right)^n \sum_{m=0}^{\infty} (-1)^m e^{-2n m\log(1 + \sqrt{2})} = \sum_{m=0}^{\infty} (-1)^m \left(\frac{4}{(1 + \sqrt{2})^{4m+2}} \right)^n. $$

The expansion is absolutely convergent, with rate of convergence

$$ (1 + \sqrt{2})^{-4n} = (17 + 12\sqrt{2})^{-n} \approx 33.97^{-n}. $$

The number of terms required for a given precision $\epsilon = 2^{-r}$ ($r = 52$ for double-precision) can be obtained solving the following equation involving the absolute value of the truncated term

$$ \left(\frac{4}{(1 + \sqrt{2})^{4m+2}} \right)^n = \epsilon, \quad m = \frac{1}{4}\left(\frac{\log(4 \epsilon^{-1/n})}{\log(1 + \sqrt{2})} - 2\right), \quad M = \lceil m \rceil. $$

The first few terms are sufficient to accurately compute $w_{n, n-1}$ for large $n$ when targeting small precision

$$ \left(\frac{4}{3 + \sqrt{8}}\right)^n - \left(\frac{4}{99 + 70\sqrt{2}}\right)^n + \left(\frac{4}{3363 + 2378\sqrt{2}}\right)^n \approx 0.686292^n - 0.020202^n + 0.000595^n. $$

In particular, for double-precision the leading term is sufficiently accurate for $n\ge 10$. Note, however, that for small $n$, a direct computation of $w_{n, n-1}$ is not problematic. To conclude, we estimate the maximum number of terms $n$ so that the last term $w_{n, n-1}$ is representable (no underflow occurs, no subnormal numbers). Consider the smallest representable number $2^{-p}$, and solve for $n$ equating with the leading term

$$ \left(\frac{4}{3 + \sqrt{8}}\right)^n = 2^{-p} \Longleftrightarrow n = \Bigg\lfloor -p \frac{\log(2)}{\log(12 - 8 \sqrt{2})} \Bigg\rfloor. $$

For example, in double-precision $p=1023$, then the maximum $n = 1883$.

Backward recurrence

Using the previous results, and inverting the recurrence for the term $b$, we obtain the backward recurrence described below

Initialize:

$$ M = \Bigg\lceil \frac{1}{4}\left(\frac{\log(4 \epsilon^{-1/n})}{\log(1 + \sqrt{2})} - 2\right)\Bigg\rceil $$ $$ b = (-1)^{n-1}\sum_{m=0}^{M} (-1)^m \left(\frac{4}{(1 + \sqrt{2})^{4m+2}} \right)^n, \quad w = b, \quad s = w \cdot a_n. $$

For $k = 0, \ldots, n-2$, repeat:

$$ b = \frac{(2(n-k)-1)(k-n)}{2(k+1)(2n-k-1)} b, \quad w = b - w, \quad s = s + w \cdot a_{n-1-k} $$

The Python mpmath code to generate the coefficients $w_{n,k}$ using the backward recurrence is

def backward_algorithm(n, digits=100):
    prec = mp.dps
    mp.dps = digits

    r = mp.one + mp.sqrt(mp.mpf('2'))
    m = (mp.log(4 * mp.eps ** (-1 / n)) / mp.log(r) - 2) / 4
    M = int(mp.ceil(m))

    bnm1 = mp.zero
    for k in range(M + 1):
        bnm1 += (-mp.one) ** k * (mp.mpf('4') / r ** (4*k + 2)) ** n

    bnm1 *= (-1) ** (n - 1)
    cnm1 = bnm1

    wk = [cnm1]
    for k in range(n-1):
        bnm1 = (2 * (n - k) - 1) * (k - n) * bnm1 / (2 * (k + mp.one) * (2 * n - k - mp.one))
        cnm1 = bnm1 - cnm1
        wk.append(cnm1)

    mp.dps = prec

    return wk[::-1]

Numerical experiments

To perform numerical experiments about the loss of precision of the forward and backward recurrence, we use as a reference a closed-form expression for the coefficients $c_{n,k}$ in terms of the hypergeometric function $_3F_2$, obtained using Mathematica

$$ c_{n,k} = (-1)^k \frac{4^{k+1} n}{(k + n + 1)} \binom{n + k + 1}{2(k+1)}\, _3F_2(1, k - n + 1, k + n + 1; k + 3/2, k + 2; -1). $$ The corresponding mpmath code is

def exact_algorithm(n, digits=100):
    prec = mp.dps
    mp.dps = digits

    d = (mp.mpf('3') + mp.sqrt(mp.mpf('8'))) ** n
    d = (d + mp.one / d) / mp.mpf('2')

    wk = []
    for k in range(n):
        c = (-mp.one) ** k * 4 ** (k + 1) * n * mp.binomial(k + n + 1, 2 * (k + 1)) / (k + n + 1)
        c *= mp.hyp3f2(1, k - n + 1, k + n + 1, k + mp.mpf('3/2'), k + 2, -1)
        wk.append(c / d)
    
    mp.dps = prec

    return wk

In addition, we include in the benchmark a pure Python implementation targeting double-precision

def backward_algorithm_float(n):
    if n >= 10:
        M = 0
    else:
        m = 10.223716117646134 / n - 0.10678014932130256
        M = int(m) + 1

    v = 0.6862915010152396
    bnm1 = v ** n
    for m in range(1, M + 1):
        v *= 0.029437251522859414
        bnm1 += (-1) ** m * v ** n
        
    bnm1 *= (-1) ** (n - 1)
    cnm1 = bnm1
    
    wk = [cnm1]
    for k in range(n-1):
        bnm1 = (2 * (n - k) - 1) * (k - n) * bnm1 / (2 * (k + 1) * (2 * n - k - 1))
        cnm1 = bnm1 - cnm1
        wk.append(cnm1)
        
    return wk[::-1]

The following plot compares the three recurrences (forward, backward and backward - float) evaluated with 15 digits of precision, against the closed-form expression computed with 300 digits of precision. Observe that the maximum relative forward recurrence grows quickly as the number of terms increases, whereas the backward recurrence (mpmath and pure Python) consistently follows the double-precision limit, $2^{-52} \approx 2.2{e-1}6$, only showing a slight increase for large $n$.

Figure 1: Maximum relative errors of the three recurrences compared (forward, backward and backward - float) against the reference value computed with 300 digits of precision.

Alternating sum implementation

To conclude, given a list of values $a$, we suggest a possible implementation of a sumalt routine using mpmath.

def sumalt(a, method='forward', digits=100):
    n = len(a)
    if method == 'forward':
        c = forward_algorithm(n, digits)
    elif method == 'exact':
        c = exact_algorithm(n, digits)
    elif method == 'backward':
        c = backward_algorithm(n, digits)

    s = mp.zero
    if method in ('forward', 'exact'):
        for k in range(n):
            s += c[k] * a[k]
    else:
        for k in range(n-1, -1, -1):
            s += c[k] * a[k]

    return s

References

H. Cohen, F. Rodriguez Villegas, and D. Zagier (2000). Convergence acceleration of alternating series. Experimental Mathematics, 9(1):3-12