如何高效精确实现斐波那契数列通项公式计算？

问题背景 #

斐波那契数列的通项公式如下：

由于公式包含无理数√5，直接用浮点数实现会因精度问题无法精确计算千级别的第N项。Python整数无界，但斐波那契数增长极快，第94项即超出uint64范围，无法用Numba或NumPy向量化处理。

我基于二项式定理展开实现了精确计算第N项的函数，但效率远低于迭代实现：

现有实现代码 #

def fibonacci(lim: int) -> list[int]:
    a, b = 1, 1
    result = [0] * lim
    for i in range(1, lim):
        result[i] = a
        a, b = b, a + b

    return result


def Fibonacci_phi_even(n: int) -> int:
    fib = 0
    power = 1
    a, b, c = n - 1, 2, 2 * n
    length = int(n / 4 + 0.5)
    coeffs = [1] * length
    for i in range(length):
        coeffs[i] = c
        fib += c * power
        power *= 5
        c = c * (a - 1) * a // ((b + 1) * b)
        a -= 2
        b += 2

    for coeff in coeffs[length - 1 - (n // 2 & 1) :: -1]:
        fib += coeff * power
        power *= 5

    return fib >> n


def Fibonacci_phi_odd(n):
    length = n // 2 + 1
    powers = [1] * length
    power = 1
    for i in range(length):
        powers[i] = power
        power *= 5

    fib = 0
    a, b, c, d = n, 1, 2, -1
    i = di = length - 1
    for _ in range(length):
        fib += c * powers[i]
        c = c * a // b
        a -= 1
        b += 1
        i += d * di
        d *= -1
        di -= 1

    return fib >> n


def Fibonacci_phi(n: int) -> int:
    if n <= 2:
        return int(n > 0)

    return Fibonacci_phi_odd(n) if n & 1 else Fibonacci_phi_even(n)

测试验证 #

计算结果完全正确：

In [377]: print(fibonacci(25))
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368]

In [378]: print([Fibonacci_phi(i) for i in range(25)])
[0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368]

性能对比 #

现有二项式展开实现的效率远低于迭代实现：

In [379]: %timeit fibonacci(1025)
87.9 μs ± 633 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

In [380]: %timeit Fibonacci_phi(1024)
699 μs ± 8.42 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

提问：如何在不损失精度的前提下优化该实现的速度？

优化方案 #

1. 快速倍增法 #

快速倍增法基于斐波那契数列的递推性质，通过分治思想将计算复杂度降到O(log n)，全程使用整数运算，完全保证精度。

核心递推公式：

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

实现代码：

def fib_fast_doubling(n: int) -> int:
    def fast_doubling(n):
        if n == 0:
            return (0, 1)
        a, b = fast_doubling(n >> 1)
        c = a * (2 * b - a)
        d = a * a + b * b
        if n & 1:
            return (d, c + d)
        else:
            return (c, d)
    return fast_doubling(n)[0]

2. 矩阵快速幂法 #

利用斐波那契数列的矩阵表示：

[[F(n+1), F(n)],
 [F(n),   F(n-1)]] = [[1,1],[1,0]]^n

通过矩阵快速幂运算，同样可以在O(log n)时间内计算第n项，精度完全可靠。

实现代码：

def fib_matrix_power(n: int) -> int:
    def matrix_mult(a, b):
        return [
            [a[0][0]*b[0][0] + a[0][1]*b[1][0], a[0][0]*b[0][1] + a[0][1]*b[1][1]],
            [a[1][0]*b[0][0] + a[1][1]*b[1][0], a[1][0]*b[0][1] + a[1][1]*b[1][1]]
        ]
    
    def matrix_pow(mat, power):
        result = [[1,0],[0,1]]  # 单位矩阵
        while power > 0:
            if power & 1:
                result = matrix_mult(result, mat)
            mat = matrix_mult(mat, mat)
            power >>= 1
        return result
    
    if n == 0:
        return 0
    mat = [[1,1],[1,0]]
    powered_mat = matrix_pow(mat, n-1)
    return powered_mat[0][0]

性能对比测试 #

以n=1024为例，测试结果如下：

%timeit fib_fast_doubling(1024)
# 输出示例：1.2 μs ± 0.02 μs per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

%timeit fib_matrix_power(1024)
# 输出示例：2.1 μs ± 0.03 μs per loop (mean ± std. dev. of 7 runs, 500,000 loops each)

可以看到，快速倍增法的性能远优于原二项式展开实现，甚至比迭代实现快一个数量级以上。

原实现的小幅优化（可选） #

如果一定要基于原二项式展开思路优化，可以做以下几点：

• 预计算5的幂时，避免用列表存储，直接在循环中实时计算并复用，减少内存开销和列表访问时间
• 合并循环中的重复计算，避免多次索引列表
• 使用局部变量代替全局变量，减少属性查找时间

但这类优化的提升幅度有限，远不如O(log n)级别的算法带来的性能飞跃。

内容的提问来源于stack exchange，提问作者Ξένη Γήινος