基於矩陣快速冪求廣義費氏數列一般項

# 寫文動機偶然在[網路上看到] ([1]) 討論求費氏數列的做法，便去找了一些資料研究並實作。 # 矩陣快速冪 **矩陣快速冪**指的是當我們想求一個矩陣M的高次方時，可透過自乘過程中得到的值加速運算，舉例來說當 k=65，我們想要求 M^{65} ，可以分解成： M^{65} = M • M^{64} = M • M^{32} • M^{32} = ... 因此我們只需要計算矩陣的以下這些次方即可。 | 1 | 2 | 4 | 8 | 16 | 32 | 64 | |----|----|----|----|----|----|----| ### 參考程式碼 ``` 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 ``` # 廣義費氏數列形如 $$H_{n+2} = a H_{n+1} + b H_n、H_0, H_1 = p, q$$ 的數列為[費氏數列](https://zh.wikipedia.org/zh-tw/%E6%96%90%E6%B3%A2%E9%82%A3%E5%A5%91%E6%95%B0)的推廣。而 **廣義費氏數列** 的一般項公式為： $$ H_n = \frac{q - p \beta}{\alpha - \beta} \alpha^n - \frac{q - p \alpha}{\alpha - \beta} \beta^n $$ $$ \text{其中 } \alpha > \beta \text{ 為特徵方程式 } x^2 - ax - b = 0 \text{ 的兩根} $$ 詳細推導可見[此連結](https://ghresource.mt.ntnu.edu.tw/uploads/1644996440806kYlpaOsb.pdf)。 # 推導我們仿效[此文](https://ccjou.wordpress.com/2012/02/24/%E8%B2%BB%E5%B8%83%E7%B4%8D%E8%A5%BF%E6%95%B8%E5%88%97%E7%9A%84%E8%A1%A8%E9%81%94%E5%BC%8F/)中的方法，令 $$ u_k = \begin{bmatrix} H_{k+2} \\ H_{k+1} \end{bmatrix} \quad , A = \begin{bmatrix} a & b \\ 1 & 0 \end{bmatrix} \quad , u_0 = \begin{bmatrix} q \\ p \end{bmatrix} \quad $$ $$ \text{並設 } A = S \Lambda S^{-1} \text{，求解 } H_{k+2} \text{ 相當於求解 } u_k。 $$ $$ \text{而 } u_k = A u_{k-1} = A (A u_{k-2}) = \cdots = A^k u_0 = (S \Lambda S^{-1})^k u_0 = S \Lambda^k S^{-1} u_0 $$ 其中 $$ S = \begin{bmatrix} \alpha & \beta \\ 1 & 1 \end{bmatrix} \quad , \Lambda = \begin{bmatrix} \alpha & 0 \\ 0 & \beta \end{bmatrix} \quad $$ 可以使用[矩陣計算器](https://zs.symbolab.com/solver/matrix-multiply-calculator/)來驗證，搭配根與係數關係 $$ \alpha + \beta = a, \alpha \beta = -b $$ ### 參考程式碼直接對矩陣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(log N)**，而將矩陣對角化可以大幅減少運算步驟，但需要判定矩陣是否可對角化以及求出其特徵根。當推廣至 3 項遞迴 $$ H_{n+3} = a H_{n+2} + b H_{n+1} + c H_n $$ $$ H_0, H_1, H_2 = p, q, r $$ 時，即**不能保證總是可對角化**，此時可以直接對原矩陣 $ A $ 求快速冪。在程式實作中，我使用了 `sympy` 和 `numpy` 來幫助求特徵根及進行矩陣運算，過程中涉及浮點數運算，需要留意求出的解是否存在誤差。 [完整程式碼] ([6])在此 [網路上看到]:<https://www.geeksforgeeks.org/program-for-nth-fibonacci-number/> [geeksforgeeks: Program for Fibonacci numbers]:<https://www.geeksforgeeks.org/program-for-nth-fibonacci-number/> [費氏數列]: <https://zh.wikipedia.org/zh-tw/%E6%96%90%E6%B3%A2%E9%82%A3%E5%A5%91%E6%95%B0> [維基百科：費波那契數]: <https://zh.wikipedia.org/zh-tw/%E6%96%90%E6%B3%A2%E9%82%A3%E5%A5%91%E6%95%B0> [此連結]: <https://ghresource.mt.ntnu.edu.tw/uploads/1644996440806kYlpaOsb.pdf> [陳建燁: 一般的二階線性遞迴數列]: <https://ghresource.mt.ntnu.edu.tw/uploads/1644996440806kYlpaOsb.pdf> [此文]: <https://ccjou.wordpress.com/2012/02/24/%E8%B2%BB%E5%B8%83%E7%B4%8D%E8%A5%BF%E6%95%B8%E5%88%97%E7%9A%84%E8%A1%A8%E9%81%94%E5%BC%8F/> [周志成: 費布納西數列的表達式]: <https://ccjou.wordpress.com/2012/02/24/%E8%B2%BB%E5%B8%83%E7%B4%8D%E8%A5%BF%E6%95%B8%E5%88%97%E7%9A%84%E8%A1%A8%E9%81%94%E5%BC%8F/> [矩陣計算器]: <https://zs.symbolab.com/solver/matrix-multiply-calculator/%5Cbegin%7Bpmatrix%7D%5Calpha%20%26%5Cbeta%5C%5C%20%201%20%261%5Cend%7Bpmatrix%7D%5E%7B%20%7D%5Cbegin%7Bpmatrix%7D%5Calpha%20%26%200%5C%5C%200%20%26%5Cbeta%20%5Cend%7Bpmatrix%7D%5Cbegin%7Bpmatrix%7D%5Calpha%20%20%26%20%5Cbeta%20%5C%5C%201%20%261%5Cend%7Bpmatrix%7D%5E%7B-1%7D?or=input> [完整程式碼]: <https://github.com/Hero0963/ithelp/tree/main/cal_gen_fib_num> [hero0963: 程式碼實作]: <https://github.com/Hero0963/ithelp/tree/main/cal_gen_fib_num> # 參考資料 [1] [geeksforgeeks: Program for Fibonacci numbers] [2] [維基百科：費波那契數] [3] [陳建燁: 一般的二階線性遞迴數列] [4] [周志成: 費布納西數列的表達式] [5] [矩陣計算器] [6] [hero0963: 程式碼實作]