快速斐波那契计算

Question

几周前我在Google+上看到了一条评论，其中有人展示了斐波那契数字的直接计算，这些数字并非基于递归而且没有使用记忆。他实际上只记得最后两个数字并不断添加它们。这是一个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。因此，预计Nth Fibonacci数字大约为N / 5位数。

当数字开始有数百万或数十亿的数字时，矩阵乘法，即使偶数乘法，也会非常慢。因此，F（100,000）计算大约0.03秒（在Python中），而F（1000,000）大约需要5秒钟。这几乎不是O（log（N））增长。我的估计是这种方法没有改进，只是将计算优化为O（（log（N））^（2.5））左右。

以这个速率计算第十亿个Fibonacci数将会非常慢（即使它只有〜1,000,000,000 / 5位数，因此很容易适应32位内存）。

有谁知道一个允许更快计算的实现或算法？也许某些东西可以计算出万亿的斐波纳契数。

而且要清楚，我不是在寻找近似值。我正在寻找精确的计算（到最后一位数）。

编辑1：我正在添加Python代码以显示我认为的O（（log N）^ 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(...))会对测试的整体运行时间做出重大贡献这一事实。将测试更改为

from math import log as log

t1 = clock()
print log(m.intpow(100000).rows[1][1])/log(10)
t2 = clock()
print t2 - t1

t1 = clock()
print log(m.intpow(1000000).rows[1][1])/log(10)
t2 = clock()
print t2 - t1

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

因为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

在.006中计算F（100,000），在.235中计算F（1,000,000），在9.51秒内计算F（10,000,000）。

这是可以预料的。对于最快的测试，它产生的结果快45％，并且预期增益应渐近接近phi /（1 + 2 * phi + phi * phi）~23.6％。

M ^ N的（0,0）元素实际上是N-2nd Fibonacci数：

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）lru_cache的solution given by Eli Korvigo below实际上加速了4倍（即75％）。因此，虽然我没有得出正式的解释，但我很想到它在N的二进制扩展中缓存1的跨度并且需要最小的乘法次数。这样就不需要找到那些范围，预先计算它们，然后在N的扩展中将它们乘以正确的点.lru_cache允许从上到下计算本来更复杂的从顶部到顶部的计算。

每当N增长10倍时，SymTwoByTwoMatrix和lru_cache-of-F（2N）-and-F（2N-1）的计算时间大约要长40倍。我认为这可能是由于Python实现了长整数的乘法。我认为大数的乘法和它们的加法应该是可并行的。因此，即使（如Daniel Fisher在评论中所述）F（N）解决方案是Theta(n)，也应该可以实现多线程子O（N）解决方案。

Answer 1

由于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教授从note推导出来的。该解决方案利用了以下事实：要计算F（2N）和F（2N-1），您只需要知道F（N）和F（N-1）。尽管如此，你仍然在处理长数量的算术，尽管开销应该小于基于矩阵的解决方案。在Python中，由于记忆和递归缓慢，你最好用命令式的方式重写它，尽管我这样写它是为了清晰的功能表达。

Answer 2

在另一个答案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

Answer 3

来自Wikipedia，

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

Answer 4

这个评论太长了，所以我会留下答案。

Aaron的回答是正确的，我也赞成你，就像你一样。我将提供相同的答案，并解释为什么它不仅正确，而且到目前为止发布的最佳答案。我们讨论的The formula是：

计算Φ是O(M(n))，其中M(n)是乘法的复杂性（currently略高于线性），而n是位数。

然后有一个幂函数，可以表示为log（O(M(n)•log(n)），乘法（O(M(n))）和exp（O(M(n)•log(n)）。

然后有一个平方根（O(M(n))），一个师（O(M(n))）和一个最后一轮（O(n)）。

这使得这个答案类似于O(n•log^2(n)•log(log(n)))用于n位。

---

我没有彻底分析除法算法，但是如果我正确读取这个，那么每个位可能需要递归（你需要将数字除以log(2^n)=n次），每次递归都需要乘法。因此，它不能比O(M(n)•n)更好，而且这种指数会更糟。