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

Posted 2023-02-21

【中文标题】n 大到 10^19 的第 N 个斐波那契数？

【英文标题】：Nth Fibonacci number for n as big as 10^19?

【发布时间】：2015-04-17 09:05:39

【问题描述】：

我正在尝试编写一个程序来查找 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

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

【问题讨论】：

你可以使用循环而不是递归。

有什么不好的？

您将达到递归深度限制。 ***.com/questions/3323001/maximum-recursion-depth

我收到无效响应。

除了创建堆栈溢出的问题之外，您可能只想存储最后两个斐波那契数，这样您就不会创建一个由 10^19 大整数组成的数组。另外，可能，请查看像 gmpy2 这样的多精度整数库。

【参考方案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 中运行的东西。

【讨论】：

OP 没有很好地描述它，但我很确定 OP 的 % 1000000007 暗示了我们只需要得到答案 mod 1000000007 的事实。矩阵形式（或减少公式，如您所愿）可能无论如何都是必要的，因为您无法为上限进行〜10 ^ 19次迭代。

@DSM 你这样做的方式是首先不进行迭代。有一种更有效的方法来计算斐波那契数。

@will：我不确定你的意思，因为我刚刚说过迭代是不可能的。使用矩阵乘法或等效的约简公式（就像我刚才所做的那样——我看到 John Pirie 刚刚发布），我可以在大约 190 ns 内得到正确的答案。

@DSM 我只是在输入一个类似这样的答案：-/

@DSM 我没有正确阅读您写的内容。我同意你的看法。

【参考方案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)

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

【讨论】：

+1，虽然这个公式对于直接计算F(n) 到n 到10^19 仍然没有希望。 （这里没有公式可以工作：结果太大而无法存储。）但是结合减少模1000000007，这会起作用。

@Mark Dickinson：在 log(n) 复杂度下，我认为这个公式需要 50 次左右的迭代，不是吗？需要计算的附属值太多？

@JohnPirie：我认为他只是指 Fib(10^19) ~ 2.2041233236015342e+2089876402499787337 的事实，所以除非我们减少，否则我们会被冲洗掉。 :-)

@DSM：啊，所以一个简单的估计同样有效；谢谢

@JohnPirie：是的，帝斯曼说的。 OP 没有直接这么说，但看起来他真正想要的是减少 F(n) 模 1000000007 而不是 F(n) 本身。 （听起来像是典型的 Project-Euler 风格的挑战问题，而不是真实世界的计算。）

【参考方案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 个斐波那契数所需的一切，没有递归。我会在明天/本周晚些时候完成它...

【讨论】：

如果你真的在认真做这件事，那么你应该对矩阵进行对角化 - 然后你可以轻松地将它提升到任意幂。

嘿@will，这对斐波那契数列有很大帮助。但是，有点离题，但我希望你能提供帮助 - 我有一个整数序列，其中包含 2n 和 2n + 1 项的自定义公式。你知道我是否可以以类似于斐波那契序列的方式解决问题并为自定义序列制作类似的 Q 矩阵吗？谢谢！

什么是递归关系？如果偏移量是固定的（即它是constant recursive sequence），那么您总是可以构造这个矩阵（它只是大小不同）。如果它是相对的（即 4th 是 4/2 = 2nd 和 4/2+1 = 3rd 的函数，20th 是 10th 和 11th 的函数，等等），那么你不能 - 但仍然有更多方法可以获得解决方案轻松 - 发布问题。

仅供参考，对于任何阅读本文的人，你们中的任何人都沿着对角化路线走，然后您可以为第 n 个斐波那契数去掉一个解析的、非递归的公式。

【参考方案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 中找到。