您现在的位置：程式碼之家 >> 編程語言 > >> C語言 >> C++ >> C++ 入門

HDU 1588 Gauss Fibonacci：等差數列索引下的斐波那契加權求和與 4×4 矩陣快速冪

編輯：C++ 入門

問題數學建模

等差數列 g(i) = k·i + b。求 Σ_{i=0}^{n-1} f(g(i)) 模 M，其中 f 為標準斐波那契數列 (f(0)=0, f(1)=1)。k, b, n, M ≤ 10^9。

關鍵數學轉化

利用斐波那契矩陣形式：[f(n+1), f(n)]^T = A^n · [1, 0]^T，其中 A = [[1,1],[1,0]]。
f(k·i + b) = [A^{k·i+b}]_{1,0} = [A^b · (A^k)^i]_{1,0}。
令 B = A^k，則 Σ f(k·i+b) = [A^b · Σ B^i]_{1,0} = [A^b · (I - B^n)(I - B)^{-1}]_{1,0} (若可逆)。

4×4 增廣矩陣統一求和

構造塊矩陣遞推同時計算 B^i 與 Σ B^i：

C = [ B  0 ]
    [ I  I ]  (4×4 分塊矩陣，每塊 2×2)

則 C^n = [ B^n 0 ] [ ΣB^i I ]。提取右上 2×2 塊即為 Σ_{i=0}^{n-1} B^i。

完整程式碼 (C++)

#include <cstdio>
#include <cstring>
using namespace std;
typedef long long LL;
struct Mat2 { LL a[2][2]; };
struct Mat4 { LL a[4][4]; };

Mat2 mul2(const Mat2& A, const Mat2& B, LL mod) {
    Mat2 C{{{0}}};
    for(int i=0;i<2;++i) for(int k=0;k<2;++k) if(A.a[i][k])
        for(int j=0;j<2;++j) C.a[i][j]=(C.a[i][j]+A.a[i][k]*B.a[k][j])%mod;
    return C;
}
Mat4 mul4(const Mat4& A, const Mat4& B, LL mod) {
    Mat4 C{{{0}}};
    for(int i=0;i<4;++i) for(int k=0;k<4;++k) if(A.a[i][k])
        for(int j=0;j<4;++j) C.a[i][j]=(C.a[i][j]+A.a[i][k]*B.a[k][j])%mod;
    return C;
}
Mat2 qpow2(Mat2 b, LL e, LL mod) {
    Mat2 r{{{1,0},{0,1}}};
    while(e){ if(e&1) r=mul2(r,b,mod); b=mul2(b,b,mod); e>>=1; }
    return r;
}
Mat4 qpow4(Mat4 b, LL e, LL mod) {
    Mat4 r{{{0}}}; for(int i=0;i<4;++i) r.a[i][i]=1;
    while(e){ if(e&1) r=mul4(r,b,mod); b=mul4(b,b,mod); e>>=1; }
    return r;
}
int main() {
    LL k,b,n,M;
    while(~scanf("%lld%lld%lld%lld",&k,&b,&n,&M)) {
        Mat2 A{{{1,1},{1,0}}};
        Mat2 Ak = qpow2(A, k, M);
        Mat2 Ab = qpow2(A, b, M);
        Mat4 C{{{0}}};
        for(int i=0;i<2;++i) for(int j=0;j<2;++j) C.a[i][j] = Ak.a[i][j];
        for(int i=2;i<4;++i) for(int j=2;j<4;++j) C.a[i][j] = (i==j);
        for(int i=0;i<2;++i) for(int j=2;j<4;++j) C.a[i][j] = (i==j-2);
        Mat4 Cn = qpow4(C, n, M);
        Mat2 Sum{{{0}}};
        for(int i=0;i<2;++i) for(int j=0;j<2;++j) Sum.a[i][j] = Cn.a[i][j+2];
        Mat2 Res = mul2(Ab, Sum, M);
        printf("%lld\n", Res.a[1][0] % M);
    }
    return 0;
}