POJ3233解题报告

Matrix Power Series

Time Limit: 3000MS Memory Limit: 131072KTotal Submissions: 16466 Accepted: 7012

Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A² + A³ + … + A^k.

Input

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

Output

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

Sample Input

2 2 40 11 1

Sample Output

1 22 3

Source

POJ Monthly--2007.06.03, Huang, Jinsong

           这道题是一道比较经典的矩阵快速幂运算的题目。在进行矩阵快速幂的过程中我们需要明确程序运行时间主要消耗在什么地方。其实在快速幂运算中最浪费时间的是取模运算。即在进行矩阵乘法时的取模运算很费时间。通过测试取模运算的处理会使得程序运行时间有10倍的差异。这一点是必须得注意的。为了防止超时，所以我们在做矩阵快速幂的题目时一定要注意在取模运算上多一些优化，防止超时。

       下面说这道题的解法，这道题是矩阵的等比求和，当然我们可以类比等比数列求和，直接利用公式求解。但是对于矩阵来说，这样会涉及到逆矩阵的计算，还有矩阵是否可逆的判断，比较麻烦，所以我们避开这种做法。比较好的方法有两种，一种是多项式处理的常见方法，即二分法对多项式的和进行处理。构造出来的递推式是n与n/2之间的递推关系式。这也说明了矩阵具有优化递推式的作用。使用矩阵处理递推关系也是很重要的方法。

      方法一：

      

           我们可以把一个矩阵看成一个元素，利用二分法+矩阵快速幂求出该矩阵的等比之和。

       参考代码：

       

#include<cstdio>#include<iostream>#include<cmath>#include<cstring>#include<algorithm>#include<string>#include<vector>#include<map>#include<set>#include<stack>#include<queue>#include<ctime>#include<cstdlib>#include<iomanip>#include<utility>#define pb push_back#define mp make_pair#define CLR(x) memset(x,0,sizeof(x))#define _CLR(x) memset(x,-1,sizeof(x))#define REP(i,n) for(int i=0;i<n;i++)#define Debug(x) cout<<#x<<"="<<x<<" "<<endl#define REP(i,l,r) for(int i=l;i<=r;i++)#define rep(i,l,r) for(int i=l;i<r;i++)#define RREP(i,l,r) for(int i=l;i>=r;i--)#define rrep(i,l,r) for(int i=1;i>r;i--)#define read(x) scanf("%d",&x)#define put(x) printf("%d\n",x)#define ll long long#define lson l,m,rt<<1#define rson m+1,r,rt<<11using namespace std;int n,k,m;struct mat{    int d[30][30];} A,E,M;mat multi(mat &a,mat &b){    mat ans;    rep(i,0,n)    {        rep(j,0,n)        {            ans.d[i][j]=0;            rep(k,0,n)                ans.d[i][j]+=a.d[i][k]*b.d[k][j];            ans.d[i][j]%=m;  //在循环外边取模这一步节省了10倍的时间，因为n最多是30,m又小于10^4,所以乘积不会超过10^8,所以乘积之和不会超过3*10^9,由于数据比较弱，所以这一层for循环求出的和可以保证在int范围内        }    }    return ans;}mat add(mat &a,mat &b){    mat ans;    rep(i,0,n)    {        rep(j,0,n)        {            ans.d[i][j]=a.d[i][j]+b.d[i][j];            ans.d[i][j]%=m;        }    }    return ans;}mat quickmulti(int n){    mat p=E,q=A;    if(n==0) return p;    if(n==1) return A;    while(n!=1)    {        if(n&1)        {            n--;            p=multi(p,q);        }        else        {            n>>=1;            q=multi(q,q);        }    }    return multi(p,q);}mat binarymulti(int k){    if(k==1) return A;    mat p=binarymulti(k/2),tmp,q;    if(k&1)          //k为奇数时sum(k)=(1+A^(k/2+1))*sum(k/2)+A^(k/2+1);    {        tmp=quickmulti(k/2+1);        q=multi(tmp,p);        p=add(p,q);        p=add(p,tmp);    }    else            //k为偶数时sum(k)=(1+A^(k/2))*sum(k/2)    {        tmp=quickmulti(k/2);        q=multi(tmp,p);        p=add(p,q);    }    return p;}int main(){    scanf("%d%d%d",&n,&k,&m);    rep(i,0,n)        rep(j,0,n)        {            read(A.d[i][j]);            A.d[i][j]%=m;   //这是比较关键的一步。题目中m的范围很小，所以提前取模，在矩阵乘法时取模次数就可以减少，防止溢出            E.d[i][j]=(i==j);        }    M=binarymulti(k);    rep(i,0,n)    {        rep(j,0,n)            printf("%d ",M.d[i][j]);        printf("\n");    }}

          方法二：

          利用好分块矩阵的性质，构造矩阵。只用一次矩阵快速幂。省掉了方法一中的二分多项式的步骤，使得运算速度有了很大提高。

         

          这个方法中需要用到分块矩阵，这样会有些麻烦。矩阵中套矩阵。不过速度比第一种方法要快一些。

          参考代码：

          

