poj 1995 Matrix Power Series 二分+矩阵快速幂

Matrix Power Series

Time Limit: 3000MS		Memory Limit: 131072K
Total Submissions: 21337		Accepted: 8944

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

题意：给出一个n*n矩阵A，求表达式S = A + A2 + A3 +
… + Ak，并对m取模。 n (n ≤
30), k (k ≤
109) and m (m <
104). 

解法：由于k太大，对于任意次幂，不能一一求出结果再相加。如果A是一个数，而不是矩阵那么必有循环节。但是A是矩阵不能从循环节入手，但是一定只需要计算部分A的次幂，容易想到计算一部分，再用这部分计算出下一部分，结果发现确实可以。

故利用二分求解即可。

#include<cstdio>
#include<string>
#include<cstring>
#include<iostream>
#include<cmath>
#include<algorithm>
#include<vector>
using namespace std;

#define all(x) (x).begin(), (x).end()
#define for0(a, n) for (int (a) = 0; (a) < (n); (a)++)
#define for1(a, n) for (int (a) = 1; (a) <= (n); (a)++)
#define mes(a,x,s)  memset(a,x,(s)*sizeof a[0])
#define mem(a,x)  memset(a,x,sizeof a)
#define ysk(x)  (1<<(x))
typedef long long ll;
typedef pair<int, int> pii;
const int INF =0x3f3f3f3f;
const int maxn= 30   ;

int n,K;
ll mod;
struct Matrix
{
ll m[maxn+3][maxn+3];
Matrix(){}
void out()
{
for0(i,n) for0(j,n)  printf("%lld%c",m[i][j],j==n-1?'\n':' ');
}
}A,E;

void init()
{
for0(i,n) for0(j,n)
{
cin>>A.m[i][j];
A.m[i][j]%=mod;
E.m[i][j]=i==j?1%mod:0;
}
}

Matrix operator*(Matrix a,Matrix b)
{
Matrix c;
for0(i,n) for0(j,n)
{
c.m[i][j]=0;
for0(k,n)
{
c.m[i][j]+=a.m[i][k]*b.m[k][j];
}
c.m[i][j]%=mod;
}
return c;
}

Matrix operator+(Matrix a,Matrix b)
{
Matrix c;
for0(i,n) for0(j,n)
{
c.m[i][j]=a.m[i][j]+b.m[i][j];
c.m[i][j]%=mod;
}
return c;
}
Matrix power(Matrix a,ll k)
{
Matrix ans=E;
while(k)
{
if(k&1) ans=ans*a;
k>>=1;
a=a*a;
}
return ans;
}

Matrix  solve(int num)
{
if(num==1)  return A;
Matrix ans;
if(num&1)
{
Matrix c=solve(num/2);
ans=c+power(A,num/2 +1 )*c+power(A,num/2+1);
}
else
{
Matrix c=  solve(num/2);
ans=  c+ power(A,num/2)*c;
}
return ans;
}
int main()
{
std::ios::sync_with_stdio(false);
while(cin>>n>>K>>mod)
{
init();
Matrix ans=solve(K);
ans.out();
}
return 0;
}