HDU 4656 So Easy! - 矩陣快速冪解共軛遞推

題目大意

給定正整數 a, b, n, m，滿足 (a-1)² < b。
求 S_n = ⌈(a+√b)ⁿ⌉ mod m。

數學推導

共軛技巧

定義 A_n = (a+√b)ⁿ, B_n = (a-√b)ⁿ。
令 C_n = A_n + B_n，兩項共軛，C_n 為整數。

關鍵不等式

條件 (a-1)² < b ⇒ 0 < a-√b < 1 ⇒ 0 < B_n < 1。
因此 C_n = ⌈A_n⌉，所以 S_n = C_n mod m。

遞推關係

C_{n+1} = (a+√b)ⁿ⁺¹ + (a-√b)ⁿ⁺¹
= (a+√b)(a-√b)ⁿ + (a-√b)(a+√b)ⁿ... 整理得：
C_{n+1} = 2a·C_n - (a²-b)·C_{n-1}

初始值

• C_0 = 2
• C_1 = 2a

矩陣快速冪

[C_{n+1}, C_n]ᵀ = Mⁿ × [C_1, C_0]ᵀ
其中 M = [[2a, -(a²-b)], [1, 0]]

完整程式碼

typedef long long LL;
LL m;

struct Matrix {
    LL a[3][3];
    Matrix() { memset(a, 0, sizeof(a)); }
};

Matrix mul(Matrix A, Matrix B) {
    Matrix C;
    for (int i = 1; i <= 2; i++)
        for (int j = 1; j <= 2; j++) {
            C.a[i][j] = 0;
            for (int k = 1; k <= 2; k++)
                C.a[i][j] = (C.a[i][j] + A.a[i][k] * B.a[k][j]) % m;
            if (C.a[i][j] < 0) C.a[i][j] += m;
        }
    return C;
}

Matrix pow(Matrix A, LL p) {
    Matrix R; R.a[1][1] = R.a[2][2] = 1;
    while (p) {
        if (p & 1) R = mul(R, A);
        A = mul(A, A);
        p >>= 1;
    }
    return R;
}

int main() {
    LL a, b, n;
    while (scanf("%lld%lld%lld%lld", &a, &b, &n, &m) == 4) {
        Matrix M;
        M.a[1][1] = 2*a % m; M.a[1][2] = (b - a*a) % m;
        M.a[2][1] = 1; M.a[2][2] = 0;
        if (M.a[1][2] < 0) M.a[1][2] += m;
        
        M = pow(M, n);
        LL ans = (M.a[1][1] * (2*a % m) + M.a[1][2] * 2) % m;
        if (ans < 0) ans += m;
        printf("%lld\n", ans);
    }
    return 0;
}

複雜度

• 矩陣乘法：O(2³) = 常數
• 快速冪：O(log n)
• 極適合 n ≤ 10¹⁸ 等大數據

題目來源

HDU 4656 / Google Code Jam 2008 Round 1A Problem C。