您的位置:首页 > 其它

poj3233(快速求幂+二分)

2013-08-06 17:27 399 查看
Matrix Power Series

Time Limit: 3000MSMemory Limit: 131072K
Total Submissions: 12322Accepted: 5255
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
快速求幂加二分,快速矩阵求幂可以套用模板;
对于二分:
当k为奇数时:S[7]=A^1+A^2+A^3……+A^7=A^4+(1+A^4)(A^1+A^2+A^3)=A^4+S[3]+A^4 *S[3];
当k为偶数时:S[6]=S[3]+S[3] *A^3;
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
const int maxn=50;
typedef struct
{
int m[maxn][maxn];
}Matrax;
Matrax a,per;
int n,m;
void inital()
{
int i,j;
for(i=0;i<n;i++)
{
for(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 a,Matrax b)
{
Matrax c;
int i,j;
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
c.m[i][j]=(a.m[i][j]+b.m[i][j])%m;
}
}
return c;
}
Matrax multi(Matrax a,Matrax b)
{
Matrax c;
int k,i,j;
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
c.m[i][j]=0;
for(k=0;k<n;k++)
c.m[i][j]+=((a.m[i][k])%m*(b.m[k][j])%m)%m;
}
c.m[i][j]%=m;
}
return c;
}
Matrax power(int k)
{
Matrax p,ans=per;
p=a;
while(k)
{
if(k&1)
{
ans=multi(ans,p);
k--;
}
else
{
k/=2;
p=multi(p,p);
}
}
return ans;
}
Matrax MatraxSum(int k)
{
if(k==1)    return a;
Matrax temp,b;
temp=MatraxSum(k/2);
if(k&1)
{
b=power(k/2+1);
temp=add(temp,multi(temp,b));
temp=add(temp,b);
}
else
{
b=power(k/2);
temp=add(temp,multi(temp,b));
}
return temp;
}
int main()
{
int k,i,j;
Matrax ans;
while(scanf("%d%d%d",&n,&k,&m)==3)
{
inital();
ans=MatraxSum(k);
for(i=0;i<n;i++)
{
for(j=0;j<n-1;j++)
printf("%d ",ans.m[i][j]);
printf("%d\n",ans.m[i][j]);
}
}
return 0;
}


内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: