Magical Fibonacci Formulae

The following Python function computes the Fibonacci sequence, without loops, recursion or floating point arithmetic:

f=lambda n:(b:=2<<n)**n*b//(b*b-b-1)%b

It really does:

>>> [f(n) for n in range(10)]
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34]

How does it work? As a teaser, look at the decimal expansions of $100 / 9899$ and $1000 / 998999$ and see if you notice anything:

$$\frac{100}{9899} = 0.0101020305081321345590463…$$ $$\frac{1000}{998999} = 0.001001002003005008013021034055089144233377610988…$$

Generating functions

We define the Fibonacci sequence as \begin{align*} f(0) &= 0,\\ f(1) &= 1,\\ f(n) &= f(n - 1) + f(n - 2). \end{align*} However, we will also define a rather strange object known as a generating function. It is an ‘infinite polynomial’ (also known as a power series) in variable $x$ whose coefficients correspond to our series of interest. This gives us $$F(x) = 0 + 1x + 1x^2 + 2x^3 + 3x^4 + 5x^5 + 8x^6 + \cdots,$$ $$F(x) = f(0)x^0 + f(1)x^1 + f(2)x^2 + \cdots = \sum_{n=0}^\infty f(n)x^n.$$

Does this function converge? Is this really allowed? All interesting questions we are going to ignore here. Instead, let’s just start doing interesting things with our new object, like taking out the first two terms from the sum:

$$F(x) = f(0)x^0 + f(1)x^1 + \sum_{n=2}^\infty f(n)x^n = x + \sum_{n=2}^\infty f(n)x^n.$$

We can also substitute $f(n) = f(n - 1) + f(n - 2)$ now, because our iteration variable starts at $2$:

\begin{align*} F(x) &= x + \sum_{n=2}^\infty \Big(f(n-1) + f(n-2) \Big)x^n\\ &= x + \sum_{n=2}^\infty f(n-1)x^n + \sum_{n=2}^\infty f(n-2)x^n. \end{align*} Now we can substitute the loop variables, re-insert the $f(0)$ term (which is just zero) into the first sum, and factor out the extra $x$ terms: \begin{align*} F(x) &= x + \sum_{n=1}^\infty f(n)x^{n+1} + \sum_{n=0}^\infty f(n)x^{n+2}\\ &= x - f(0)x^{1} + \sum_{n=0}^\infty f(n)x^{n+1} + \sum_{n=0}^\infty f(n)x^{n+2}\\ &= x + x\sum_{n=0}^\infty f(n)x^{n} + x^2\sum_{n=0}^\infty f(n)x^{n}. \end{align*}

And now using the crucial observation that $F(x) = \sum_{n=0}^\infty f(n)x^{n}$: \begin{align*} F(x) &= x + xF(x) + x^2F(x),\\ (1 - x - x^2)F(x) &= x,\\ F(x) &= \frac{x}{1 - x - x^2}. \end{align*} Wow. Somehow the simple expression $x / (1 - x - x^2)$ ‘contains’ the entire Fibonacci sequence. If you substitute $x = 10^{-3}$ in $F(x)$ you will retrieve our earlier value $$\frac{1000}{998999} = 0.001001002003005008013021034055089144233377610988…$$

An interlude by Binet

Solving $1 - x - x^2 = 0$ with the quadratic formula gives us roots $-(\sqrt{5} \pm 1)/2$, which one might recognize as the (negative) golden ratio $\phi$ and its inverse $\phi^{-1}$: $$\phi = \frac{\sqrt{5} + 1}{2}, \quad \phi^{-1} = \frac{2}{\sqrt{5} + 1} = \frac{2\left(\sqrt{5} - 1\right)}{\left(\sqrt{5} + 1\right)\left(\sqrt{5} - 1\right)} = \frac{\sqrt{5} - 1}{2}.$$

This allows us to factor and use partial fraction decomposition on $F(x)$: $$\frac{x}{1 - x - x^2} = \frac{x}{(1 - \phi x)(1 + \phi^{-1} x)} = \frac{1}{\phi + \phi^{-1}}\left(\frac{1}{1 - \phi x} - \frac{1}{1 + \phi^{-1} x}\right).$$ This is a rather arduous (but strictly elementary) algebraic process so it is much easier to verify by expanding that the identity holds than following along. To verify use the fact that $\phi \cdot \phi^{-1} = \phi - \phi^{-1} = 1$.

