TZOJ 7929: Matrix Power Series

描述

Given a n × n matrix A and a positive integer k, find the sum S = A + A2 + A3 + … + Ak.

输入

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 109) and m (m < 104). Then follow n lines each containing n nonnegative integers below 32,768, giving A’s elements in row-major order.

输出

Output the elements of S modulo m in the same way as A is given.

样例输入

2 2 4
0 1
1 1

样例输出

1 2
2 3

解题思路：

首先 A^k通过矩阵快速幂求，以A^(k/2)为分界，

令     S(k)= A + A^2 + A^3 + … + A^k ，

则     S(n)=S(n-1)*A+A     把S(n)当成一个数，

然后利用矩阵快速幂求解，讨论k的奇数情况和偶数情况：

当k是偶数时

当k是奇数时

 

最后就是递推求解

思路理解快，代码写起来就很慢，看着样例数据出错误结果，真的心死了

AC代码

#include<bits/stdc++.h>
using namespace std;
#define IOS ios::sync_with_stdio(false);cin.tie(NULL);cout.tie(NULL)
typedef long long ll;
const int N=40,INF=0x3f3f3f3f;
int n,k,m;
ll ans[N][N],a[N][N];
void sp(ll x[][N]){
    for(int i=1;i<=n;i++)
        x[i][i]=1;
}
void march(ll x[][N],ll y[][N]){
    for(int i=1;i<=n;++i)
        for(int j=1;j<=n;++j)
            x[i][j]=y[i][j];
}
void mulity(ll x[][N],ll y[][N]){
    ll z[N][N]={0};
    for(int i=1;i<=n;++i){
        for(int j=1;j<=n;++j){
            for(int k=1;k<=n;k++)
                z[i][j]=(z[i][j]+x[i][k]*y[k][j]%m)%m;
        }
    }
    march(x,z);
}
void qiumi(ll x[][N],int kk){
    ll b[N][N]={0};
    int t=kk;
    march(b,a);
    while(t){
        if(t&1) mulity(x,b);
        mulity(b,b);
        t>>=1;//二进制求下一位 1 
    }
}
void add(ll x[][N],ll y[][N]){
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            x[i][j]=(x[i][j]+y[i][j])%m;
}
void solve(int kk,ll t[][N]){
    if(kk==0) return;
    ll f[N][N]={0},ff[N][N]={0};
    sp(f),sp(ff);
    qiumi(ff,kk/2);
    solve(kk/2,f);
    if(kk!=1){
        mulity(t,f);
        mulity(t,ff);
        add(t,f);
    }
    if(kk%2){//奇数 
        memset(ff,0,sizeof ff);
        sp(ff);
        qiumi(ff,kk);
        if(kk==1)
            march(t,ff);
        else
            add(t,ff);
    }
}
int main()
{    
    cin>>n>>k>>m;
    for(int i=1;i<=n;++i)
        for(int j=1;j<=n;++j)
            cin>>a[i][j];
    sp(ans);
    solve(k,ans);
    for(int i=1;i<=n;++i){
        for(int j=1;j<=n;++j){
            if(j>1) cout<<" ";
            cout<<ans[i][j];
        }
        if(i<n) cout<<'\n';
    }
    return 0;
}