[ACM] POJ 3233 Matrix Power Series （求矩阵A+A^2+A^3...+A^k，二分求和或者矩阵转化）

Matrix Power Series

Time Limit: 3000MS Memory Limit: 131072KTotal Submissions: 15417 Accepted: 6602

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

解题思路：

第一种方法：

题意为给定矩阵A（下面代码中用ori表示)，以及k, mod ,求 A+A^2+A^3+......A^k    的和对mod取余。

一开始用循环k次，递推的做法，超时。。。

看了解题报告，求和的时候要用到二分求和。

所求的和用s(k)表示。

当k为偶数时：

比如 k=6,那么  A+A^2+A^3+A^4+A^5+A^6=    A+A^2+A^3+   A^3*（A+A^2+A^3）   

s(k)=s(k/2)+A^(n/2) * s(k/2) 即s(k)=(E+A^(n/2))*s(n/2)  （E为单位矩阵）

当k为奇数时：

s(k)=s(k-1)+A^k ,  那么k-1为偶数，可以按照上面的二分

PS：代码要写仔细啊，否则一个小错误查半天.....计算两个矩阵相乘时ret.arr[i][j]+=a.arr[i][k]*b.arr[k][j]; 竟然写成了ret.arr[i][j]+=a.arr[i][k]*a.arr[k][j]; T T

代码：

#include <iostream>#include <stdio.h>#include <string.h>using namespace std;const int maxn=31;int n,k,mod;struct mat{    int arr[maxn][maxn];    mat()    {        memset(arr,0,sizeof(arr));    }};mat mul(mat a,mat b){    mat ret;    for(int i=0;i<n;i++)        for(int k=0;k<n;k++)    {        if(a.arr[i][k])            for(int j=0;j<n;j++)        {            ret.arr[i][j]+=a.arr[i][k]*b.arr[k][j];            if(ret.arr[i][j]>=mod)                ret.arr[i][j]%=mod;        }    }    return ret;}mat add(mat a,mat b){    mat an;    for(int i=0;i<n;i++)        for(int j=0;j<n;j++)        {            an.arr[i][j]=a.arr[i][j]+b.arr[i][j];            if(an.arr[i][j]>=mod)                an.arr[i][j]%=mod;        }    return an;}mat power(mat p,int k){    if(k==1) return p;    mat e;    for(int i=0;i<n;i++)        e.arr[i][i]=1;    if(k==0) return e;    while(k)    {        if(k&1)            e=mul(p,e);        p=mul(p,p);        k>>=1;    }    return e;}void output(mat ans){    for(int i=0;i<n;i++)        for(int j=0;j<n;j++)    {        if(j==n-1)            cout<<ans.arr[i][j]<<endl;        else            cout<<ans.arr[i][j]<<" ";    }}mat  cal(mat ori,int k){    if(k==1) return ori;    if(k&1)        return add(cal(ori,k-1),power(ori,k));//当k为奇数时，减1变为偶数 S(K)=S(K-1)+ori^K    else        return mul(add(power(ori,0),power(ori,k>>1)),cal(ori,k>>1));        //当K为偶数时,S(K)=(1+ori^(K/2))*S(K/2)}int main(){    while(scanf("%d%d%d",&n,&k,&mod)!=EOF)    {        mat ori,ans;        for(int i=0;i<n;i++)            for(int j=0;j<n;j++)        {            scanf("%d",&ori.arr[i][j]);            if(ori.arr[i][j]>=mod)                ori.arr[i][j]%=mod;        }        ans=cal(ori,k);        output(ans);    }    return 0;}

第二种方法：

解题思路转载自：http://www.cnblogs.com/jiangjing/archive/2013/05/28/3103336.html

分析：求a^1+..a^n这是矩阵乘法中关于等比矩阵的求法：

|A  E|

|0  E|

其中的A为m阶矩阵，E是单位矩阵，0是零矩阵。而我们要求的是：                                                                             

A^1+A^2+..A^L，由等比矩阵的性质

|A  ,  1|                 |A^n , 1+A^1+A^2+....+A^(n-1)|

|0  ,  1| 的n次方等于|0     ,                   1                   |

所以我们只需要将A矩阵扩大四倍，变成如上形式的矩阵B，然后开L+1次方就可以得到1+A^1+A^2+....+A^L。由于多了一个1，所以最后得到的答案我们还要减去1。同理我们把矩阵A变成B：

                                                          |A  E|

                                                          |0  E|

然后我们就是求B的n+1次幂之后得到的矩阵为|x1   x2|

                                                          |x3   x4|

右上角的矩阵x2减去单位矩阵E，得到就是要求的矩阵了！

#include <iostream>#include <string.h>using namespace std;const int maxn=62;int n,k,mod;struct mat{    int arr[maxn][maxn];    mat()    {        memset(arr,0,sizeof(arr));    }};mat mul(mat a,mat b){    mat ans;    for(int i=0;i<n;i++)        for(int k=0;k<n;k++)    {        if(a.arr[i][k])            for(int j=0;j<n;j++)        {            ans.arr[i][j]+=a.arr[i][k]*b.arr[k][j];            if(ans.arr[i][j]>=mod)                ans.arr[i][j]%=mod;        }    }    return ans;}mat power(mat p,int k){    if(k==1) return p;    mat e;    for(int i=0;i<n;i++)        e.arr[i][i]=1;    if(k==0) return e;    while(k)    {        if(k&1)            e=mul(e,p);        p=mul(p,p);        k>>=1;    }    return e;}void output(mat ans){    for(int i=0;i<n;i++)        for(int j=0;j<n;j++)    {        if(i==j)        {            ans.arr[i][j+n]--;            while(ans.arr[i][j+n]<0)                ans.arr[i][j+n]+=mod;        }        if(j==n-1)            cout<<ans.arr[i][j+n]<<endl;        else            cout<<ans.arr[i][j+n]<<" ";    }}int main(){    while(cin>>n>>k>>mod)    {        mat ori,ans;        for(int i=0;i<n;i++)            for(int j=0;j<n;j++)        {            cin>>ori.arr[i][j];            if(i==j)            {                ori.arr[i][j+n]=1;                ori.arr[i+n][j+n]=1;            }        }        n*=2;        ans=power(ori,k+1);        n/=2;        output(ans);    }    return 0;}