Wednesday, June 3, 2009

Exploring the Fibonacci Numbers

Recently I decided to teach myself Haskell by using it to solve Project Euler problems. This kills two birds with one stone, since I've been wanting to tackle both. Also, it sounded like a lot of fun.

Project Euler problem #2 is simple enough:
Each new term in the Fibonacci sequence is generated by adding the previous two terms. By starting with 1 and 2, the first 10 terms will be:

1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...

Find the sum of all the even-valued terms in the sequence which do not exceed four million.
But, since I'm doing this for personal enrichment I decided to investigate in more detail. I did not notice this before solving it the first time (asking a stupid question in the process), but it turns out that every third fibonacci number is even. Furthermore, there is a surprising relationship between every third fibonacci number: f(n) = 4*f(n-3) + f(n-6)[1]. So it is possible to solve the problem efficiently by enumerating just the even fibonacci numbers. In this post I'll use g(n) to denote the nth even fibonacci number. The sequence can be implemented in Haskell, using list comprehension syntax, like so:

-- g(n) = 4*g(n-1) + g(n-2)
g :: [Integer]
g = 2 : 8 : [ a + 4*b | (a,b) <- zip g (tail fibs2)]

But can we do better? I remembered one of Raganwald's homoiconic posts about computing fibonacci numbers efficiently using a matrix multiplication technique (pardon the ascii-art notation):

m^n, where m = (1 1) , is equivalent to: ( f(n) f(n-1) )
(1 0) ( f(n-1) f(n-2) )

So the nth fibonacci number is the upper-left element in the matrix obtained by multiplying m by itself n times[2].

Also interesting is that a source article of the post at homoiconic gives an efficient algorithm for calculating the sum of fibonacci numbers using the matrix representation of the sequence. Now we're getting somewhere! It should be possible to modify the matrix approach so that it generates only the even fibonacci numbers, using the formula described above[3]. And if we can do that, we can calculate their sum efficiently with the provided algorithm.

You can do exactly that using these two matrices:

b = (8 2) m = (4 1)
(2 0) (1 0)
where b is the base case and m is the multiplier. Notice that if you multiply a generic fibonacci matrix (where a,b,c are consecutive even fibonaccis) by m, you get:

(a b) * (4 1) = (4a+b, a)
(b c) (1 0) (4b+c, b)
Which generates the matrix containing the next consecutive even fibonacci number (recall that the nth even fib g(n) is equal to 4*g(n-1) + g(n-2), which is 4a+b if a and b are the two previous even fibs). Using the matrix b to bootstrap the multiplication (because it contains 2 and 8, the first two even fibs), it follows that b*m should produce the 3rd even fibonacci number, b*m*m should produce the fourth, b * m^3 the fifth, and in general:

g(1) = 2
g(n) = b * m^(n-2) [4]
So, now we have an equation for using matrix multiplication to compute even fibonacci numbers. How can we use it to efficiently compute the sum of the first n numbers in that sequence? Let's use a concrete example, n=7. We want to calculate

g(1) + g(2) + ... + g(7) = 2 + b + b*m + b*m^2 + ... + b*m^5
factoring out the common b term, we have

2 + b(m^0 + m^1 + m^2 + m^3 + m^4 + m^5)
where m^0 evaluates to the identity matrix. For convenience, let's define a function sum(n,m) that computes the sum m^0 + ... + m^n, so the above equation can be written 2 + b*sum(5,m) We can see that the bulk of the computation for large values of n will be in calculating each of the matrices m^2, m^3,...,m^n. But notice the following equivalence:

m^0 + m^1 + ... + m^5 = (m^0 + m^1 + m^2) + m^3(m^0 + m^1 + m^2)
by expanding the equation in this manner, we can calculate sum(5,m) with considerably less work. Notice the recursive structure of the equation: sum(5,m) = sum(2,m) + m^3 * sum(2,m). In general pseudocode:

sum(n,m) =
if n == 0 then return i //the identity matrix
if odd(n) then let s = sum(floor(n/2), m)
return s + m^ceil(n/2) * s
if even(n) then let s = sum(n/2, m)
return s + m^(n/2) * s + m^n

the 'even' case is slightly more complicated because it needs to take into account that the sequence m^0 + ... m^n is not evenly divisible into two pieces. The algorithm can be further optimized by having our sum function return m^n as well as the current running sum. This unfortunately obfuscates the algorithm with extra bookkeeping, but it ensures that we aren't wasting any work from one call to the next, and sum(n,m) must calculate m^n or m^n/2 anyway, so it is not much overhead. The 'odd' step from the above algorithm would then look something like this:

let s, m_exp = sum(floor(n/2), m) //m_exp is m^floor(n/2)
return s + (m_exp * m) * s, m_exp*m_exp*n //return both the sum and m^n
One final optimization. Raganwald noted in his original post that there is redundancy in the matrix representation, and we can easily re-write our matrices as tuples of 3 elements:

b = (8, 2, 0) m = (4, 1, 0)
Multiplication on these tuples is essentially the same as matrix multiplication (details at homoiconic). Putting it all together, we have the final, executable haskell:

The naive implementation is faster for small n, but once you get above n=5-10k this optimized version overtakes it and never looks back. The full code, including multiple solutions to the project Euler question and benchmarking of the naive implementation against this one, is available here.

[1] This is easy to verify yourself. Start with 4*f(n-3) + f(n-6) and keep simplifying until you get f(n)
[2] For convenience, define the first two fibonacci numbers to be 1,2. Take some time to convince yourself of this, if you like. The computation can be done in O(log n) multiplications (instead of n multiplications) by multiplying in a clever order. See the source articles for details.
[3] g(n) = 4*g(n-1) + g(n-2)
[4] I'm abusing notation slightly here. Technically this is the 2x2 matrix containing g(n), not g(n) itself. Also, note that when n=2, m^0 evaluates to the identity matrix. b * m^0 = b * i = b, so no special case is needed for 8.


Frank said...

(Note, I'm using F[0] = 0, F[1] = 1, F[2] = 1, F[3] = 2, ...
Thus indicies which are multiples of three give even Fibonacci numbers.)

By the recursion relation:
Sum[F[3*i],{i,0,n}] = Sum[F[3*i-1] + F[3*i-2],{i,1,n}] = Sum[F[i],{i,0,3*n}] / 2.

From the analytical form, F[n] = (phi^n - (-1/phi)^n)/sqrt(5), you can derive:
Sum[F[i],{i,0,n}] = F[n+2] - 1.

Then for large enough arguments, you can approximate F[n] as phi^n / sqrt(5):

phi = (1 + sqrt(5))/2
f = 4000000
n = floor(log(sqrt(5)*f)/log(phi))
n = n - n mod 3
ans = (floor(phi^(n+2) / sqrt(5) + 1/2) - 1) / 2

Anyway, thanks for the interesting post. I'm going to have to give Haskell a shot. Hope you have fun with Project Euler.

Anonymous said...

(floor (fromIntegral n/2))

should be written as

(n `div` 2)