目的:快速计算矩阵A的b次幂,即计算A^b;
代码如下:
typedef struct
{
int m[MAXN][MAXN];
}Matrix;
Matrix a,per;
int n,M;
void init()///给矩阵赋初值
{
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);///如果i==j,返回1,否则,返回0.
}///per.m[][]为对角线元素为1,其余为0的矩阵
}
Matrix multi(Matrix a,Matrix b)
{
Matrix c;
int i,j,k;
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]*b.m[k][j];
}
c.m[i][j]%=M;
}
return c;
}
Matrix power(int k)///矩阵的快速幂
{
Matrix c,p;
Matrix ans=per;///p第一次与ans相乘后,得到的结果还是p本身不变
p=a;
while(k)///二分法
{
if(k&1)
{
ans=multi(ans,p);
k--;
}
k>>=1;
p=multi(p,p);
}
return ans;
}
例题一 Matrix Power Series(pku 3233 )
给定一个n*n的矩阵和一个整数k,要求计算S=A+A^2+A^3.......+A^k;
输入:输入仅包含一组测试数据,第一行包含三个整数n(n<=30),k(k<=10^9),m(m<10^4);接下来的n行每行包含n个n个小于32 768的非负整数位矩阵n*n个元素
输出:S%m;
代码如下:
#include<stdio.h>///把G++改成C++就AC
#include<iostream>
#include<string.h>
#define MAXN 101
using namespace std;
typedef struct///自定义结构体
{
int m[MAXN][MAXN];
}Matrix;
Matrix a,per;
int n,m;
void init()///初始化
{
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);
}
}
Matrix add(Matrix a,Matrix b)///矩阵相加
{
Matrix 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;
}
Matrix multi(Matrix a,Matrix b)///矩阵相乘
{
Matrix 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]*b.m[k][j];
}
c.m[i][j]%=m;
}
return c;
}
Matrix power(int k)///矩阵的k次幂
{
Matrix c,p,ans=per;
p=a;
while(k)
{
if(k&1)
{
ans=multi(ans,p);
k--;
}
k>>=1;
p=multi(p,p);
}
return ans;
}
Matrix MatrixSum(int k)///矩阵之和,递归求解
{ ///求的是A+A^2+A^3...A^k
if(k==1)
return a;
Matrix temp,b;
temp=MatrixSum(k/2);
if(k&1)///若k为奇数.
{
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 i,j,k;
while(scanf("%d%d%d",&n,&k,&m)!=EOF)
{
init();
Matrix ans=MatrixSum(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;
}