This code, somewhat surprisingly, generates Fibonacci numbers.
def fib(n):
return (4 << n*(3+n)) // ((4 << 2*n) - (2 << n) - 1) & ((2 << n) - 1)
In this blog post, I’ll explain where it comes from and how it works.
Before getting to explaining, I’ll give a whirlwind background overview of Fibonacci numbers and how to compute them. If you’re already a maths whiz, you can skip most of the introduction, quickly skim the section “Generating functions” and then read “An integer formula”.
Overview
The Fibonacci numbers are a well-known sequence of numbers:
$$1, 1, 2, 3, 5, 8, 13, 21, 34, 55, \ldots$$
The $n$th number in the sequence is defined to be the sum of the previous two, or formally
by this recurrence relation:
$$\begin{eqnarray} \mathrm{Fib}(0) &=& 1 \\ \mathrm{Fib}(1) &=& 1 \\ \mathrm{Fib}(n) &=& \mathrm{Fib}(n - 1) + \mathrm{Fib}(n - 2) \end{eqnarray}$$
I’ve chosen to start the sequence at index 0 rather than the more usual 1.
Computing Fibonacci numbers
There’s a few different reasonably well-known ways of computing the sequence. The obvious recursive implementation is slow:
def fib_recursive(n):
if n < 2: return 1
return fib_recursive(n - 1) + fib_recursive(n - 2)
An iterative implementation works in $O(n)$ operations:
def fib_iter(n):
a, b = 1, 1
for _ in xrange(n):
a, b = a + b, a
return b
And a slightly less well-known matrix power implementation works in $O(\mathrm{log}\ n)$ operations.
def fib_matpow(n):
m = numpy.matrix('1 1 ; 1 0') ** n
return m.item(0)
The last method works by considering the a and b in fib_iter
as sequences, and noting that
$$ \left(\begin{matrix} a_{n+1} \\ b_{n+1} \\ \end{matrix}\right) = \left(\begin{matrix} 1 & 1 \\ 1 & 0 \\ \end{matrix}\right) \left(\begin{matrix} a_n \\ b_n \\ \end{matrix}\right) $$
From which follows
$$ \left( \begin{array}{c} a_{n} \\ b_{n} \end{array} \right) = \left( \begin{array}{cc} 1 & 1 \\ 1 & 0 \end{array} \right) ^ n \left( \begin{array}{c} 1 \\ 1 \end{array} \right) $$
and so if
$m = \left(\begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix} \right)^n$ then $b_n = m_{11}$ (noting that unlike Python, matrix indexes are usually 1-based).
It’s $O(\mathrm{log}\ n)$ based on the assumption that numpy’s matrix power does something like exponentation by squaring.
Another method is to find a closed form for the solution of the recurrence relation. This leads to the real-valued formula: $\mathrm{Fib}(n) = (\phi^{n+1} - \psi^{n+1}) / \sqrt{5}$ where $\phi = (1 + \sqrt{5}) / 2$ and $\psi = (1 - \sqrt{5}) / 2$. The practical flaw in this method is that it requires arbitrary precision real-valued arithmetic, but it works for small $n$.
def fib_phi(n):
phi = (1 + math.sqrt(5)) / 2.0
psi = (1 - math.sqrt(5)) / 2.0
return int((phi ** (n+1) - psi ** (n+1)) / math.sqrt(5))
Generating Functions
A generating function for an arbitrary sequence $a_n$ is the infinite sum $\Sigma_n a_nx^n$. In the specific case of the Fibonacci numbers, that means $\Sigma_n \mathrm{Fib}(n)x^n$. In words, it’s an infinite power series, with the coefficient of $x^n$ being the $n$th Fibonacci number.
Now, $\mathrm{Fib}(n+2) = \mathrm{Fib}(n+1) + \mathrm{Fib}(n)$
Multiplying by $x^{n+2}$ and summing over all $n$, we get:
$$\Sigma_n\mathrm{Fib}(n+2)x^{n+2} = \Sigma_n\mathrm{Fib}(n+1)x^{n+2} + \Sigma_n\mathrm{Fib}(n)x^{n+2}$$
If we let $F(x)$ to be the generating function of $\mathrm{Fib}$, which is defined to be $\Sigma_n\mathrm{Fib}(n)x^n$ then this
equation can be simplified:
$$F(x) - x - 1 = x(F(x) - 1) + x^2F(x)$$
and simplifying,
$$F(x) = xF(x) + x^2F(x) + 1$$
We can solve this equation for $F$ to get
$$F(x) = \frac{1}{1 - x - x^2}$$
It’s surprising that we’ve managed to find a small and simple formula which captures all of the Fibonacci numbers, but it’s not yet obvious how we can use it. We’ll get to that in the next section.
A technical aside is that we’re going to want to evaluate $F$ at some values of $x$, and we’d like the power series to converge. We know the Fibonacci numbers grow like $\phi^n$ and that geometric series $\Sigma_n a^n$ converge if $|a|<1$, so we know that if $|x| < 1/\phi \simeq 0.618$ then the power series converges.
An integer formula
Now we’re ready to start understanding the Python code.
To get the intuition behind the formula, we’ll evaluate the generating function $F$ at $10^{-3}$.
$$ \begin{align*} & F(x) = \frac{1}{1 - x - x^2} \\ & F(10^{-3}) = \frac{1}{1 - 10^{-3} - 10^{-6}} \\ & = 1.001\,002\,003\,005\,008\,013\,021\,034\,055\,089\,144\,233\,377\,610\,988\,599\,588\,187\,\ldots \end{align*}$$
Interestingly, we can see some Fibonacci numbers in this decimal expansion: $1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89$. That seems magical and surprising, but it’s because $F(10^{-3}) = \mathrm{Fib}(0) + \mathrm{Fib}(1)/10^3 + \mathrm{Fib}(2)/10^6 + \mathrm{Fib}(3)/10^9 + \ldots$.
In this example, the Fibonacci numbers are spaced out at multiples of $1/1000$, which means once they start getting bigger that 1000 they’ll start interfering with their neighbours. We can see that starting at 988 in the computation of $F(10^{-3})$ above: the correct Fibonacci number is 987, but there’s a 1 overflowed from the next number in the sequence causing an off-by-one error. This breaks the pattern from then on.
But, for any value of $n$, we can arrange for the negative power of 10 to be large enough that overflows don’t disturb the $n$th Fibonacci. For now, we’ll just assume that there’s some $k$ which makes $10^{-k}$ sufficient, and we’ll come back to picking it later.
Also, since we’d like to use integer maths (because it’s easier to code), let’s multiply by $10^{kn}$, which also puts the $n$th Fibonacci
number just to the left of the decimal point, and simplify the equation.
$$ 10^{kn} F(10^{-k}) = \frac{10^{kn}}{1 - 10^{-k} - 10^{-2k}} = \frac{10^{kn+2k}}{10^{2k} - 10^{k} - 1} $$
If we take this result modulo $10^k$, we’ll get the $n$th Fibonacci number (again, assuming we’ve picked $k$ large enough).
Before proceeding, let’s switch to base 2 rather than base 10, which changes nothing but will make it easier to program.
$$ 2^{kn} F(2^{-k}) = \frac{2^{k(n+2)}}{2^{2k} - 2^{k} - 1} $$
Now all that’s left is to pick a value of $k$ large enough so that $\mathrm{Fib}(n+1) < 2^k$. We know that the Fibonacci numbers grow like $\phi^n$, and $\phi < 2$, so $k = n+1$ is safe.
So! Putting that together:
$$ \begin{align*} \mathrm{Fib}(n) & \equiv 2^{(n+1)n}F(2^{-(n+1)})\ \mathrm{mod}\ 2^{n+1}\\ & \equiv \frac{2^{(n+1)(n+2)}}{(2^{n+1})^2 - 2^{n+1} - 1}\ \mathrm{mod}\ 2^{n+1} \\ & \equiv \frac{2^{(n+1)(n+2)}}{2^{2n+2} - 2^{n+1} - 1}\ \mathrm{mod}\ 2^{n+1} \end{align*} $$
If we use left-shift notation that’s available in python, where $a\ mathrm{«}\ k = a \cdot 2^k$ then we can write this as:
$$ \mathrm{Fib}(n) \equiv \frac{4\ \mathrm{<<}\ n(3+n)}{(4\ \mathrm{<<}\ 2n) - (2\ \mathrm{<<}\ n) - 1} \mathrm{mod}\ (2\ \mathrm{<<}\ n) $$
Observing that $\mathrm{mod}\ (2\ \mathrm{«}\ n)$ can be expressed as the bitwise and (&) of $(2\ \mathrm{«}\ n) - 1$, we reconstruct our original Python program:
def fib(n):
return (4 << n*(3+n)) // ((4 << 2*n) - (2 << n) - 1) & ((2 << n) - 1)
Although it’s curious to find a non-iterative, closed-form solution, this isn’t a practical method at all. We’re doing integer arithmetic with integers of size $O(n^2)$ bits, and in fact, before performing the final bit-wise and, we’ve got an integer that is the first $n$ Fibonacci numbers concatenated together!
