# Generalized Project Euler 2: A sledgehammer to crack a nut

I have previously asked about how optimizing Project Euler 1. See Generalized Project Euler 1: A sledgehammer to crack a nut. In this question I wanted to explore the next problem: Project Euler 2: Even Fibonacci numbers.

Even Fibonacci numbers

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, \ldots$$

By considering the terms in the Fibonacci sequence whose values do not exceed four million, find the sum of the even-valued terms.

## My attempt

As in standard literature I define the fibonacci sequence as

$$F_n = F_{n-1} + F_{n-2}$$

with initial values $F_0 = 0$, $F_1 = 1$. It can be shown that every even fibonacci number is on the form $F_{3k}$, where $k$ is some positive integer. Let $E_n = F_{3n}$, denote the $n$th even Fibonacci number. We then have the following reccurence

$$E_n = 4 E_{n-1} + E_{n-2}$$

with initial values $E_0 = 0$ and $E_1 = 2$. The sum of the first $n$ even fibonacci numbers has a closed form, and can be written as \begin{align*} F_0 + F_3 + \cdots + F_{3n-1} + F_{3n} & = \frac{F_{3n+2}-1}{2} \\ E_0 + E_1 + \cdots + E_{n-1} + E_{n} & = \frac{E_{n} + E_{n+1}-2}{4} \end{align*} using the notation defined earlier. The two expressions above are equal, and just differ in notation used. Next we can find the index of the largest even Fibonacci number not exceeding $M$.

$$F_{3n} = \text{round}\left[\frac{\phi^{3n}}{\sqrt{5}}\right] = M \quad \Rightarrow \quad n = \frac{1}{6} \log_\phi 5 + \frac{1}{3} \log_\phi M\,,$$

where $M$ is the golden ratio. So $n$ rounded down gives us how many terms we need to iterate over.

## The struggle

So the problem boils down to either computing $F_{3n+2}$ fast, or computing $E_n + E_n+2$ fast. I found the following paper A fast algorithm for computing large Fibonacci numbers - Daisuke Takahashi. This explains a method to compute the $n$th Fibonacci number in logarithmic time. Using this I was able to write a quick and fast solution to the problem. An alternative is to use matrix multiplication to find the $n$th. However as my benchmarks show this is slightly slower than Daisuke Takahashi optimized algorithm. A third option is to use the relationship between the Lucas numbers and Fibonacci to generate the n'th term quickly. I will use that idea, but use the more generall Lucas sequences instead.

Another solution is to use the relationship between the Fibonacci numbers and the Lucas

Another solution I thought of was to use the Lucas sequences.

\begin{align} U_0(P,Q)&=0, \\ U_1(P,Q)&=1, \\ U_n(P,Q)&=P\cdot U_{n-1}(P,Q)-Q\cdot U_{n-2}(P,Q) \mbox{ for }n>1, \end{align}

and

\begin{align} V_0(P,Q)&=2, \\ V_1(P,Q)&=P, \\ V_n(P,Q)&=P\cdot V_{n-1}(P,Q)-Q\cdot V_{n-2}(P,Q) \mbox{ for } n>1. \end{align}

For example the Fibonacci sequence can be represented as $F_n = U_n(1,-1)$. The idea is to introduce a second sequence $V_n(P, Q)$ along side $U_n(P, Q)$ such that we can compute $U_n$ faster, than by the standard recursion: $U_n = P U_{n-1} - Q U_{n-2}$.

First, we can double the subscript from $k$ to $2k$ in one step using the recurrence relations

\begin{align*} U_{2k} & = U_k\cdot V_k \\ V_{2k} & = V_k^2-2Q^k . \end{align*} Next, we can increase the subscript by $1$ using the recurrences \begin{align*} U_{2k+1} & = (P\cdot U_{2k} + V_{2k})/2 \\ V_{2k+1} & = (D\cdot U_{2k} + P\cdot V_{2k})/2. \end{align*}

Using this square and add algorithm, we can evaluate $F_{44}$ with the following computations: \begin{align*} F_1, \, F_2, \, F_4, \, F_5, \, F_{10}, \, F_{11}, \, F_{22}, \, F_{44} \end{align*} The benefit of using the generalized lucas sequences instead of the lucas numbers is that we can now iterate over $E_n = 2\cdot U_{n}(4,-1)$ instead of the slower $F_{3n+2} = U_{n}(1,-1)$.

## Implementation and question

Using the theory above I have implemented three functions to solve this problem

• Daisuke Takahashi - Finding $F_{3n+2}$
• Lucas sequence finding $F_{3n+2}$
• Lucas sequence finding $E_n$ with $2 \cdot U_{4,-1}$ and then use this to find the next term $E_{n+1}$.

I have done all these implentations in the code below, and I do have some questions. Right now I do not have time to include benchmarks, but there is something strange with the running times of the three implementations.

• Why is Takadashi faster than both my implentations? I can fully understand why it beats finding $F_{3n+2}$ using the generalized Lucas sequence. However finding $E_n$ should be quite a bit quicker, why is it not?
• I am not happy with how I have implemented the generalized lucas sequence.
• Is there a faster way to find $E_n$ and $E_{n+1}$ than with the generalized Lucas sequences?

