Matrix Power Series
Time Limit: 3000MS | Memory Limit: 131072K | |
Total Submissions: 24316 | Accepted: 10110 |
Description
Given a n × n matrix A and a positive integer k, find the sum S = A + A2 + A3 + … + Ak.
Input
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
Output the elements of S modulo m in the same way as A is given.
Sample Input
2 2 4 0 1 1 1
Sample Output
1 2 2 3
Source
POJ Monthly--2007.06.03, Huang, Jinsong
解析:
用矩阵快速幂求A^k,然后把每一项矩阵相加即可,但这里的k很大,会超时,需要二分求解。
k为偶数:sum(k) = (1+A^(k/2)) *( A+A^2+……+A^(k/2)) = (1+A^(k/2)) * sum(k/2)
k为奇数:sum(k) = (1+A^((k-1)/2)) * sum(k/2) + A^k
k为奇数:sum(k) = (1+A^((k-1)/2)) * sum(k/2) + A^k
#include <iostream> #include<cstdio> #include<string.h> #include<cmath> using namespace std; typedef struct Node { int m[32][32]; }Matrax; Matrax a,per; int n,M; void init() { for(int i=0;i<n;i++) for(int j=0;j<n;j++) { scanf("%d",&a.m[i][j]); a.m[i][j]%=M;//不处理会超时 per.m[i][j]=(i==j);//单位矩阵,任何矩阵*单位矩阵=原矩阵 } } Matrax add(Matrax b,Matrax d)//两矩阵之和 { Matrax c; for(int i=0;i<n;i++) for(int j=0;j<n;j++) { c.m[i][j]=(b.m[i][j]+d.m[i][j])%M; } return c; } Matrax multi(Matrax b,Matrax d)//两矩阵相乘 { Matrax c; for(int i=0;i<n;i++) { for(int j=0;j<n;j++) { c.m[i][j]=0; for(int k=0;k<n;k++) { c.m[i][j]+=b.m[i][k]*d.m[k][j]; } c.m[i][j]%=M; } } return c; } Matrax pow(int k)//用矩阵快速幂求矩阵的k次方 { Matrax p,ans=per; p=a; while(k) { if(k&1) ans=multi(ans,p); k>>=1; p=multi(p,p); } return ans; } Matrax MatraxSum(int k)//二分法 { if(k==1) return a; Matrax temp,b; temp=MatraxSum(k/2); b=pow(k/2); temp=add(temp,multi(temp,b)); if(k&1) { temp=add(temp,pow(k)); } return temp; } int main() { int k; while(~scanf("%d%d%d",&n,&k,&M)) { init(); Matrax c=MatraxSum(k); for(int i=0;i<n;i++) { for(int j=0;j<n-1;j++) { printf("%d ",c.m[i][j]); } printf("%d\n",c.m[i][n-1]); } } return 0; }