矩阵快速幂和普通的快速幂差不多，只不过写起来比较麻烦一点，需要重载*运算符。
模板：
struct mat
{
int m[maxn][maxn];
}unit; mat operator * (mat a,mat b)
{
mat ret;
ll x;
for(int i=;i < n;i++)
for(int j=;j < n;j++)
{
x = ;
for(int k=;k < n;k++)
x += mod((ll)a.m[i][k]*b.m[k][j]);
ret.m[i][j] = mod(x);
}
return ret;
} void init_unit() //初始化单位矩阵
{
for(int i=;i < maxn;i++)
unit.m[i][i] = ;
return;
} mat pow_mat(mat a,ll n)
{
mat ret = unit;
while(n)
{
if(n&) ret = ret*a;
a = a*a;
n >>= ;
}
return ret;
}
例题：POJ3070
用矩阵快速幂求fib并取模10000
有这个定理就很好求了：
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int INF = 0x3f3f3f3f;
const int moder = ;
const int maxn = ;
#define mod(x) ((x)%moder)
int n = ; struct mat
{
int m[maxn][maxn];
}unit; mat operator * (mat a,mat b)
{
mat ret;
ll x;
for(int i=;i < n;i++)
for(int j=;j < n;j++)
{
x = ;
for(int k=;k < n;k++)
x += mod((ll)a.m[i][k]*b.m[k][j]);
ret.m[i][j] = mod(x);
}
return ret;
} void init_unit()
{
for(int i=;i < maxn;i++)
unit.m[i][i] = ;
return;
} mat pow_mat(mat a,ll n)
{
mat ret = unit;
while(n)
{
if(n&) ret = ret*a;
a = a*a;
n >>= ;
}
return ret;
} int main()
{
ll p;
init_unit();
while(cin >> p)
{
if(p == -) break;
mat a;
a.m[][] = ;
a.m[][] = ;
a.m[][] = ;
a.m[][] = ;
a = pow_mat(a,p);
cout << a.m[][] << endl;
}
return ;
}
要注意的是maxn开小一点，不然本地会炸。