用矩阵计算斐波那契数列

我们都知道斐波那契数列存在递推关系 $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{matrix} a_{n} \\ a_{n - 1} \end{matrix}) = A \cdot (\begin{matrix} a_{n - 1} \\ a_{n - 2} \end{matrix})

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

(\begin{matrix} a_{n - 1} + a_{n - 2} \\ a_{n - 1} \end{matrix}) = A \cdot (\begin{matrix} a_{n - 1} \\ a_{n - 2} \end{matrix})

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

A = (\begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix})

矩阵形式的递推公式就是

(\begin{matrix} a_{n} \\ a_{n - 1} \end{matrix}) = (\begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix}) \cdot (\begin{matrix} a_{n - 1} \\ a_{n - 2} \end{matrix})

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

(\begin{matrix} a_{n} \\ a_{n - 1} \end{matrix}) = {(\begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix})}^{n - 1} \cdot (\begin{matrix} a_{1} \\ a_{0} \end{matrix})

其中 ${(\begin{matrix} a_{1} & a_{0} \end{matrix})}^{T}$ 应该是 ${(\begin{matrix} 1 & 0 \end{matrix})}^{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 (l o g_{2} n)$ 。假设我们想要计算 $3^{7}$ ，我们需要计算 $3 \times 3 \times 3 \times 3 \times 3 \times 3 \times 3$ ，需要计算 $7$ 次但是转换一下思路，我们可以根据结合律计算 $3^{6} \times 3$ ，而 $3^{6}$ 又可以写成 $3^{2} \times 3^{2} \times 3^{2}$ 需要计算4次。可以看到，当 $n$ 比较大的时候，这样做减少计算过程。推广一下，当计算 $a^{n}$ 时，得到递推公式

a^{n} = {\begin{cases} a^{n / 2} \times a^{n / 2} & n 是 偶 数 \\ a^{\frac{(n - 1)}{2}} \times a^{\frac{(n - 1)}{2}} \times a & n 是 奇 数 \\ 1 & 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} \dots 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]