斐波那契数列的矩阵算法及 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 sympy.Integer(0)
n -= 1
bin_n = bin(n)[2:]
bool_n = list(map(bool, map(int, bin_n)))
bool_n.reverse()
m = [np.array([1, 1, 1, 0], dtype=dtype).reshape(2, 2)]
for _ in range(len(bool_n) - 1):
t = m[-1]
m.append(np.matmul(t, t))
needed = np.array(m)[bool_n]
result = reduce(np.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)\)。
但是考虑到使用大数字 sympy.Integer
进行乘法有与数字大小相关的开销。一般来说,这个开销约为 \(n^\alpha\),其中 \(\alpha\) 约为 \(1.5\) 到 \(2\),\(n\) 为相乘的两个数字的位数之和。
\(F_i = \frac{\phi^i-\hat\phi^i}{\sqrt5}\),其中 \(\phi = \frac{1+\sqrt5}{2}, \hat\phi = \frac{1-\sqrt5}{2}\),考虑到 \({\vert\hat\phi\vert\over\sqrt5}<\frac12\),得到 \(F_i=\left\lfloor{\phi^i\over\sqrt5}+\frac12\right\rfloor\),即斐波那契数以指数形式增长。
因此实际上的时间复杂度为 \(O(n^\alpha)\),这个 \(\alpha\) 和选择的大数字乘法算法有关。
将数字转为二进制和计算矩阵的过程融合。
def fib2(n, dtype=sympy.Integer):
m = np.array([1, 1, 1, 0], dtype=dtype).reshape(2, 2)
result = np.eye(2, dtype=dtype)
while n:
if n & 1:
result = result @ m
m = m @ m
n >>= 1
return result[1, 0]
也可以不采用矩阵的形式。这表示我们要保存上述代码中 m
和 result
的值。
考虑矩阵
所以形如 \(\begin{bmatrix}a+b&a\\a&b\end{bmatrix}\) 的对称矩阵与 \(\begin{bmatrix}1&1\\1&0\\\end{bmatrix}\) 相乘仍然是对称矩阵。
因此由于矩阵的特殊性,对于 m
我们只需要保存两个变量 \(p,q\),对于 result
也只需要保存两个变量 \(a,b\)。
初始时有 \(m=\begin{bmatrix}1&1\\1&0\\\end{bmatrix},result=\begin{bmatrix}1&0\\0&1\\\end{bmatrix}\),所以初始化 \(a=0,b=1,p=1,q=0\)。
def fib3(n, dtype=sympy.Integer):
one, zero = dtype(1), dtype(0)
a, b, p, q = zero, one, one, zero
while n:
if n & 1:
ap = a * p
a, b = ap + a*q + b*p, ap + b*q
pp = p * p
p, q = pp + 2*p*q, pp + q*q
n >>= 1
return a
通项公式
在我前面的文章算法导论 Introduction to Algorithms #算法基础中,也证明了 \(F_i = \frac{\phi^i-\hat\phi^i}{\sqrt5}\),其中 \(\phi = \frac{1+\sqrt5}{2}, \hat\phi = \frac{1-\sqrt5}{2}\)。
def fib4(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, 17))
f = [fib1, fib2, fib3, fib4, fibonacci]
t = [[] for _ in range(len(f))]
for i in n:
times = 2 ** i
for j in range(len(f)):
a = %timeit -o f[j](times)
t[j].append(a.average)
plt.figure(figsize=(12, 8))
for j in range(len(f)):
plt.subplot(2, 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()
函数。在计算的时候使用的数据类型不是 sympy.Integer
,而是其他更接近底层的数据类型,所以速度会更快。
使用 np.log2()
取 \(\lg\),再使用 np.polyfit(n, t, 1)
进行直线拟合,得到曲线和系数
得到 \(\alpha\approx1.58\),符合前面的时间复杂度的计算。
如果使用 np.int32
作为数据类型,那么就没有大数字的开销了,曲线为对数。