矩阵的幂其实就是多次相乘，但是用最普通的方法求往往是超时的这就用到了快速幂的思想
还是用了二分的思想：
下面Down了网友的模板（可以结合者后面的具体题目理解~）
//HOJ 3493
/*===================================*/
|| 快速幂（quickpow）模板
|| P 为等比，I 为单位矩阵
|| MAX 要初始化！！！！
||
/*===================================*/
/*****************************************************/
#include <cstdio>
const int MAX = 3;
typedef struct{
int m[MAX][MAX];
} Matrix;
Matrix P = {5,-7,4,
1,0,0,
0,1,0,
};
Matrix I = {1,0,0,
0,1,0,
0,0,1,
};
Matrix matrixmul(Matrix a,Matrix b) //矩阵乘法
{
int i,j,k;
Matrix c;
for (i = 0 ; i < MAX; i++)
for (j = 0; j < MAX;j++)
{
c.m[i][j] = 0;
for (k = 0; k < MAX; k++)
c.m[i][j] += (a.m[i][k] * b.m[k][j])%9997;
c.m[i][j] %= 9997;
}
return c;
}
Matrix quickpow(long long n)
{
Matrix m = P, b = I;
while (n >= 1)
{
if (n & 1)
b = matrixmul(b,m);
n = n >> 1;
m = matrixmul(m,m);
}
return b;
}
/*************************************/
int main()
{
Matrix re;
int f[3] = {2,6,19};
long long n;
while (scanf("%I64d",&n) && n != 0)
{
if (n == 1)
printf("1\n");
else if (n <= 4)
printf("%d\n",f[n-2]);
else {
re = quickpow(n - 4);
printf("%d\n",(((re.m[0][0]*f[2])
+ (re.m[0][1]*f[1]) + (re.m[0][2]*f[0])) %9997 + 9997) % 9997);
}
}
return 0;
}
下面举一个用快速矩阵幂求解的题目：
Matrix Power Series
题目描述：
Description
Given a n × n matrix A and a positive integer k, find the sumS = A + A2 + A3 + … + Ak.
Input
The input contains exactly one test case. The first line of input contains three positive integersn (n ≤ 30), k (k ≤ 109) and m (m < 104). Then follown lines each containing n nonnegative integers below 32,768, givingA’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
矩阵快速幂。首先我们知道 A^x 可以用矩阵快速幂求出来（具体可见poj 3070）。其次可以对k进行二分，每次将规模减半，分k为奇偶两种情况，如当k = 6和k = 7时有：
偶数: F[n]=F[n/2]+F[n/2]*An/2
（是不是和求数的快速幂很像？？？）
#include<iostream>
#include <cstdio>
using namespace std;
const int MAX = 32;
//int count=0;
struct Matrix //用于相同的类型的数组的开辟，比用单独的名称好
{
int v[MAX][MAX];
};
int n, k, M;
Matrix mtAdd(Matrix A, Matrix B) // 求矩阵 A + B
{
int i, j;
Matrix C;
for(i=0; i<n; i++)
for(j=0; j<n; j++)
C.v[i][j]=(A.v[i][j] + B.v[i][j]) % M;
// count++;
return C;
}
Matrix mtMul(Matrix A, Matrix B) // 求矩阵 A*B
{
int i, j, k;
Matrix C;
for(i = 0; i < n; i ++)
for(j = 0; j < n; j ++)
{
C.v[i][j] = 0;
for(k = 0; k < n; k ++)
C.v[i][j] = (A.v[i][k] * B.v[k][j] + C.v[i][j]) % M;
// count++;
}
return C;
}
Matrix mtPow(Matrix A, int k) // 求矩阵 A^k
{
if(k == 0) {
memset(A.v, 0, sizeof(A.v));
for(int i = 0; i < n; i ++)
A.v[i][i]=1;
// count++;
return A;
}
//count++;
if(k == 1) return A;
Matrix C = mtPow(A, k / 2);
if(k % 2 == 0)
return mtMul(C, C);
else
return mtMul(mtMul(C, C), A);
}
Matrix mtCal(Matrix A, int k) // 求S (k) = A + A2 + A3 + … + Ak
{
//count++;
if(k == 1) return A;
Matrix B = mtPow(A, (k+1) / 2);
Matrix C = mtCal(A, k / 2);
if(k % 2 == 0)
return mtMul( mtAdd(mtPow(A, 0), B), C ); // 如S(6) = ( 1 （单位矩阵） + A^3 ) * S(3)。
else
return mtAdd( A, mtMul(mtAdd(A, B), C) ); // 如S(7) = A + (A + A^4) * S(3)
}
int main()
{
int i, j;
Matrix A;//输入的原矩阵
cin>>n>>k>>M;
for(i=0; i<n; i++)
for(j=0;j<n;j++)
scanf("%d",&A.v[i][j]);
A=mtCal(A, k);
for(i=0; i<n; i++){
for(j=0; j<n; j++){
printf("%d ",A.v[i][j]);
}
printf("\n");
}
//cout<<count;
return 0;
}
还有一个，估计其实质是一样的但是现在还不是很懂，以后弄懂~~~
因为S可以看成S=A(I+A(I+A(I+...A(I+A)))) (I是单位矩阵)
拿k=3举例S=A(I+A(I+A))
那么我们想，可不可以构造一个矩阵T使得T*T（因为是k次幂）这样乘下去每次可以得到A*(A+I)
那么肯定T有个两个元素就是A与I
那么假设：T=A I
I I
那么T=T*T=A*A+I*I A*I+I*I
A*I+I*I I*I+I*I
这样存在一个I*（A+I）的式子 ，当T再乘以T的时候会出现A（A+I）
这个时候我们可以简化将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次幂就可以了
而k如此庞大所以需要二分来对T求k+1次幂
对T求快速幂：
首先我们知道A^19=A^16*A^2*A^1，因为19=B(10011)
那么就这样，拿A总是去和本身去乘，那么就可以取到 A A^2 A^4 A^8 A^16...
然后问A A^2 A^4 A^8 A^16 ...这么多项取与不取的抉择
其实就是一个求二进制的这样一个过程
19%2=1 那么A取进来
19/2=9
9%2=1 那么A^2取进来
9/2=4
4%2=0 那么A^4不用取进来
4/2=2
2%2=0 那么A^8不用去进来
2/2=1
1%2=1 那么A^16取进来
1/2=0 计算完毕
这个就是用模板做出来的~~~
抄了一点别人的代码~~~呵呵~~~
#include <iostream>
#include <string.h>
using namespace std;
#define MAXV 66
struct{
int r,c; //c行r列
int mat[MAXV][MAXV];///矩阵中再放入矩阵
}Matrix;
Matrix ans,cnt;
int n,k,m;
void Input(){
int i,j;
/*构造B矩阵*/
memset(cnt.mat,0,sizeof(cnt.mat));
memset(ans.mat,0,sizeof(ans.mat));
for(i=0;i<n;i++){
for(j=0;j<n;j++)
scanf("%d",&cnt.mat[i][j]);
}
for(i=0;i<n;i++){
cnt.mat[i+n][i+n]=cnt.mat[i][i+n]=1;
ans.mat[i][i]=ans.mat[i+n][i+n]=1;
}
/*对矩阵B和B^（k+1）初始化*/
cnt.c=cnt.r=2*n;
ans.c=ans.r=2*n;
}
Matrix MatrixMul(Matrix x,Matrix y){ //矩阵乘法
Matrix t;
int i,j,v;
memset(t.mat,0,sizeof(t.mat));
t.c=x.c;
t.r=y.r;
for(i=0;i<t.c;i++)
for(j=0;j<t.r;j++){
for(v=0;v<x.r;v++)
t.mat[i][j]=t.mat[i][j]+((x.mat[i][v]*y.mat[v][j])%m);
t.mat[i][j]=t.mat[i][j]%m;
}
return t;
}
void Binary(){
//二分快速幂
k++;
while(k){
if(k & 1)
ans=MatrixMul(ans,cnt);
cnt=MatrixMul(cnt,cnt);
k=k>>1;
}
}
void Output(){
/*输出的时候要减去一个单位矩阵*/
int i,j;
for(i=0;i<n;i++){
for(j=0;j<n;j++)
if(i!=j)
printf("%d ",ans.mat[i][j+n]);
else
printf("%d ",ans.mat[i][j+n]-1);
printf("\n");
}
}
int main(){
while(~scanf("%d%d%d",&n,&k,&m)){
Input();
Binary();
Output();
}
return 0;
}
另外，本体还可以用类似于数列的前N项和来求解，这样的话用一次对幂求解一套公式就出来了，别说现在前N项和记不起来了~~~
A(1-A^q)/(1-q) (这里的1应该是单位向量了)
