Project Euler 2（1）

Problem 2
フィボナッチ数列の新しい項は、前の2つの項を足して生成される。最初を1,2とすると、最初の10項は、
1, 2, 3, 5, 8, 13, 21, 34, 55, 89
400万を超えない偶数の項の総和を求めよ。
http://projecteuler.net/index.php?section=problems&id=2

そのまま書くと、

from itertools import ifilter, takewhile
def gen_Fibonacci():

    a, b = 1, 2

    while True:

        yield a

        a, b = b, a + b
N = 4 * 10 ** 6

print sum(ifilter(lambda n: n % 2 == 0,

        takewhile(lambda n: n <= N, gen_Fibonacci())))

しかし、この手の数列は一般項を求める方法があるので、それを使いましょう。計算すると、

Fn+2 = Fn+1 + Fn

F0 = 1

F1 = 2

α = (1 - √5) / 2

β = (1 + √5) / 2

Fn = (βn(2 - α) - αn(2 - β)) / (β - α)

ところで、この数列の偶奇は、

奇, 偶, 奇, 奇, 偶, 奇, 奇, …

なので、3項ごとに偶数が現れます。この式でそれはより明確になります。

F_n+3 = 2F_n+1 + F_n

偶数の項は次のように表せます。

a_n = F_3n+1 = (β³ⁿ⁺¹(2 - α) - α³ⁿ⁺¹(2 - β)) / (β - α)

第2項は小さいので、

N = 4000000
(β³ⁿ⁺¹(2 - α)) / (β - α) < N + 0.5

両辺の対数を取って、最大のnを求めることができます。
そのnまでの総和は、

1 + x + ... + xⁿ = (xⁿ⁺¹ - 1) / (x - 1)

を使って求めます。

from math import sqrt, log
N = 4 * 10 ** 6

a = (1 - sqrt(5)) / 2

b = (1 + sqrt(5)) / 2

f0, f1 = 1, 2

c = b * (f1 - f0 * a) / (b - a)

m_max = int(log( (N + 0.5) / c) / (3 * log(b)))

print int(c * (b ** (3 * m_max + 3) - 1) / (b ** 3 - 1) + 0.5)