poj 3233 Matrix Power Series （矩阵快速幂）

这里介绍两种做法：

一：思路：

做这题时用的是类似二分的方法做，仔细观察可以发现

Sk = A + A2 + A3 + … + Ak   

    =(1+Ak/2)*(A + A2 + A3 + … + Ak/2  )+{Ak}

    =(1+Ak/2)*(Sk/2 )+{Ak}

当k为偶数时不要大括号里面的数；

代码：//800ms
#include<iostream>
#include<cstdlib>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
int n , m ;
struct matrix {//矩阵
int a[30][30];
void init( )
{//初始化
int i ;
memset ( a, 0 ,sizeof ( a ) ) ;
for ( i = 0 ; i < 30 ; i ++ )
a[i][i] = 1 ;
}
};
void print ( matrix s )
{//打印一个矩阵
int i , j ;
for ( i = 0 ; i < n ; i ++ )
{
for ( j = 0 ; j < n ; j ++ )
{
if ( j )
cout<<' ';
cout<<s.a[i][j]%m;
}
cout<<endl;
}
}
matrix matrixadd ( matrix a , matrix b )
{//矩阵加法，计算时记得要取模，避免超int
int i , j ;
matrix c ;
for ( i = 0 ; i < n ; i ++ )
for ( j = 0 ; j < n ; j ++ )
c.a[i][j]=((a.a[i][j]+b.a[i][j])%m);
return c ;
}
matrix matrixmul ( matrix a , matrix b )
{//矩阵乘法，计算时记得要取模，避免超int
int i , j , k ;
matrix c ;
for ( i = 0 ; i < n ; i ++ )
{
for ( j = 0 ; j < n ; j ++ )
{
c.a[i][j]=0;
for ( k = 0 ; k < n ; k ++ )
c.a[i][j] +=((a.a[i][k]*b.a[k][j])%m) ;
c.a[i][j] %= m ;
}
}
return c ;
}
matrix mul ( matrix s , int k )
{//矩阵的k次方，快速幂
matrix ans ;
ans .init () ;
while ( k >= 1 )
{
if ( k & 1 )
ans = matrixmul ( ans , s ) ;
k = k >> 1 ;
s = matrixmul ( s , s ) ;
}
return ans ;
}
matrix sum ( matrix s , int k )
{//举证前k项求和
if ( k == 1 )
return s ;
matrix tmp ;//用来保存答案
tmp.init();//初始化
tmp = matrixadd ( tmp , mul ( s , k >> 1 ) ); //计算1+A^(k/2)
tmp = matrixmul ( tmp , sum ( s , k >> 1 ) ) ; //计算(1+A^(k/2))*(A + A^2 + A^3 + … + A^(k/2) )
if ( k&1 )//考虑是否要+A^k
tmp = matrixadd ( tmp , mul ( s , k ) ) ;
return tmp ;//返回前n项的值
}

int main()
{
int k ;
while ( cin >> n >> k >> m )
{
int i , j ;
matrix s ;
for ( i = 0 ; i < n ; i ++ )
for ( j = 0 ; j < n ; j ++ )
cin >> s.a[i][j] ;
s = sum ( s , k ) ;
print(s);
}
return 0;
}二：
构造一个矩阵T

T={A  I}  

      0   I

所以T*T={A*A   A*I+I*I}

                   0           I

那么容易得到T^(K+1)={A^(K+1)           I+A+A^2+A^3+...+A^K}

                                         0                                I

这样我们只需要算T的k+1次幂就可以了；

代码：#include <cstdio>
#include <cstring>
#include <iostream>
#include <cmath>
#include <algorithm>
const int maxn = 66;
using namespace std;

struct mat{
int v[maxn][maxn];
mat() {
memset(v, 0, sizeof(v));
}
};
int n, m;

mat matrix_mul(mat a, mat b) {
mat c;
int i, j, k;
for(i = 0; i < 2*n; i++)
for(j = 0; j < 2*n; j++)
for(k = 0; k < 2*n; k++) {
c.v[i][k] += (a.v[i][j]*b.v[j][k])%m;
c.v[i][k] %= m;
}
return c;
}

mat matrix_mi(mat p, int k) {
mat t;
for(int i = 0; i < 2*n; i++)
t.v[i][i] = 1;
while(k) {
if (k & 1)
t = matrix_mul(t, p);
k >>= 1;
p = matrix_mul(p, p);
}
return t;
}

int main() {
int k;
cin>>n>>k>>m;
mat p, t;
for(int i = 0; i < n; i++) {
p.v[i][i+n] = p.v[i+n][i+n] = 1;
for(int j = 0; j < n; j++)
cin>>p.v[i][j];
}
t = matrix_mi(p, k + 1);
for(int i = 0; i < n; i++) {
for(int j = 0; j < n; j++)
if (i != j)
printf("%d ", t.v[i][j+n]);
else printf("%d ", (t.v[i][j+n] + m - 1)%m);
printf("\n");
}
return 0;
}