#include<cstdio>#include<iostream>#include<cmath>#include<cstring>#include<algorithm>#include<string>#include<vector>#include<map>#include<set>#include<stack>#include<queue>#include<ctime>#include<cstdlib>#include<iomanip>#include<utility>#define pb push_back#define mp make_pair#define CLR(x) memset(x,0,sizeof(x))#define _CLR(x) memset(x,-1,sizeof(x))#define REP(i,n) for(int i=0;i<n;i++)#define Debug(x) cout<<#x<<"="<<x<<" "<<endl#define REP(i,l,r) for(int i=l;i<=r;i++)#define rep(i,l,r) for(int i=l;i<r;i++)#define RREP(i,l,r) for(int i=l;i>=r;i--)#define rrep(i,l,r) for(int i=1;i>r;i--)#define read(x) scanf("%d",&x)#define put(x) printf("%d\n",x)#define ll long long#define lson l,m,rt<<1#define rson m+1,r,rt<<11using namespace std;int n,k,M;struct mat{    int d[30][30];};struct Mat{    mat m[2][2];} A,E;mat add(mat &a,mat &b){    mat ans;    rep(i,0,n)    {        rep(j,0,n)            ans.d[i][j]=(a.d[i][j]+b.d[i][j])%M;    }    return ans;}mat multi(mat &a,mat &b){    mat ans;    rep(i,0,n)    {        rep(j,0,n)        {            ans.d[i][j]=0;            rep(k,0,n)            {                if(a.d[i][k]&&b.d[k][j])                    ans.d[i][j]+=a.d[i][k]*b.d[k][j];            }            ans.d[i][j]%=M;        }    }    return ans;}Mat Multi(Mat &a,Mat &b){    Mat p;    CLR(p.m);    rep(i,0,2)    {        rep(j,0,2)        {            rep(k,0,2)            {                mat q=multi(a.m[i][k],b.m[k][j]);                p.m[i][j]=add(p.m[i][j],q);            }        }    }    return p;}Mat quickmulti(int n){    Mat p=E,q=A;    if(n==0) return p;    if(n==1) return A;    while(n)    {        if(n&1)        {            n--;            p=Multi(p,q);        }        else        {            n>>=1;            q=Multi(q,q);        }    }    return p;}int main(){    CLR(E.m);    rep(i,0,2)        rep(j,0,30)    E.m[i][i].d[j][j]=1;    scanf("%d%d%d",&n,&k,&M);    rep(i,0,30)        A.m[0][0].d[i][i]=1;    rep(i,0,n)    {        rep(j,0,n)        {            read(A.m[0][1].d[i][j]);            A.m[0][1].d[i][j]%=M;            A.m[1][1].d[i][j]=A.m[0][1].d[i][j];        }    }    Mat m1=quickmulti(k);    rep(i,0,n)    {        rep(j,0,n)            printf("%d ",m1.m[0][1].d[i][j]);        printf("\n");    }}