快速斐波那契计算

Question

几周前，我在Google+上看到一条评论，其中有人演示了一种直接计算Fibonacci数的方法，这种方法不基于递归，也不使用记忆。他实际上只是记住了最后两个数字，并不断地将它们相加。这是一个O（n）算法，但他实现得非常干净。所以我很快指出，一个更快的方法是利用它们可以计算为[[0,1]，[1,1]]矩阵的幂，而且只需要O（log（N））计算。在

当然，问题是，这远远不是超过某一点的最佳值。只要数值不太大，它是有效的，但长度以N*log（phi）/log（10）的速率增长，其中N是第N个斐波纳契数，phi是黄金比率（1+sqrt（5））/2~1.6。事实证明，log（phi）/log（10）非常接近1/5。所以第N个Fibonacci数大概有N/5位数。在

矩阵乘法，甚至是数字乘法，当数字开始有数百万或数十亿位数时，变得非常缓慢。所以F（100000）计算大约需要0.03秒（在Python中），而F（1000000）大约需要5秒。这几乎不是O（logn）增长。我的估计是，这种方法在没有改进的情况下，只将计算优化为O（（log（N））^（2.5））左右。在

以这样的速度计算第十亿个斐波纳契数会非常慢（尽管它只有大约1000000000/5位数字，所以很容易放入32位内存中）。在

有没有人知道一种实现或算法可以使计算速度更快？也许可以计算出万亿次的斐波纳契数。在

为了清楚起见，我不是在寻找一个近似值。我在寻找精确的计算结果（精确到最后一位）。在

编辑1:我正在添加Python代码，以显示我认为是O（（logn）^2.5））算法。在

from operator import mul as mul
from time import clock

class TwoByTwoMatrix:
    __slots__ = "rows"

    def __init__(self, m):
        self.rows = m

    def __imul__(self, other):
        self.rows = [[sum(map(mul, my_row, oth_col)) for oth_col in zip(*other.rows)] for my_row in self.rows]
        return self

    def intpow(self, i):
        i = int(i)
        result = TwoByTwoMatrix([[long(1),long(0)],[long(0),long(1)]])
        if i <= 0:
            return result
        k = 0
        while i % 2 == 0:
            k +=1
            i >>= 1
        multiplier = TwoByTwoMatrix(self.rows)
        while i > 0:
            if i & 1:
                result *= multiplier
            multiplier *= multiplier # square it
            i >>= 1
        for j in xrange(k):
            result *= result
        return result


m = TwoByTwoMatrix([[0,1],[1,1]])

t1 = clock()
print len(str(m.intpow(100000).rows[1][1]))
t2 = clock()
print t2 - t1

t1 = clock()
print len(str(m.intpow(1000000).rows[1][1]))
t2 = clock()
print t2 - t1

编辑2: 看起来我没有考虑到len(str(...))会对测试的整个运行时做出重大贡献这一事实。将测试更改为

^{pr2}$

将运行时间缩短到.008秒和.31秒（当使用len(str(...))时，从.03秒和5秒缩短到.008秒和.31秒）。在

