python - 快速斐波那契计算

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。所以第 N 个斐波那契数可以预期有大约 N/5 位数字。

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

以这种速度计算十亿分之一的斐波那契数会非常慢(即使它只有 ~ 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

将运行时间缩短到 0.008 秒和 0.31 秒(从使用 len(str(...)) 时的 0.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

计算 F(100,000) in .006, F(1,000,000) in .235 and F(10,000,000) in 9.51 seconds.

这是可以预料的。对于最快的测试，它产生结果的速度快了 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% 的加速。但是 lru_cache F(2N) 和 F(2N-1) solution given by Eli Korvigo below实际上给出了 加速4倍 (即，75%)。因此，虽然我还没有做出正式的解释，但我很想认为它在 N 的二进制扩展中缓存了 1 的跨度，并进行了最少数量的乘法运算。这样就不需要找到这些范围，预先计算它们，然后在 N 展开式中的正确点乘以它们。 lru_cache允许从上到下计算本来会更复杂的从下到上计算。

两者 SymTwoByTwoMatrix当 N 增长 10 倍时，lru_cache-of-F(2N)-and-F(2N-1) 的计算时间大约要长 40 倍。我认为这可能是由于 Python 实现了 long int 的乘法。我认为大数的乘法及其加法应该是可并行化的。因此，即使(正如 Daniel Fisher 在评论中所说)F(N) 解决方案是 Theta(n)，多线程 sub-O(N) 解决方案也应该是可能的。 .

Answer 1

由于斐波那契数列是线性递推，其成员可以以封闭形式计算。这涉及计算幂，这可以在 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

这是基于 note 的推导Edsger Dijkstra 教授。该解决方案利用了这样一个事实，即要计算 F(2N) 和 F(2N-1)，您只需要知道 F(N) 和 F(N-1)。尽管如此，您仍在处理长数算术，尽管开销应该小于基于矩阵的解决方案的开销。在 Python 中，由于记忆化和递归速度较慢，您最好以命令式风格重写它，尽管我这样编写是为了使函数公式清晰。