3233:Matrix Power Series

keyword

多項式　行列　分割統治　C++

概要

行列A(30*30以下)とk(<10^9)とm(<10^4)が与えられる。S=A+A^2+…+A^kをmod mで求める問題。
kが大きいのでO(log k)は必須。まさか公式一発ということもないだろう。そう考えると分割統治の発見は容易で、{A+…+A^(k/2)} + A^(k/2){A+…+A^(k/2)}で求まりそう。分割統治はなかなか書けないからちょっと嬉しかった。手元で最大ケースだと1.8秒くらいで、制限時間は3秒だから多分大丈夫だと思ってsubmitしたらTLE。行列を全部参照渡しに書き直して再送信したら無事ACがとれた。

#define matrix vector< vector<int> >

int n, m;
matrix a(30, vector<int>(30,0)), e(30, vector<int>(30,0));

inline matrix prod(matrix &x, matrix &y){
    matrix ret(30, vector<int>(30,0));
    int i, j, k;
    REP(i,n)REP(j,n)REP(k,n) ret[i][j] += x[i][k]*y[k][j]%m;
    REP(i,n)REP(j,n) ret[i][j] %= m;
    return ret;
}

inline matrix add(matrix &x, matrix &y){
    int i,j;
    vector<vector<int> > ret(30, vector<int>(30,0));
    REP(i,n)REP(j,n) ret[i][j] = (x[i][j] + y[i][j])%m;
    return ret;
}

inline matrix power(int k, matrix &x){
    if(k==1) return x;
    if(k&1){
        matrix t1 = power(k-1,x);
        return prod(x, t1);
    }
    matrix y = prod(x,x);
    return power(k/2, y);
}

inline matrix solve(int k){
    if(k==1) return a;
    matrix t1 = solve(k/2), t2 = power(k/2,a);
    t2 = add(e,t2);
    matrix ret = prod(t1, t2);
    if(k&1){
        t1 = power(k,a);
        ret = add(ret, t1);
    }
    return ret;
}

int main(){
    int k, i, j;
    scanf("%d%d%d",&n,&k,&m);
    REP(i,n)REP(j,n){
        if(i==j) e[i][j] = 1;
        scanf("%d",&a[i][j]);
    }
    matrix ans = solve(k);
    REP(i,n){
        REP(j,n){
            printf("%d",ans[i][j]%m);
            if(j==n-1)printf("\n");
            else printf(" ");
        }
    }
    return 0;
}