问题

求斐波那契数列第$10^{19}$项的值,由于结果巨大,输出$mod$ $1000000007$的结果即可。

定义

前两项为1,从第三项开始都满足下面的公式
$$k_n = k_{n-1} + k_{n-2}$$ $$1,1,2,3,5,8,13\ldots$$

简单的python实现

def fib(n):
    a = 0
    b = 1
    for _ in range(n):
        a, b = b, a + b
    return a

可以看出,上面算法的时间复杂度是$O(n)$,可以说是比较优的方式了,但是当我在我的电脑上运行上面的代码求第$10^{19}$次项的值时,它很久都没有给我运行结果。

需要使用到斐波那契Q矩阵来减少运算次数

先介绍一下快速幂运算的概念

$x^2 = x × x$
$x^4 = x^2 × x^2$
$x^{100} = x^{64} × x^{32} × x^4$

观察上面的等式,$x^4$的计算如果用一次项乘4次需要计算3次,但是如果先计算2次项的平方则只需要计算2次。

同理,求$x^n$只需要计算$log2^n$次(实际操作是将低次项的结果保留起来,计算高次项时利用低次项的结果作为基础 )

同时, $mod$运算也满足 $(A×B)\%C=(A\%C)×(B\%C)$

使用幂运算的规则,所以只需要找到一个数或者表达式 $A$ 满足:

$$k_n = A × k_{n-1}$$

那么求第$n$项$fibonacci$值的问题就可以转化成求 $A^{n-1}$的问题了

但是通过$fibonacci$数列的定义可以知道,仅已知$k_{n-1}$无法确定$k_n$的值,必须要同时已知$k_{n-1}$和$k_{n-2}$才能确定$k_n$的值,即无法找到到 $A$ 满足$k_n = A × k_{n-1}$

重新定义序列

这里需要将原序列的相邻两项看作是一项,间隔为原序列的一项,组成新的序列满足

$$s_n = (k_{n-1}, k_n)$$

即:常量或常量表达式$A$ 需要满足:

$$s_n = A × s_{n-1}$$

也即为:

$$(k_{n-1}, k_n) = A × (k_{n-2}, k_{n-1})$$

又因为:

$$k_n = k_{n-1} + k_{n-2}$$

所以:

$$(k_{n-1}, k_{n-2}+k_{n-1}) = A × (k_{n-2}, k_{n-1})$$

可以联想到矩阵的乘法,于是将上面最后的结果用矩阵乘法表达:
$$ \begin{bmatrix} k_{n-2} & k_{n-1} \end{bmatrix} × A = \begin{bmatrix} k_{n-1} & k_{n-2} + k_{n-1} \end{bmatrix} $$

可以推导出:
$$A=\begin{bmatrix} 0&1\\ 1&1 \end{bmatrix}$$

这样的 $A$ 矩阵也被称为斐波那契Q矩阵,可以应用于解决很多问题。

下面是使用上述方式的python实现

def mat_sqrt(m, mod):
        return [
            [
                (m[0][0]*m[0][0] + m[0][1] * m[1][0]) % mod,
                (m[0][0]*m[0][1] + m[0][1] * m[1][1]) % mod,
            ],
            [
                (m[1][0]*m[0][0] + m[1][1] * m[1][0]) % mod,
                (m[1][0]*m[0][1] + m[1][1] * m[1][1]) % mod,
            ]
        ]

def mat_mul(m, n, mod):
    return [
        (m[0]*n[0][0] + m[1] * n[1][0]) % mod,
        (m[0]*n[0][1] + m[1] * n[1][1]) % mod
    ]

def fib(n):
    mod = 1000000007
    if n < 0:
        return
    elif n <= 2:
        return 1
    mat = [1, 1]

    q_matrix = [
        [0, 1],
        [1, 1]
    ]
    bin_n = bin(n - 2)[2:]
    matrixs = []
    
    for i in range(len(bin_n)):
        matrixs.insert(0, q_matrix)
        q_matrix = mat_sqrt(q_matrix, mod)

    for j, sbin in enumerate(bin_n):
        if sbin == '1':
            mat = mat_mul(mat, matrixs[j], mod)
    return mat[1]

注1:如果使用python,numpy已经提供了矩阵计算的方法,这里自己定义仅仅是为了演示,本文不提倡重复造轮子
注2:纯python int原则上可以存储无限大如果使用numpy 需要注意精度溢出问题