斐波那契数的矩阵算法及 python 实现

import numpy as np
import matplotlib.pyplot as plt
from functools import reduce
from sympy import sqrt, simplify, fibonacci
import sympy

斐波那契数的矩阵形式

对于一个整数 \(n\)，可以写为二进制 \(n=\sum_{i=0}^{\lfloor\lg n\rfloor}2^in_i\)，其中 \(n_i\) 为 \(n\) 的二进制的各位数，\(n_i\in\lbrace0,1\rbrace\)。于是把上述公式写为

可以先把 \(\begin{bmatrix}1&1\\1&0\\\end{bmatrix}^{2^i}\) 计算出来，再按照 \(m_i\) 相乘即可。

def fib1(n, dtype=sympy.Integer):
    if n == 0:
        return 0

    # n_i
    n -= 1  # 公式中的 n-1
    bin_n = bin(n)[2:]
    bool_n = list(map(bool, map(int, bin_n)))  # 转为 bool 数组，方便 numpy 索引
    bool_n.reverse()

    # [[1,1],[1,0]]^{2^i}
    m = [np.array([1, 1, 1, 0], dtype=dtype).reshape(2, 2)]
    for _ in range(len(bool_n) - 1):
        t = m[-1]
        m.append(t @ t)

    # produce
    needed = np.array(m)[bool_n]
    result = reduce(np.ndarray.__matmul__, needed, np.eye(2, dtype=dtype))  # 连续矩阵乘法
    return result[0, 0]

采用 sympy.Integer 作为数据类型可以防止溢出。

分析时空复杂度：m 的长度为 \(\lfloor\lg(n-1)\rfloor+1\)，需要进行 \(\lfloor\lg(n-1)\rfloor\) 次 \(2\times2\) 的矩阵乘法；needed 的长度由 n-1 的二进制中的 1 的个数决定，但也不会超过 m 的长度。故时空复杂度约为 \(O(\lfloor\lg(n-1)\rfloor+1)=O(\lg n)\)。

在我前面的文章算法导论 Introduction to Algorithms #算法基础中，也证明了 \(F_i = \frac{\phi^i-\hat\phi^i}{\sqrt5}\)$，其中 $$\phi = \frac{1+\sqrt5}{2}, \hat\phi = \frac{1-\sqrt5}{2}$

def fib2(n):
    p = (1 + sqrt(5)) / 2
    q = (1 - sqrt(5)) / 2
    return simplify((p ** n - q ** n) / sqrt(5))

尝试用 Jupyter Notebook 测试算法速度，其中 fibonacci() 为 sympy 自带的斐波那契函数，三者返回的结果都为 sympy.Integer。

n = list(range(8, 16))
f = [fib1, fib2, fibonacci]
t = [[] for _ in range(3)]

for i in n:
    times = 2 ** i
    for j in range(3):
        a = %timeit -o f[j](times)
        t[j].append(a.average)

plt.figure(figsize=(16, 3))
for j in range(3):
    plt.subplot(1, 3, j + 1)
    plt.plot(n, t[j], label=f[j].__name__)
    plt.legend()
plt.show()

虽然用公式 \(F_i = \frac{\phi^i-\hat\phi^i}{\sqrt5}\) 看上去简洁明了，但是由于求指数没有经过优化、比较花时间。

sympy 内置的 fibonacci() 实际上调用了 mpmath.libmp 库中的 ifib() 函数，是一个相当底层的库，所以很快。

原文地址：http://www.cnblogs.com/violeshnv/p/16928787.html