import numpy as np
import matplotlib.pyplot as plt
from functools import reduce
from sympy import sqrt, simplify, fibonacci
import sympy
斐波那契数的矩阵形式
\[\begin{align*} \begin{bmatrix} F_n\\ F_{n-1}\\ \end{bmatrix}&= \begin{bmatrix} 1&1\\ 1&0\\ \end{bmatrix}\cdot \begin{bmatrix} F_{n-1}\\ F_{n-2}\\ \end{bmatrix}\\ &=\begin{bmatrix} 1&1\\ 1&0\\ \end{bmatrix}^{n-1}\cdot \begin{bmatrix} F_1\\ F_0\\ \end{bmatrix}\\ \end{align*} \]对于一个整数 \(n\)，可以写为二进制 \(n=\sum_{i=0}^{\lfloor\lg n\rfloor}2^in_i\)，其中 \(n_i\) 为 \(n\) 的二进制的各位数，\(n_i\in\lbrace0,1\rbrace\)。于是把上述公式写为
\[\begin{align*} \begin{bmatrix} F_n\\ F_{n-1}\\ \end{bmatrix}&\xlongequal{m\leftarrow n-1}\prod_{i=0}^{\lfloor\lg m\rfloor} \begin{bmatrix} 1&1\\ 1&0\\ \end{bmatrix}^{2^im_i}\cdot \begin{bmatrix} F_1\\ F_0\\ \end{bmatrix}\\ \end{align*} \]可以先把 \(\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() 函数，是一个相当底层的库，所以很快。