快速斐波那契计算

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秒(从.03秒开始,len(str(...))使用时间为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%的加速.但是Eli Korvigo给出lru_cache的F(2N)和F(2N-1)解决方案实际上给出了4倍的速度(即75%).因此,虽然我没有得出正式的解释,但我很想到它在N的二进制扩展中缓存1的跨度并且需要最小的乘法次数.这样就不需要找到那些范围,预先计算它们然后在N的扩展中的正确点处将它们相乘. lru_cache允许从上到下计算本来更复杂的从顶部到顶部的计算.

SymTwoByTwoMatrix每当N增长10倍时,两者和lru_cache-of-F(2N)-and-F(2N-1)的计算时间大约要长40倍.我认为这可能是由于Python实现了长整数的乘法.我认为大数的乘法和它们的加法应该是可并行的.因此,即使(如Daniel Fisher在评论中所述)F(N)解决方案,也应该可以实现多线程子O(N)解决方案Theta(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教授的一份说明的推导.该解决方案利用了以下事实:要计算F(2N)和F(2N-1),您只需要知道F(N)和F(N-1).尽管如此,你仍然在处理长数量的算术,尽管开销应该小于基于矩阵的解决方案.在Python中,由于记忆和递归缓慢,你最好用命令式的方式重写它,尽管我这样写它是为了清晰的功能表达.

Answer 2

来自维基百科，

\n\n

对于所有 n \xe2\x89\xa5 0，数字 Fn 是最接近 phi^n/sqrt(5) 的整数，其中 phi 是黄金比例。因此，可以通过四舍五入的方式求得，即利用最接近的整数函数

\n