The task at handI 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 numbersEach 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 attemptAs 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 struggleSo 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 questionUsing the theory above I have implemented three functions to solve this problemDaisuke 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.Codefrom math import logPHI = (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, vdef sum_fibonacci_with_lucas(limit):    n = largest_even_fib_under_n(limit)    U = lucas_sequence(1, -1, 3*n+2)    return (U[0] - 1)/2def 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)/4def 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))    mask = 2**exp    for i in xrange (exp - 1):        mask = mask >> 1        F2   = F**2        FL2  = (F + L)**2        F    = ((FL2 - 6*F2) >> 1) - sign        if n & mask:            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)/2if __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} mssum using lucas  n   used {:10.7f} mssum using takahashi  used {:10.7f} mstakahashi was {:10.7f} times faster than lucas  n takahashi was {:10.7f} times faster than lucas 3n+2lucas  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))
Generalized Project Euler 2: A sledgehammer to crack a nut
python;performance;python 2.7;mathematics