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: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:
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.
-- 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)
(a b) * (4 1) = (4a+b, a)
(b c) (1 0) (4b+c, b)
g(1) = 2
g(n) = b * m^(n-2) [4]
g(1) + g(2) + ... + g(7) = 2 + b + b*m + b*m^2 + ... + b*m^5
2 + b(m^0 + m^1 + m^2 + m^3 + m^4 + m^5)
m^0 + m^1 + ... + m^5 = (m^0 + m^1 + m^2) + m^3(m^0 + m^1 + m^2)
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
b = (8, 2, 0) m = (4, 1, 0)
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.
