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, ...

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



We could do this one easily by brute force, but I sense there's a closed-form solution.

I remember there's a closed form expression for the $$n^\text{th}$$ Fibonacci number: http://en.wikipedia.org/wiki/Fibonacci_number#Closed-form_expression

Also note that the Fibonacci numbers go $$o, e, o, o, e, o, o, e, o, o, ...$$

where o/e means odd/even. So we want the ones at positions 2, 5, 8, 11, (...)--so it'll look something like:

$$ \sum_{i=1}^{u} F_{3i - 1} $$

where $u$ is some upper limit.

The closed form $i^{\text{th}}$ Fibonacci number is:

$$F_{i} = \frac{(\frac{1+\sqrt{5}}{2})^{i} - (\frac{1-\sqrt{5}}{2})^{i}}{\sqrt{5}}$$

In [3]:
# The values of a and b:
a = (1.0 + sqrt(5.0))/2.0
b = (1.0 - sqrt(5.0))/2.0
print "a =",a, "b =",b, "\n"

def fib_ith(i):
    return (a**i - b**i)/sqrt(5.0)

#for i in range(1,10):
#    print fib_ith(i)
def fib_even_sum_brute(N):
    s = 0
    i = 3
    i_f = fib_ith(i)
    while i_f < N:
        s += i_f
        #print "i =", i, "F_i =", i_f, "sum =", s
        i += 3
        i_f = fib_ith(i)
    #print "At i =", i-3, "Sum: ",s
    return floor(s)
fib_even_sum_brute(N=4e6)


a = 1.61803398875 b = -0.61803398875 

Out[3]:
4613732.0

Which is the right answer for this challenge. But let's find the closed form.

For sums of powers (finite geometric series), we have: $$\sum_{i=1}^{n} c^{i} = \frac{c (c^n-1)}{c-1} $$

The Fibonacci sequence as in closed form by $F_{i}$ above starts with 0, 1; so we actually want the positions 3,6,9,...

I.e. we want: $$ \sum_{i=1}^{u} F_{3i} \\ = \sum_{i=1}^{u} \frac{(\frac{1+\sqrt{5}}{2})^{3i} - (\frac{1-\sqrt{5}}{2})^{3i}}{\sqrt{5}} \\ = \frac{1}{\sqrt{5}}\sum_{i=1}^{u} (\frac{1+\sqrt{5}}{2})^{3i} - \frac{1}{\sqrt{5}}\sum_{i=1}^{U} (\frac{1-\sqrt{5}}{2})^{3i} $$

Now we have to do a little trickery as the formula above for a finite geometric series isn't quite in the form we want.

$$\sum_{i=1}^{u} c^{3i} = (1 - c^{3})\sum_{i=1}^{u} c^{3i} \\ = \sum_{i=1}^{u} c^{3i} - (c^{3u + 3} + \sum_{i=2}^{u} c^{3i}) \\ = c^{3} - c^{3u + 3} $$

and so:

$$ \sum_{i=1}^{u} c^{3i} = \frac{c^{3u + 3} - c^{3}}{c^{3} - 1} $$

Thus we can now express the Fibonacci sum as:

$$ \frac{1}{\sqrt{5}}(\frac{a^{3u + 3} - a^{3}}{a^{3} - 1} - \frac{b^{3u + 3} - b^{3}}{b^{3} - 1}) $$

where $a = \frac{1+\sqrt{5}}{2}$ and $b = \frac{1-\sqrt{5}}{2}$

The next question is what that upper limit $u$ is.


In [19]:
import numpy as np
import matplotlib.pyplot as plt

t = np.arange(1., 5., 1)

# red dashes, blue squares and green triangles
plt.plot(t, a**t, 'r--')
plt.plot(t, b**t, 'b--')


Out[19]:
[<matplotlib.lines.Line2D at 0x48ed6d0>]

From the above we can see that the second term involving $b$ goes to zero pretty quickly--so for a large $u$ we can get lucky most likely and there's an easy solution.

$$\frac{(\frac{1+\sqrt{5}}{2})^{i} - (\frac{1-\sqrt{5}}{2})^{i}}{\sqrt{5}} \\ \approx \frac{1}{\sqrt{5}}(\frac{1+\sqrt{5}}{2})^{i} $$

And so $u = \lfloor \frac{1}{3} \log_{a} N + \frac{1}{3} \log_{a}{\sqrt 5} \rfloor = \lfloor \frac{1}{3}\frac{\log N}{\log a} + \frac{1}{3} \frac{\log \sqrt 5}{\log a}\rfloor$ ought to do the trick.


In [4]:
def fib_even_sum(N):
    a = (1.0 + sqrt(5.0))/2.0
    b = (1.0 - sqrt(5.0))/2.0
    
    u = floor(log(N)/(3.0*log(a)) + log(sqrt(5))/(log(a)*3))
    
    return floor(((a**(3.0*(u+1)) - a**3.0)/(a**3.0 - 1) -
            (b**(3.0*(u+1)) - b**3.0)/(b**3.0 - 1))/sqrt(5.0))

In [3]:
fib_even_sum(N=4e6)


Out[3]:
4613732.0

In [17]:
fib_even_sum_brute(N=4e6)


Out[17]:
4613732.0

I.e. it works. Let's do some timing just to halfway rationalize that extra neural expenditure:


In [14]:
import timeit
print "Time for 10,000:", timeit.timeit('fib_even_sum(N=4e6)', number=10000, setup="from __main__ import fib_even_sum"), "s"


Time for 10,000: 0.266054153442 s

In [13]:
import timeit
print "Time for 10,000:", timeit.timeit('fib_even_sum_brute(N=4e6)', number=10000, setup="from __main__ import fib_even_sum_brute"), "s"


Time for 10,000: 0.9137840271 s

Woo!