n 大到 10^19 的第 N 个斐波那契数？

Question

我正在尝试编写一个程序来查找 1

这是我使用动态编程的代码。

memo = {}
def fib(n):
    if n in memo:
        return memo[n]
    if n <= 2:
        f = 1
    else:
        f = fib(n-1) + fib(n-2)
    memo[n]=f
    return f
print fib(input()) % 1000000007

我的代码似乎不适用于大量数字。我收到无效响应错误。 有什么建议吗？

Answer 1

Python 的默认 recursion limit 为 1000（通常）。要了解您的系统的确切限制：

>>> import sys
>>> sys.getrecursionlimit()

首先，如果您想递归地编写此代码并且您使用的是 Python 3.2 及更高版本（从 print 语句来看，它看起来不像您），那么您可以像这样使用 @functools.lru_cache(maxsize=128, typed=False)：

import functools

@functools.lru_cache()
def fib(n):
    if n <= 2:
        return 1
    else:
        return fib(n-1) + fib(n-2)

话虽如此，但对于大量数字来说，这仍然不会很快。更好的方法是编写一个迭代解决方案，在任何给定时间，您需要“记忆”的只是最后 2 个数字。

您当然可以使用matrix form 以获得更好的性能。

最终，因为n 与10**19 一样大，如果不给你OverflowError，你将很难编写任何在Python 中运行的东西。

Answer 2

以 O(n) 的效率，你永远不会到达那里。与代码无关，但Dijkstra's note "In honor of Fibonacci" 描述了一种以 O(log(n)) 效率求 F(n) 的方法。

F(2n-1) = F(n-1)^2 + F(n)^2

F(2n) = (2*F(n-1)+F(n))*F(n)

你不仅可以这样做，而且还可以递归地这样做。

Answer 3

在 N 为 10^19 时获得第 N 个斐波那契数，如果你以天真的方式去做（至少我猜它不会起作用），那是行不通的。

有一个很多更好的方法来做到这一点。这种技术适用于很多这样的系列。它被称为Fibonacci Q Matrix。

在哪里

这样想：

您有一些矩阵将向量 A 转换为 B：

填写这些条目很容易。特殊之处在于它现在是一个矩阵运算符，所以如果我们想要第 1000 个斐波那契数，我们只需要进行矩阵乘法即可。

你可以用一个循环来做到这一点，但是你需要很长时间才能达到 10^19，并且做 10^19 矩阵乘法（即使它们很小）也需要也有一段时间。

相反，我们采用了另一条捷径。 x^N 可以重写为它们总和为 N 的幂的乘积，即

x**100 == x**90 * x**10

所以目标是在不进行大量计算的情况下在索引中获得大量数字：

x**2 和x*x 一样难——它们花费的时间相同。但是x*x*x*x 给出了与(x**2)**2 相同的答案，但需要额外的乘法。当您获得更高的权力时，收益会变得更多。因此，如果您将指数分解为 2 的幂（任何幂都有效，但这是最简单的情况），

X**100 == X**64 * X**32 * X**4

即

X**100 == (((((X**2)**2)**2)**2)**2)**2 + ...

所以你要做的就是计算出你想要达到的总功率的二的幂，然后取 Q 矩阵的二的幂的乘积。

这似乎对我有用：

fib_matrix = [[1,1],
              [1,0]]

def matrix_square(A, mod):
    return mat_mult(A,A,mod)


def mat_mult(A,B, mod):
  if mod is not None:
    return [[(A[0][0]*B[0][0] + A[0][1]*B[1][0])%mod, (A[0][0]*B[0][1] + A[0][1]*B[1][1])%mod],
            [(A[1][0]*B[0][0] + A[1][1]*B[1][0])%mod, (A[1][0]*B[0][1] + A[1][1]*B[1][1])%mod]]


def matrix_pow(M, power, mod):
    #Special definition for power=0:
    if power <= 0:
      return M

    powers =  list(reversed([True if i=="1" else False for i in bin(power)[2:]])) #Order is 1,2,4,8,16,...

    matrices = [None for _ in powers]
    matrices[0] = M

    for i in range(1,len(powers)):
        matrices[i] = matrix_square(matrices[i-1], mod)


    result = None

    for matrix, power in zip(matrices, powers):
        if power:
            if result is None:
                result = matrix
            else:
                result = mat_mult(result, matrix, mod)

    return result

print matrix_pow(fib_matrix, 10**19, 1000000007)[0][1]

然后，你可以更进一步 - 它只是一个 2x2 矩阵，所以我们可以对角化它，然后得到第 n 个斐波那契数的公式，就像 n 的函数一样 - 没有递归。像这样：

如上所述，我们计算将我们从一个步骤带到下一步的矩阵：

然后是从一组数字到下一组数字的关系：

我们可以在哪里链接这些矩阵乘法：

没有什么可以阻止我们回到第一个斐波那契数字：

现在游戏变成了“我们如何将矩阵提升到 n 次方”——这正是上面代码中所做的。但是有比我上面提出的解决方案更好的方法。我们可以将 Q 矩阵分解为特征值和向量，这样写：

其中 U 是一个酉矩阵，包含 Q 的特征值，而 Λ是对应特征值的矩阵。这些特征值和向量是：

然后你使用这种分解方式的标准优势之一，当你将它提升到一个幂时，相邻的 U 矩阵和它的逆组合得到酉矩阵，留下一个 U 并且它是逆的在末端，中间有一系列对角矩阵，将它们提升到幂是微不足道的：

所以现在我们有了用一个公式来写第 n 个斐波那契数所需的一切，没有递归。我会在明天/本周晚些时候完成它...

Answer 4

我不认为你可以达到 1E19，但这里是你如何避免双重溢出和递归深度限制：

import decimal
import operator


def decimal_range(start, stop, step=1):
    """Provides an alternative to `xrange` for very high numbers."""
    proceed = operator.lt
    while proceed(start, stop):
        yield start
        start += step


def fib(n):
    """
    Computes Fibonacci numbers using decimal.Decimal for high 
    precision and without recursion
    """
    a, b = decimal.Decimal(0), decimal.Decimal(1)
    for i in decimal_range(0, n):
        a, b = b, a + b
    return a

在我的机器上，计算1E6用了26.5秒，但我不能保证结果的正确性：

In [26]: %time f2(n)
CPU times: user 26.4 s, sys: 130 ms, total: 26.5 s
Wall time: 26.5 s
Out[26]: Decimal('1.953282128707757731632014830E+208987')

迭代器取自this SO thread，改动很小，而fib 函数可以在in this other thread 中找到。