Please no comments on the Takadashi implementation. I included it to benchmark my use of the generalized lucas sequences.

## Code

from math import log

PHI = (1 + 5**0.5)/float(2)
LOG_5 = log(5, PHI)/float(6)

def largest_even_fib_under_n(limit):
'''
The nth even fibonacci number can be written as
F_{3n} = round( phi^3k / sqrt 5 )
This function solves the equation F_{3n} = limit for n.
'''
return int(LOG_5 + log(limit, PHI)/float(3))

def lucas_sequence(p, q, n):
'''
Calculates the n'th term of the lucas sequences U_n and V_n.
https://en.wikipedia.org/wiki/Lucas_sequence#Algebraic_relations
U_0 = 0
U_1 = 1
U_n = P * U_(n-1) - Q * U_(n-2)
and
V_0 = 2
V_1 = P
V_n = P * V_(n-1) - Q * V_(n-2)
Uses the doubling relations to compute it in log n time.
https://en.wikipedia.org/wiki/Lucas_pseudoprime#Implementing_a_Lucas_probable_prime_test
U_2k = U_k * V_k
V_2k = V_k^2 - 2*Q^k
and
U_(2k+1) = (P * U_2k + V_2k)/2
V_(2k+1) = (D * U_2k + P * V_2k)/2
If the index is odd, we can make it even using the last relation, and then square.
'''
d = p * p - 4 * q

un, vn, qn = 1, p, q
u = 0 if n % 2 == 0 else 1
v = 2 if n % 2 == 0 else p
k = 1 if n % 2 == 0 else q

n = n // 2
while n > 0:

u2 = un * vn
v2 = vn * vn - 2 * qn
q2 = qn * qn

n2 = n // 2
if n % 2 == 1:

uu = (u * v2 + u2 * v)/2
vv = (v * v2 + d * u * u2)/2
u, v, k = uu, vv, k * q2

un, vn, qn, n = u2, v2, q2, n2

return u, v

def sum_fibonacci_with_lucas(limit):
n = largest_even_fib_under_n(limit)
U = lucas_sequence(1, -1, 3*n+2)
return (U[0] - 1)/2

def sum_even_fibonacci_with_lucas(limit):
n = largest_even_fib_under_n(limit)
p = 4
q = -1

u0, v0 = lucas_sequence(p, q, n)

u0 *= 2
u1 = (p*u0 + 2*v0)/2

return (u0+u1-2)/4

def fib_takahashi(n):
'''
A fast algorithm for computing large Fibonacci numbers by Daisuke Takahashi
http://www.math.tamu.edu/~snpolloc/math491_691/Takahashi00.pdf
'''
if n == 0:
return n
F, L, sign, exp = 1, 1, -2, int(log (n,2))
for i in xrange (exp - 1):
F2   = F**2
FL2  = (F + L)**2
F    = ((FL2 - 6*F2) >> 1) - sign
temp = (FL2 >> 2 ) + F2
L     = temp + (F << 1)
F     = temp
else:
L    = 5*F2 + sign
sign = -2 if n & mask else 2
if n & (mask >> 1) == 0:
return F * L
else:
return ((F + L) >> 1) * L - (sign >> 1)

def sum_fibonacci_takahashi(limit):
n = largest_even_fib_under_n(limit)
return (fib_takahashi(3*n+2) - 1)/2

if __name__ == '__main__':
limit = 4*10**6
print sum_fibonacci_with_lucas(limit)
print sum_even_fibonacci_with_lucas(limit)
print sum_fibonacci_takahashi(limit)

import timeit
times = 100
t1 = timeit.timeit("sum_fibonacci_with_lucas(4*10**(10**6))",
setup="from __main__ import sum_fibonacci_with_lucas", number=times)/float(times)
t2 = timeit.timeit("sum_even_fibonacci_with_lucas(4*10**(10**6))",
setup="from __main__ import sum_even_fibonacci_with_lucas", number=times)/float(times)
t3 = timeit.timeit("sum_fibonacci_takahashi(4*10**(10**6))",
setup="from __main__ import sum_fibonacci_takahashi", number=times)/float(times)
print'''
sum using lucas 3n+2 used {:10.7f} ms
sum using lucas  n   used {:10.7f} ms
sum using takahashi  used {:10.7f} ms

takahashi was {:10.7f} times faster than lucas  n
takahashi was {:10.7f} times faster than lucas 3n+2
lucas  n  was {:10.7f} times faster than lucas 3n+2
'''.format(1000*t1, 1000*t2, 1000*t3, t1/float(t3), t2/float(t3), t1/float(t2))

• Did you ignore the PEP8 on purpose? – Mast Jul 16 '16 at 12:08

A few small things:

• u = 0 if n % 2 == 0 else 1 => u = n % 2
• k = 1 if n % 2 == 0 else q => k = 1 + n%2 * (q-1)
• v = 2 if n % 2 == 0 else p => v = 2 + n%2 * (p-2)
• u1 = (p*u0 + 2*v0)/2 => u1 = p*u0/2 + v0

These are all quite small optimizations, but they are both purely faster.

I just realized that you are missing the obvious answer. As long as x < 70, floats are accurate enough that:

def quick_fib(n):
return round(PHI**(n)/5**.5)


works correctly. This method is much faster, and more importantly (arguably) much simpler.