If we recall the formula for the geometric series,

$$\frac{1}{1 - r} = \sum_{n=0}^\infty r^n,$$

and apply it to our above expression (once again completely ignoring convergence) while noticing that $\phi + \phi^{-1} = \sqrt{5}$ we find

$$F(x) = \frac{1}{\sqrt{5}} \left( \sum_{n=0}^\infty \phi^n x^n - \sum_{n=0}^\infty {(-\phi})^{-n} x^n \right),$$ $$F(x) = \sum_{n=0}^\infty \frac{1}{\sqrt{5}} \left( \phi^n - {(-\phi})^{-n} \right) x^n.$$

And now for the true magic, recall that we defined $F(x) = \sum_{n=0}^\infty f(n)x^n$, and thus we conclude $$f(n) = \frac{1}{\sqrt{5}}\left(\phi^n - {(-\phi})^{-n} \right).$$

We have recovered Binet’s formula for the Fibonacci numbers, a closed form. Unfortunately evaluating it in Python would eventually fail due to the use of floating-point numbers, which is why this is only an interlude. But it is nevertheless cool:

>>> phi = (1 + 5**0.5) / 2
>>> [(phi**n - (-phi)**-n) / 5**0.5 for n in range(10)]
[0.0, 1.0, 1.0, 2.0, 3.0000000000000004, 5.000000000000001, 8.000000000000002, 13.000000000000002, 21.000000000000004, 34.00000000000001]

Evaluating generating functions

Take another look at $F(10^{-3})$:

$$\frac{1000}{998999} = 0.001001002003005008013021034055089144233377610988…$$

Each next integer in the series starts three places shifted back from the previous one. This makes sense, because $F(10^{-3})$ is the sum of $f(n)10^{-3n}$ for all $n$.

This also means that eventually the method fails, when numbers outgrow three digits and overflow into the previous one. If we ignore that for now, we can study $10^{3n} F(10^{-3})$:

n = 0                   0.001001002003005008013021034055...
n = 1                   1.001002003005008013021034055089...
n = 2                1001.002003005008013021034055089144...
n = 3             1001002.003005008013021034055089144233...
n = 4          1001002003.005008013021034055089144233377...
n = 5       1001002003005.008013021034055089144233377610...
n = 6    1001002003005008.013021034055089144233377610988...
                      ^^^

If we could extract just the highlighted column we have our Fibonacci numbers. Luckily, we can. By flooring we can remove the entire fractional part, and with modulo $10^3$ we can ignore everything after the third digit.

We still have the overflowing issue however. But this is easily fixed by choosing a number $k$ instead of $3$, large enough such that the next Fibonacci number doesn’t overflow into our number of interest. For example the choice $k = n$ works, as the $n$th Fibonacci number certainly won’t overflow $n$ decimal digits, giving formula

$$f(n) = \lfloor 10^{n^2} \cdot F(10^{-n}) \rfloor \bmod 10^{n}.$$

We can generalize much more however. Our choice of $10^3$ and $10^n$ as bases was rather arbitrary. You can use any base $b$, as long as $b$ is large enough. This gives:

$$f(n) = \left\lfloor b^n \cdot F(b^{-1}) \right\rfloor \bmod b,$$ $$f(n) = \left\lfloor b^n \cdot \frac{b^{-1}}{1 - b^{-1} - b^{-2}} \right\rfloor \bmod b,$$ $$f(n) = \left\lfloor \frac{b^{n+1}}{b^2 - b - 1} \right\rfloor \bmod b.$$

This is actually a closed form we can evaluate without the use of floating point arithmetic, as it simply consists of the division of two integers. I have experimentally found that choosing $b = 3^n$ suffices to not overflow for computing $f(n)$, giving our magical formula:

>>> f = lambda n: 3**(n*(n+1)) // (3**(2*n) - 3**n - 1) % 3**n
>>> [f(n) for n in range(10)]
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34]

We can golf this quite a bit by using Python’s “walrus operator” (also called assignment expression by boring people) introduced in Python 3.8. If you write (foo := bar) in an expression the parentheses take on value bar as well as storing bar in a new variable foo available in the rest of the expression. Finally as a bit of flair and efficiency, ${b = 2^{n+1}}$ also works, which can be computed as 2 << n:

>>> f=lambda n:(b:=2<<n)**n*b//(b*b-b-1)%b
>>> [f(n) for n in range(10)]
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34]