寫文動機

偶然在網路上看到 ([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: 程式碼實作