因为M=[[0,1]，[1,1]]的幂N是[[F（N-2），F（N-1）]，[F（N-1），F（N）]， 另一个明显的效率低下的原因是计算矩阵的（0,1）和（1,0）元素，好像它们是不同的一样。这个（我换成了Python3，但是Python2.7倍是相似的）：

class SymTwoByTwoMatrix():
    # elments (0,0), (0,1), (1,1) of a symmetric 2x2 matrix are a, b, c.
    # b is also the (1,0) element because the matrix is symmetric

    def __init__(self, a, b, c):
        self.a = a
        self.b = b
        self.c = c

    def __imul__(self, other):
        # this multiplication does work correctly because we 
        # are multiplying powers of the same symmetric matrix
        self.a, self.b, self.c = \
            self.a * other.a + self.b * other.b, \
            self.a * other.b + self.b * other.c, \
            self.b * other.b + self.c * other.c
        return self

    def intpow(self, i):
        i = int(i)
        result = SymTwoByTwoMatrix(1, 0, 1)
        if i <= 0:
            return result
        k = 0
        while i % 2 == 0:
            k +=1
            i >>= 1
        multiplier = SymTwoByTwoMatrix(self.a, self.b, self.c)
        while i > 0:
            if i & 1:
                result *= multiplier
            multiplier *= multiplier # square it
            i >>= 1
        for j in range(k):
            result *= result
        return result

在9.51秒内计算出F（100000）in.006，F（1000000）in.235和F（10000000）。在

这是意料之中的。它产生的结果比最快的测试快45%，预期增益应该渐近接近 功率因数/（1+2*功率因数+功率因数*功率因数）~23.6%。在

M^N的（0,0）元素实际上是N-2nd斐波纳契数：

for i in range(15):
    x = m.intpow(i)
    print([x.a,x.b,x.c])

给予

[1, 0, 1]
[0, 1, 1]
[1, 1, 2]
[1, 2, 3]
[2, 3, 5]
[3, 5, 8]
[5, 8, 13]
[8, 13, 21]
[13, 21, 34]
[21, 34, 55]
[34, 55, 89]
[55, 89, 144]
[89, 144, 233]
[144, 233, 377]
[233, 377, 610]

我认为不必计算元素（0,0）将产生额外的1/（1+phi+phi*phi）~19%的加速。但F（2N）和F（2N-1）solution given by Eli Korvigo below的lru_cache实际上给出了4倍的加速（即75%）。所以，虽然我还没有找到一个正式的解释，但我很想认为它在N的二进制展开中缓存了1的跨距，并且做了必要的最小数量的乘法。这就避免了寻找这些范围的需要，预先计算它们，然后在N的展开中的正确点乘以它们。lru_cache允许从上到下的计算，这将是一个更复杂的按钮到顶部的计算。在

当N增长10倍时，SymTwoByTwoMatrix和lru缓存-of-F（2N）-和-F（2N-1）都要花费大约40倍的时间来计算。我认为这可能是由于Python实现了长整型的乘法。我认为大数的乘法和加法应该是可并行化。因此，多线程sub-O（N）解决方案应该是可能的，即使（正如Daniel Fisher在评论中所说）F（N）解决方案是Theta(n)。在

Answer 1

从Wikipedia

对于所有n≥0，Fn是最接近phi^n/sqrt（5）的整数，其中phi是黄金比率。因此，可以通过四舍五入，即使用最近的整数函数来找到它

Answer 2

由于Fibonacci序列是一个线性递归序列，它的成员可以用封闭形式来计算。这涉及到计算一个幂，这可以在O（logn）中完成，类似于矩阵乘法解决方案，但是常量开销应该更低。这是我知道的最快的算法。在

编辑

对不起，我错过了“确切”部分。矩阵乘法的另一个精确的O（log（n））备选方案可计算如下

from functools import lru_cache

@lru_cache(None)
def fib(n):
    if n in (0, 1):
        return 1
    if n & 1:  # if n is odd, it's faster than checking with modulo
        return fib((n+1)//2 - 1) * (2*fib((n+1)//2) - fib((n+1)//2 - 1))
    a, b = fib(n//2 - 1), fib(n//2)
    return a**2 + b**2

这是基于Edsger Dijkstra教授从a note派生出来的。该解决方案利用了这样一个事实：要同时计算F（2N）和F（2N-1），只需要知道F（N）和F（N-1）。尽管如此，您仍然在处理长数算法，尽管开销应该小于基于矩阵的解决方案的开销。在Python中，您最好用命令式的方式重写它，因为它的记忆和递归速度很慢，不过我这样写是为了使函数公式更清晰。在

Answer 3

使用另一个答案中奇怪的平方根方程closed form fibo你可以精确地计算出第k个斐波纳契数。这是因为$\sqrt（5）$最终落空。你只需安排好乘法运算就可以同时跟踪它。在

def rootiply(a1,b1,a2,b2,c):
    ''' multipy a1+b1*sqrt(c) and a2+b2*sqrt(c)... return a,b'''
    return a1*a2 + b1*b2*c, a1*b2 + a2*b1

def rootipower(a,b,c,n):
    ''' raise a + b * sqrt(c) to the nth power... returns the new a,b and c of the result in the same format'''
    ar,br = 1,0
    while n != 0:
        if n%2:
            ar,br = rootiply(ar,br,a,b,c)
        a,b = rootiply(a,b,a,b,c)
        n /= 2
    return ar,br

def fib(k):
    ''' the kth fibonacci number'''
    a1,b1 = rootipower(1,1,5,k)
    a2,b2 = rootipower(1,-1,5,k)
    a = a1-a2
    b = b1-b2
    a,b = rootiply(0,1,a,b,5)
    # b should be 0!
    assert b == 0
    return a/2**k/5

if __name__ == "__main__":
    assert rootipower(1,2,3,3) == (37,30) # 1+2sqrt(3) **3 => 13 + 4sqrt(3) => 39 + 30sqrt(3)
    assert fib(10)==55