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