首页 文章

快速斐波那契计算

提问于
浏览 次
4

几周前我在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位内存） .

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

而且要清楚，我不是在寻找近似值 . 我正在寻找 exact computation （到最后一位数） .

Edit 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

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

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

python algorithm performance fibonacci

4 回答

  • 1

    由于Fibonacci序列是线性递归，因此可以以封闭形式评估其成员 . 这涉及计算功率，其可以与矩阵乘法解决方案类似地在O（logn）中完成，但是恒定开销应该更低 . 这是我所知道的最快的算法 .

    fib

    编辑

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

    fib2

    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中，由于记忆和递归的缓慢，你最好用命令式的方式重写它，尽管我这样做是为了清晰的功能表达 .

    回复于
  • 0

    在另一个答案中使用奇怪的平方根方程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
    回复于
  • 0

    从Wikipedia，

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

    回复于
  • 5

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

    Aaron的答案是正确的，我也赞成你，你应该这样做 . 我将提供相同的答案，并解释为什么它不仅正确，而且是迄今为止发布的最佳答案 . The formula我们正在讨论的是：

    formula

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

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

    然后有一个平方根（O(M(n))），一个分区（O(M(n))）和一个最后一轮（ O(n) ） .

    对于 n 位，这使得这个回答类似 O(n•log^2(n)•log(log(n))) .

    我没有't thoroughly analyzed the division algorithm, but if I'读这个权利，每个位可能需要一个递归（你需要将数字除以 log(2^n)=n 次），每个递归需要一个乘法 . 因此它不能比 O(M(n)•n) 更好，而且这种指数会更糟 .

    回复于

相关问题