Project Euler(2)

プロジェクトオイラー
http://projecteuler.net/index.php

Q2.
フィボナッチ数列の4000000を超えない偶数項の和。


from itertools import takewhile

def fib():
a, b = 1, 2
while True:
if a & 1 == 0:
yield a
a, b = b, a + b

N = 4000000
print reduce(lambda x, y: x + y, takewhile(lambda x: x <= N, fib()))

ここでのフィボナッチ数の一般項は、


α = (1 - √5) / 2
β = (1 + √5) / 2
fn = (βn - αn) / √5

この一般項を、


fn = a0a1n + a2a3n

と書く。
偶数項は、f2、f5のように3で割って2余る項だから、


gn = f3n-1

だから、これも同じように一般項を書ける。
これを使って、和も書ける。


from math import sqrt, log

alpha = (1 - sqrt(5)) / 2
beta = (1 + sqrt(5)) / 2
a = [ 1 / sqrt(5) * beta, beta, -1 / sqrt(5) * alpha, alpha ]
b = [ a[0] / a[1], a[1] ** 3, a[2] / a[3], a[3] ** 3 ]

def item(a, n):
return int(a[0] * a[1] ** n + a[2] * a[3] ** n + 0.5)

def f(n):
return item(b, n)

def max_n(f, N):
n = int(log(N / b[0]) / log(b[1]))
if f(n) <= N:
while True:
n += 1
if f(n) > N:
return n - 1
else:
while True:
n -= 1
if f(n) <= N:
return n

def sum(b, n):
n1 = b[0] * b[1] * (b[1] ** n - 1) / (b[1] - 1)
n2 = b[2] * b[3] * (b[3] ** n - 1) / (b[3] - 1)
return int(n1 + n2)

N = 4000000
n = max_n(f, N)
print sum(b, n)