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

目录

• 矩阵算法
• 通项公式
• 测试

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

矩阵算法

斐波那契数的矩阵形式

$$\begin{array}{rl}\left[\begin{array}{c}{F}_n\\ F_{n-1}\end{array}\right]& =\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]\cdot \left[\begin{array}{c}{F}_{n-1}\\ F_{n-2}\end{array}\right]\\ & ={\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]}^{n-1}\cdot \left[\begin{array}{c}{F}_1\\ F_0\end{array}\right]\end{array}$$

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

$$\begin{array}{rl}\left[\begin{array}{c}{F}_n\\ F_{n-1}\end{array}\right]& \stackrel{m\leftarrow n-1}{=}\prod _{i=0}^{\lfloor \lg m\rfloor }{\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]}^{2^i{m}_i}\cdot \left[\begin{array}{c}{F}_1\\ F_0\end{array}\right]\end{array}$$

可以先把 ${\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]}^{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\times 2$ 的矩阵乘法；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{{\varphi }^i-{\hat{\varphi }}^i}{\sqrt{5}}$，其中 $\varphi =\frac{1+\sqrt{5}}{2},\hat{\varphi }=\frac{1-\sqrt{5}}{2}$，考虑到 $\frac{|\hat{\varphi }|}{\sqrt{5}}<\frac{1}{2}$，得到 $F_i=\lfloor \frac{{\varphi }^i}{\sqrt{5}}+\frac{1}{2}\rfloor$，即斐波那契数以指数形式增长。

$$\begin{array}{rl}m[i]& =\left[\begin{array}{cc}{F}_{2^i+1}& F_{2^i}\\ F_{2^i}& F_{2^i-1}\end{array}\right]\\ m[i]\cdot m[i]& =\left[\begin{array}{cc}{F}_{2^i+1}^2+F_{2^i}^2& F_{2^i+1}{F}_{2^i}+F_{2^i}{F}_{2^i-1}\\ F_{2^i+1}{F}_{2^i}+F_{2^i}{F}_{2^i-1}& F_{2^i}^2+F_{2^i-1}^2\end{array}\right]\\ T(m[i]\cdot m[i])& \approx 5c{n}^{\alpha },\text{ where }n=2\lfloor \lg F_{2^i}\rfloor =2\lg \lfloor \frac{{\varphi }^{2^i}}{\sqrt{5}}+\frac{1}{2}\rfloor \approx 2^{i+1}\lg \varphi -\lg 5\\ & \approx 5c(2^{i+1}\lg \varphi {)}^{\alpha }\\ T(n)& =\sum _{i=0}^{\lfloor \lg (n-1)\rfloor -1}T(m[i]\cdot m[i])\\ & \approx 5c{\lg }^{\alpha }\varphi \sum _{i=0}^{\lfloor \lg (n-1)\rfloor -1}{2}^{\alpha (i+1)}\\ & =c^{\prime }(2^{\alpha (\lfloor \lg (n-1)\rfloor +1)}-1)\\ & =\mathrm{\Theta }(2^{\alpha \lg n})=\mathrm{\Theta }(n^{\alpha })\end{array}$$

因此实际上的时间复杂度为 $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{array}{rl}\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]\cdot \left[\begin{array}{cc}a+b& a\\ a& b\end{array}\right]& =\left[\begin{array}{cc}2a+b& a+b\\ a+b& a\end{array}\right]\\ & \stackrel{a^{\prime }\leftarrow a+b,b^{\prime }\leftarrow a}{=}\left[\begin{array}{cc}{a}^{\prime }+b^{\prime }& a^{\prime }\\ a^{\prime }& b^{\prime }\end{array}\right]\end{array}$$

所以形如 $\left[\begin{array}{cc}a+b& a\\ a& b\end{array}\right]$ 的对称矩阵与 $\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right]$ 相乘仍然是对称矩阵。

$$\begin{array}{rl}\left[\begin{array}{cc}p+q& p\\ p& q\end{array}\right]\cdot \left[\begin{array}{cc}p+q& p\\ p& q\end{array}\right]& =\left[\begin{array}{cc}(p+q{)}^2+p^2& p^2+2pq\\ p^2+2pq& p^2+q^2\end{array}\right]\\ & \stackrel{p^{\prime }\leftarrow p^2+2pq,q^{\prime }\leftarrow p^2+q^2}{=}\left[\begin{array}{cc}{p}^{\prime }+q^{\prime }& p^{\prime }\\ p^{\prime }& q^{\prime }\end{array}\right]\\ \left[\begin{array}{cc}p+q& p\\ p& q\end{array}\right]\cdot \left[\begin{array}{cc}a+b& a\\ a& b\end{array}\right]& =\left[\begin{array}{cc}(p+q)(a+b)+ap& a(p+q)+bp\\ a(p+q)+bp& ap+bq\end{array}\right]\\ & \stackrel{a^{\prime }\leftarrow a(p+q)+bp,b^{\prime }\leftarrow ap+bq}{=}\left[\begin{array}{cc}{a}^{\prime }+b^{\prime }& a^{\prime }\\ a^{\prime }& b^{\prime }\end{array}\right]\end{array}$$

因此由于矩阵的特殊性，对于 m 我们只需要保存两个变量 $p,q$，对于 result 也只需要保存两个变量 $a,b$。

初始时有 $m=\left[\begin{array}{cc}1& 1\\ 1& 0\end{array}\right],result=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$，所以初始化 $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{{\varphi }^i-{\hat{\varphi }}^i}{\sqrt{5}}$，其中 $\varphi =\frac{1+\sqrt{5}}{2},\hat{\varphi }=\frac{1-\sqrt{5}}{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{{\varphi }^i-{\hat{\varphi }}^i}{\sqrt{5}}$ 看上去简洁明了，但是由于求指数没有经过优化、比较花时间。

sympy 内置的 fibonacci() 实际上调用了 mpmath.libmp 库中的 ifib() 函数。在计算的时候使用的数据类型不是 sympy.Integer，而是其他更接近底层的数据类型，所以速度会更快。

使用 np.log2() 取 $\lg$，再使用 np.polyfit(n, t, 1) 进行直线拟合，得到曲线和系数

得到 $\alpha \approx 1.58$，符合前面的时间复杂度的计算。

如果使用 np.int32 作为数据类型，那么就没有大数字的开销了，曲线为对数。