フィボナッチ数列 - inamori’s diary

フィボナッチ数列F_nは次のように定義されます。

F₁ = 1
F₂ = 1
F_n+2 = F_n+1 + F_n

F₂ = 2とする流儀もあるようです。
フィボナッチ数列はいろいろなところに出てくるためか、Project Eulerでもしばしば取り上げられます。

フィボナッチ数列をPythonで生成するには、F₀ = 0と考えて、

def gen_fib():
    a, b = 0, 1
    while True:
        yield b
        a, b = b, a + b

n項を求めるにこれを使うと、

N = 100000
print next(islice(gen_fib(), N - 1, None))

n項を直接求めるときは行列を使うと速くなります。定義から、

$A = \begin{pmatrix} 0 & 1 \\ 1 & 1 \end{pmatrix}$
$\begin{pmatrix} F_{n+1} \\ F_{n+2} \end{pmatrix} = A\begin{pmatrix} F_n \\ F_{n+1} \end{pmatrix}$
$\begin{pmatrix} F_{n-1} \\ F_{n} \end{pmatrix} = A^{n-1}\begin{pmatrix} 0 \\ 1 \end{pmatrix}$

となるので、A^n-1をバイナリ法で求めればよいです。

def mul(A, B):
    def f(i, k):
        return sum(A[i][j] * B[j][k] for j in xrange(nb))
    
    na = len(A)
    nb = len(B)
    mb = len(B[1])
    return [ [ f(i, k) for k in xrange(mb) ] for i in xrange(na) ]

def square(A):
    a, b, d = A[0][0], A[0][1], A[1][1]
    a_ = a * a + b * b
    b_ = b * (a + d)
    return [ [ a_, b_ ], [ b_, a_ + b_ ] ]

def pow(A, e):
    if e == 0:
        return [ [ 1 if i == k else 0 for k in xrange(len(A)) ]
                                            for i in xrange(len(A)) ]
    
    B = pow(A, e / 2)
    C = square(B)
    return C if e % 2 == 0 else mul(C, A)

def fib(n):
    A = [ [ 0, 1 ], [ 1, 1] ]
    B = pow(A, n / 2)
    if n % 2 == 1:
        return B[0][1] * B[1][0] + B[1][1] * B[1][1]
    else:
        return B[1][0] * (B[0][0] + B[1][1])

N = 1000000
print fib(N) & 1023

実際には、最後に使う行列の要素は一つだけでよいので、そこは余分な計算をしないように工夫しています。また、2×2行列の乗算は要素の乗算が8回要りますが、自乗の場合対称行列あることを考慮すると4回になります。あと、行列の要素にa₂₂ = a₁₁ + a₁₂という関係があることを考慮すると、3回になります。計算時間は大きな桁の乗算が支配的なので、最後の乗算の計算時間によることになります。nが奇数の場合は2回計算しているので、偶数のときより時間がかかります。

一般式を使う方法もあります。

F_n = (βⁿ - αⁿ) / √5
β = (1 + √5) / 2
α = (1 - √5) / 2

ここで、(a + b√5) / 2を(a, b)という整数の組合せと考えて、

(a, b) * (c, d) = ((ac + 5bd) / 2, (ad + bc) / 2)

という乗算を定義してべき乗を計算します。

def mul(a, b):
    x1, y1 = a
    x2, y2 = b
    return ((x1 * x2 + 5 * y1 * y2) / 2, (x1 * y2 + x2 * y1) / 2)

def square(a):
    x, y = a
    return ((x * x + 5 * y * y) / 2, x * y)

def pow(a, e):
    if e == 0:
        return (2, 0)
    
    b = pow(a, e / 2)
    c = square(b)
    return c if e % 2 == 0 else mul(c, a)

def fib(n):
    a = (1, 1)
    b = pow(a, n / 2)
    if n % 2 == 0:
        return b[0] * b[1]
    else:
        b = pow(a, n)
        return b[1]

N = 1000000
print fib(N) & 1023

偶数ならさきほどより少し速いですが、奇数なら少し遅いという結果になりました。
自乗を求めるときに乗算を3回行っていましたが、

def square(a):
    x, y = a
    return ((x * x + 5 * y * y) / 2, x * y)

よく考えると2回で済みますね。

def square(a):
    x, y = a
    p = (x + 5 * y) * (x + y)
    q = x * y
    return ((p - 6 * q) / 2, q)

これで少し速くなります。

Project Eulerに必要な用語集