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

#python #費氏數列 #矩陣快速冪 #對角化

ahero0963 2022-08-10 15:26:39 ‧ 1490 瀏覽

分享至

寫文動機

偶然在網路上看到討論求費氏數列的做法，便去找了一些資料研究並實作。

寫文動機

偶然在網路上看到 ([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$$ 的數列為費氏數列的推廣。

而 廣義費氏數列 的一般項公式為：
$$
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{ 的兩根} $$

詳細推導可見此連結。

推導

我們仿效此文中的方法，令
$$
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
$$

可以使用矩陣計算器來驗證，搭配根與係數關係

$$
\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