写文动机

偶然在网路上看到 ([1]) 讨论求费氏数列的做法，便去找了一些资料研究并实作。

矩阵快速幂

矩阵快速幂指的是当我们想求一个矩阵M的高次方时，可透过自乘过程中得到的值加速运算，举例来说当k=65，我们想要求 ，65次方可拆分成 1+64，而 64=32x2, 32=16x2, ...通过自乘我们可得到: ，因此我们只需要计算这些值即可。

参考程式码

def matrix_fast_power(mat_a: np.array, k) -> np.array:    if k == 1:        return mat_a    if k % 2 == 0:        h = matrix_fast_power(mat_a, k // 2)        return h @ h    else:        h = matrix_fast_power(mat_a, (k - 1) // 2)        return h @ h @ mat_a

广义费氏数列

形如 , 的数列为费氏数列 ([2])之推广。
而广义费氏数列一般项公式为: ，其中 为特徵方程式 之两根，详细推导可见此连结 ([3])。

推导

我们仿效此文 ([4])中的方法，令 ，并设 ，求解 相当于求解 。而 ，其中 ，可使用矩阵计算器 ([5])验证，搭配根与係数关係 。

参考程式码

直接对矩阵A计算:

def gen_fib_3(n: int) -> int:    # set up coefficients    a, b, p, q = Constant.a, Constant.b, Constant.p, Constant.q    if n == 0:        return p    if n == 1:        return q    mat_a = np.array([[a, b],                      [1, 0], ])    u_0 = np.array([[q],                    [p], ])    u_k = matrix_fast_power(mat_a, n - 1) @ u_0    return int(u_k[0])

先将矩阵A做对角化再计算:

def gen_fib_5(n: int) -> int:    # set up coefficients    a, b, p, q = Constant.a, Constant.b, Constant.p, Constant.q    if n == 0:        return p    if n == 1:        return q    u_0 = np.array([[q],                    [p], ])    roots = get_roots_of_c_polynomial(a, b)    alpha, beta = roots[0], roots[1]    mat_s = np.array([[alpha, beta],                      [1, 1], ])    alpha_power_k = num_fast_power(alpha, n - 1)    beta_power_k = num_fast_power(beta, n - 1)    mat_lambda_power_k = np.array([[alpha_power_k, 0],                                   [0, beta_power_k], ])    mat_inv_s = np.linalg.inv(mat_s)    u_k = mat_s @ mat_lambda_power_k @ mat_inv_s @ u_0    return int(np.round(u_k[0]))     def get_roots_of_c_polynomial(a: int, b: int) -> List[float]:    mat_a = np.array([[a, b],                      [1, 0], ])    roots = np.roots(np.poly(mat_a))    # sort in descending order    roots = -np.sort(-roots)    return roots

小结

通过矩阵快速幂求费氏数列第N项之时间複杂度为 O(logN)，而将矩阵对角化可大幅减少运算步骤，但需要判定是否可对角化以及求出特徵根。当推广至3项递迴 时，即不能保证总是可对角化，此时可直接对原矩阵A求快速幂。
我在程式实作中 import sympy,numpy 帮助求特徵根以及矩阵运算，过程中有转为浮点数进行运算，需留意求出的解是否有误差。

完整程式码 ([6])在此

参考资料

[1] geeksforgeeks: Program for Fibonacci numbers
[2] 维基百科：费波那契数
[3] 陈建烨: 一般的二阶线性递迴数列
[4] 周志成: 费布纳西数列的表达式
[5] 矩阵计算器
[6] hero0963: 程式码实作