用矩阵计算斐波那契数列

我们都知道斐波那契数列存在递推关系 $a_{n}=a_{n-1}+a_{n-2}$ 所以写出计算第N项的值，用递归很容易实现

def fibo(n):
    if n <= 2:
        return 1
    return fibo(n-1) + fibo(n-2)

但是要这么写，计算 $n$ 的值稍微大一些的时候就有可能爆栈了。如果不用递归，我们考虑迭代的方法，也非常直接

def fibo(n):
    cur = 1
    prev = 1
    if n <= 2:
        return cur
    for i in range(2, n):
        res = cur + prev
        prev = cur
        cur = res
    return cur

但是如果我们执意要用递归呢？如果Python支持的话，还有尾递归这个选项。

def fibo(n):
    def __fibo(n, prev, cur):
        if n <= 2:
            return cur
        return __fibo(n-1, cur, cur+prev)
    return __fibo(n, 1, 1)

尾递归的形式很像迭代的写法，只不过，将循环演化成了对自身的调用，cur和 prev这两个值放在了递归调用的参数上。除了以上三种常规写法，还有使用矩阵进行运算的，类似动态规划的算法。我们可以把递推关系推广到矩阵形式

\begin{pmatrix} a_{n} \\ a_{n-1} \end{pmatrix} = A \cdot \begin{pmatrix} a_{n-1} \\ a_{n-2} \end{pmatrix}

容易得到 $A$ 是一个 $2\times2$ 的矩阵,通过递推关系,可以得到

\begin{pmatrix} a_{n-1} + a_{n-2} \\ a_{n-1} \end{pmatrix} = A \cdot \begin{pmatrix} a_{n-1} \\ a_{n-2} \end{pmatrix}

那么 $A$ 的值就可以求得

A = \begin{pmatrix} 1 & 1\\ 1 & 0 \end{pmatrix}

矩阵形式的递推公式就是

\begin{pmatrix} a_{n} \\ a_{n-1} \end{pmatrix} = \begin{pmatrix} 1 & 1\\ 1 & 0 \end{pmatrix} \cdot \begin{pmatrix} a_{n-1} \\ a_{n-2} \end{pmatrix}

如果我们从头开始算起，那么递推公式可以表示成

\begin{pmatrix} a_{n} \\ a_{n-1} \end{pmatrix} = \begin{pmatrix} 1 & 1\\ 1 & 0 \end{pmatrix}^{n-1} \cdot \begin{pmatrix} a_{1} \\ a_{0} \end{pmatrix}

其中 $\begin{pmatrix}a_{1} & a_{0}\end{pmatrix}^T$ 应该是 $\begin{pmatrix}1 & 0\end{pmatrix}^T$ 使用 numpy来在 python中处理矩阵的乘法

import numpy as np
def fibo(n):
    matrix = np.array([[1, 1], [1, 0]])
    origin = np.array([[1], [0]])
    for i in range(n):
        origin = np.dot(matrix, origin)
    return origin[0][0]

矩阵的乘法，可以使用快速幂算法来减小时间复杂度快速幂可以把复杂度将为 $O(log_{2}n)$ 。假设我们想要计算 $3^{7}$ ，我们需要计算 $3\times3\times3\times3\times3\times3\times3$ ，需要计算 $7$ 次但是转换一下思路，我们可以根据结合律计算 $3^{6}\times3$ ，而 $3^{6}$ 又可以写成 $3^2\times3^2\times3^2$ 需要计算4次。可以看到，当 $n$ 比较大的时候，这样做减少计算过程。推广一下，当计算 $a^n$ 时，得到递推公式

a^n = \begin{cases}a^{n/2} \times a^{n/2} & \text n是偶数 \\ a^{\frac{(n-1)}{2}} \times a^{\frac{(n-1)}{2}} \times a & \text n是奇数 \\ 1 & \text n \end{cases}

因此我们可以根据它来写出Python递归函数计算 $a^n$

def fast_pow(a, n):
    if n <= 0:
        return 1
    if n % 2 == 0:
        return fast_pow(a, n/2) * fast_pow(a, n/2)
    return a * fast_pow(a, (n-1)/2) * fast_pow(a, (n-1)/2)

进一步考虑迭代的做法， $3^7=3^2 \times 3^4 \times 3^1$ ，从二进制的角度可以看到 $3^{(111)_{2}}=3^{2^{0}} \times 3^{2^{1}} \times 3^{2^{2}}$

也就是说，如果用二进制表示 $n = (t_1 t_2 t_3 \cdots t_n)_2$ 那么 $a^n=\prod_{i=0}^n a^{t_i \times 2^i}$

代码如下

def fast_pow(a, n):
    res = 1
    while n > 0:
        if n & 1 != 0:
            res = res * a
        a *= a
        n = n // 2
    return res

尾递归的代码也可以根据迭代的思路写出来

def fast_pow(a, n):
    def __fast_pow(res, a,  n):
        if n <= 0:
            return res
        if n & 1 == 0:
            return __fast_pow(res, a*a, n>>1)
        return __fast_pow(res*a, a*a, n>>1)
    return __fast_pow(1, a, n)

不过以上都是针对整数的快速幂的应用，那么该如何用快速幂来求解矩阵呢？并没有太大的区别，只不过将参与运算的整数和运算符号替换成矩阵的就可以了。

python代码如下

import numpy as np


def fast_pow(origin_matrix, n):
    def __fast_pow(res, matrix, n):
        if n <= 0:
            return res
        if n & 1 == 0:
            return __fast_pow(res, np.dot(matrix, matrix), n >> 1)
        return __fast_pow(np.dot(res, matrix), np.dot(matrix, matrix), n >> 1)
    unit_matrix = np.array([[1, 0], [0, 1]])
    return __fast_pow(unit_matrix, origin_matrix, n)


def fibo(n):
    matrix = np.array([[1, 1], [1, 0]])
    origin = np.array([[1], [0]])
    matrix = fast_pow(matrix, n-1)
    origin = np.dot(matrix, origin)
    return origin[0][0]

这里使用的尾递归的方式来计算快速幂，写的稍微好看一点

def fast_pow(times):
    def __fast_pow(init, el, n):
        def u_fast_pow(res, el, n):
            if n <= 0:
                return res
            if n & 1 == 0:
                return u_fast_pow(res, times(el, el), n >> 1)
            return u_fast_pow(times(res, el), times(el, el), n >> 1)

        return u_fast_pow(init, el, n)
    return __fast_pow




def fibo(n):
    matrix = np.array([[1, 1], [1, 0]])
    origin = np.array([[1], [0]])
    unit_matrix = np.array([[1, 0], [0, 1]])
    matrix = fast_pow(np.dot)(unit_matrix, matrix, n-1)
    origin = np.dot(matrix, origin)
    return origin[0][0]