poj 3233 --- Matrix Power Series (二分,矩阵)

Matrix Power Series

Time Limit: 3000MS Memory Limit: 131072KTotal Submissions: 9393 Accepted: 4018

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

S_k = A + A² + A³ + … + A^k  

    =(1+A^k/2)*(A + A² + A³ + … + A^k/2  )+{A^k}

    =(1+A^k/2)*(S_k/2 )+{A^k}// k为偶数时无 {A^k}

A^k   可用二分迭代求出

因此,只要求出 上面的三部分就可以求出 S_k 

//800MS#include <iostream>#include <cstdio>#include <cstring>using namespace std;int m,n,K;int a[30][30];class Matrix{public:    int num[30][30];    Matrix(bool is=true)//初始化    {        memset(num,0,sizeof(num));if(is)for(int i=0;i<n;i++)num[i][i]=1;    }    void print()//输出函数    {        for(int i=0;i<n;++i)        {            printf("%d",num[i][0]);            for(int j=1;j<n;++j)                printf(" %d",num[i][j]);            printf("\n");        }    }//重载乘法运算   friend Matrix& operator *(const Matrix& max1,const Matrix& max2){Matrix tmp(false);//注意这里是false,即初始化的矩阵不是单位矩阵Ifor(int i=0;i<n;++i)for(int j=0;j<n;++j){for(int k=0;k<n;++k)tmp.num[i][j]+=(max1.num[i][k]*max2.num[k][j])%m;            tmp.num[i][j]%=m;}return tmp;}   //重载+=运算   Matrix& operator +=(const Matrix& max)   {   for(int i=0;i<n;++i)   for(int j=0;j<n;++j)   num[i][j]=(num[i][j]+max.num[i][j])%m;return *this;}}ans;Matrix mul(Matrix A,int k)//求A^K{if(k==1)return A;Matrix tmp ;while(k){if(k&1)tmp = tmp * A;k>>=1;A = A*A;}return tmp;}Matrix S(Matrix A,int k)//求 S[k]{if(k==1)return A;Matrix tmp ;tmp += mul(A,k>>1);//求 (I + A^(k/2) )tmp = tmp*S(A,k>>1);//求 (I + A^(k/2) )*S[k/2]if(k&1)//判断时候要加上 A^ktmp+= mul(A,k);//S[k] = (I+A^(k/2)) * S[k/2] + {A^k}return tmp;} int main(){int i,j,k;    scanf("%d %d %d",&n,&K,&m);    for( i=0;i<n;++i)        for( j=0;j<n;++j)            scanf("%d",&ans.num[i][j]);ans = S(ans,K);ans.print();return 0;}

效率更高的代码1:

//266MS#include <iostream>#include <cstdio>#include <cstring>using namespace std;int m,n,K;int a[30][30];class Matrix{    public:    int num[61][61];    Matrix()    {        memset(num,0,sizeof(num));    }    void print()    {        for(int i=0;i<n;++i)        {            printf("%d",num[i][0]);            for(int j=1;j<n;++j)                printf(" %d",num[i][j]);            printf("\n");        }        printf("\n");    }    friend Matrix& operator *(Matrix max1,Matrix max2);};Matrix& operator *(Matrix max1,Matrix max2){    Matrix tmp;    for(int i=0;i<2*n;++i)        for(int j=0;j<2*n;++j)        {            for(int k=0;k<2*n;++k)                tmp.num[i][j]+=(max1.num[i][k]*max2.num[k][j])%m;            tmp.num[i][j]%=m;        }    return tmp;}int main(){    scanf("%d %d %d",&n,&K,&m);    for(int i=0;i<n;++i)        for(int j=0;j<n;++j)            scanf("%d",&a[i][j]);    Matrix b,c;    for(int i=0;i<n;++i)        for(int j=0;j<n;++j)            c.num[i+n][j+n]=b.num[i][j+n]=a[i][j];    for(int k=0;k<2;++k)        for(int i=0;i<n;++i)            c.num[i+k*n][i]=1;    while(K)    {        if(K&1)            b=b*c;        K>>=1;        c=c*c;    }    b.print();    return 0;}

效率更高的代码2:

#include<cstdio>#include<cstring>using namespace std;int i,j,k,n,l,m,r;bool s[32];struct matrix{    int a[30][30];    void pr()    {        for(i=0; i<n;printf("%d/n",a[i++][j]))            for(j=0;j<n-1;++j)printf("%d ",a[i][j]);    }    matrix& operator+=(const matrix&x)    {        for(int i=0; i<n; ++i)            for(int j=0; j<n; ++j)            if(x.a[i][j])            {                a[i][j]+=x.a[i][j];                if(a[i][j]>=m)a[i][j]%=m;            }        return *this;    }    matrix& operator*=(const matrix&x)    {        matrix t;        int i,j,k;        for(i=0; i<n; ++i)            for(j=0; j<n; ++j)t.a[i][j]=0;        for(i=0; i<n; ++i)            for(k=0; k<n; ++k)                if(a[i][k])                    for(j=0; j<n; ++j)                    {                        t.a[i][j]+=a[i][k]*x.a[k][j];                        if(t.a[i][j]>=m)t.a[i][j]%=m;                    }        return memcpy(a,t.a,sizeof(a)),*this;    }} t,a,unit,ans,exp;int& fun(int &a){    return a;}int main(){    for(scanf("%d%d%d",&n,&l,&m),i=0; i<30; ++i)unit.a[i][i]=1;    for(i=0; i<n; ++i)        for(j=0; j<n; ++j)            scanf("%d",a.a[i]+j);    for(r=-1;l;l>>=1)s[++r]=l&1;    for(ans=exp=a; --r>=0; exp*=exp,s[r]?exp*=a:exp)        t=exp,s[r]?ans*=(t*=a)+=unit,ans+=(t=exp)*=a:ans*=t+=unit;    ans.pr();    return 0;}