The reciprocal Fibonacci constant

November 15, 2023

The reciprocal Fibonacci constant $$C = \sum_{n=1}^{\infty} \frac{1}{F_n} = 3.359885...$$ isn't really useful for anything; its only claim to fame is to be proven irrational (André-Jeannin, 1989).

Computationally, this constant turns out to be somewhat interesting though. The convergence rate of the defining series is decent enough ($\approx 0.2$ digits per term) but unlike series like $e = \sum_{n=0}^{\infty} 1 / n!$ or $\log(2) = \sum_{n=1}^{\infty} (\tfrac{1}{2})^n / n$ which can be computed to $N$ digits in $N^{1+o(1)}$ (quasilinear) time using binary splitting, the reciprocal Fibonacci series ostensibly requires $N^{2+o(1)}$ time; indeed, the LCM of $F_1, \ldots, F_N$ has $O(N^2)$ digits, so we cannot hope to extract a small common denominator for the partial sums.

Theta-like sums

We can achieve $N^{1.5+o(1)}$ time if we express $C$ using theta-like series with $q^{n^2}$-type terms. For example (see the MathWorld entry for some references), we can write $$C = \sqrt{5} \, (L(\phi^{-2}) - L(\phi^{-4})) + \frac{1}{4} \sqrt{5} \, \theta_2^2(\varphi^{-2})$$ where $\phi$ is the golden ratio and $L(q)$ is the Lambert series $$L(q) = \sum_{n=1}^{\infty} \frac{q^n}{1 - q^n} = \sum_{n=1}^{\infty} q^{n^2} \frac{1+q^n}{1-q^n}.$$ I actually know how to compute the Jacobi theta function $\theta_2$ in quasilinear time using the AGM, but I don't know how to do this for $L(q)$. Various Lambert-like series can be expressed in terms of Jacobi theta functions, but I have not seen such a formula for $L(q)$ itself. If anyone has a quasilinear algorithm for $L(q)$, I'd be interested in it.

Quasilinear algorithm

Surprisingly, $C$ can actually be computed in $N^{1+o(1)}$ time using an algorithm that I discovered just yesterday (quite possibly someone thought of this before me but didn't write it down). The idea is to apply binary splitting to Gosper's series $$C = \sum_{n=0}^{\infty} \frac{(-1)^{n(n-1)/2} (F_{4n+3} + (-1)^n F_{2n+2})}{F_{2n+1} F_{2n+2} L_1 L_3 \cdots L_{2n+1}}$$ in which $L_n$ denotes a Lucas number (Gosper p. 66; Arndt).

This is a really weird series; I don't think I've ever seen anything else like it. It has theta-like convergence and a denominator which is a factorial-like product taken over a holonomic sequence. For $N$ digits, you end up multiplying $O(\sqrt N)$ matrices with $O(\sqrt N)$ digits each. Normally in binary splitting based algorithms you have $O(N)$ or $O(N / \log N)$ matrices with $O(\log N)$-bit entries.

Here is a non-optimised reference implementation in Python:

from mpmath import *

def recfib():
    with extraprec(20):
        prec = mp.prec
        N = int(1.200175024907849 * prec**0.5 + 2)
        def bs(a, b):
            if b - a == 1:
                f1, f2 = fib(2*a), fib(2*a+1)
                f1f2 = f1*f2; f22 = f2**2; f4 = f1**2 + 2*(f1f2+f22); q = 2*f1+f2
                return q, (-1)**((a&3)>>1)*(f4+(-1)**a*(f1+f2)), (f1f2+f22)*q
            else:
                m = a + (b - a) // 2
                P1, R1, Q1 = bs(a, m); P2, R2, Q2 = bs(m, b)
                Q2 *= P1; P1 *= P2; R1 *= Q2
                R2 *= Q1; R1 += R2; Q1 *= Q2
                return P1, R1, Q1
        P, R, Q = bs(0, N)
    return R / Q

mp.dps = 1000000
C = recfib()

This takes 0.8 seconds on my laptop (using the GMPY backend for mpmath).

Just because it's an amusing benchmark problem, I've put a more optimised implementation in FLINT. With this I can compute 1 million digits in 0.1 seconds and 1 billion digits in 12 minutes (only using 2 cores for the latter because otherwise I run out of memory):

~/src/flint$ build/examples/pi 1000005 -constant reciprocal_fibonacci -threads 8
precision = 3321949 bits... cpu/wall(s): 0.558 0.111
virt/peak/res/peak(MB): 548.91 923.86 55.47 55.47
[3.35988566624317755317{...999964 digits...}31740420902017255257 +/- 3.28e-1000005]

~/src/flint$ build/examples/pi 1000000005 -constant reciprocal_fibonacci -threads 2
precision = 3321928116 bits... cpu/wall(s): 1009.31 718.785
virt/peak/res/peak(MB): 16534.32 18638.34 15483.07 18098.75
[3.35988566624317755317{...999999964 digits...}02927932123820488401 +/- 2.91e-1000000005]

In practice computing $C$ is about two times slower than $\pi$ and two times faster than $\log(2)$.

I'm now curious to what extent Gosper's series can be generalised. If I'm reading Gosper's notes correctly, this particular formula should be possible to generalise to certain other Fibonacci-like sequences, but maybe there are also entirely different interesting numbers which admit a series representation with similar overall form, i.e. something like $\sum_n A(n) / (B(n) \prod_i C(i))$ with $A$, $B$, $C$ holonomic?

Addendum: I believe evaluating the Lambert series and theta series using binary splitting with arithmetic in $\mathbb{Q}(\phi)$ also leads to a quasilinear algorithm. It would be interesting to compare this to the Gosper series.

fredrikj.net | Blog index | RSS feed | Follow me on Mastodon | Become a sponsor