矩阵快速幂 poj3070 3233 3735 3150

一、矩阵的基础知识

1.结合性 （AB)C=A(BC）．

2.对加法的分配性 （A+B)C=AC+BC，C(A+B）=CA+CB ．

3.对数乘的结合性 k(AB）=（kA)B =A(kB）．

4.关于转置 (AB)'=B'A'．

一个矩阵就是一个二维数组，为了方便声明多个矩阵，我们一般会将矩阵封装一个类或定义一个矩阵的结构体，我采用的是后者。

最特殊的矩阵应该就是单位矩阵e了，它的对角线的元素为1，非对角线元素为0。一个n*n的矩阵的0次幂就是单位矩阵。

若A为n×k矩阵，B为k×m矩阵，则它们的乘积AB(有时记做A·B)将是一个n×m矩阵。其乘积矩阵AB的第i行第j列的元素为：

一般矩阵乘法采用朴素的O(n^3)的算法，但是对于一些比较稀疏的矩阵（就是矩阵中0比较多），对于这样的矩阵我们可以采用矩阵的优化，这个算法也适用于一般的矩阵，0特别多时，复杂度可能会降低到O(n^2)，实现如下：

还要注意的是，我们要尽可能的减少取模运算，因为取模的复杂度很高，这样我们就可以节约时间了。

矩阵加法就是简单地将对应的位置的两个矩阵的元素相加。

我们一般考虑的是n阶方阵之间的乘法以及n阶方阵与n维向量（把向量看成n×1的矩阵）的乘法。矩阵乘法最重要的性质就是满足结合律，同时它另一个很重要的性质就是不满足交换率，这保证了矩阵的幂运算满足快速幂取模（A^k % MOD）算法，矩阵快速幂其实就是二分指数，避免重复的计算。我们可以采用递归的方式很容易的写出来，但是当指数比较大，或者矩阵比较大得时候，我们就会出现栈溢出的状况，不断RE（我就被坑过）。所以还是写成迭代的方式比较好。

制作矩阵图一般要遵循以下几个步骤： 

1、列出质量因素： 

2、把成对对因素排列成行和列，表示其对应关系 

3、选择合适的矩阵图类型 

4、在成对因素交点处表示其关系程度，一般凭经验进行定性判断，可分为三种：关系密切、关系较密切、关系一般（或可能有关系)，并用不同符号表示 

5、根据关系程度确定必须控制的重点因素 

6、针对重点因素作对策表。

二、矩阵快速幂的应用

7、poj3070 是求解菲波那切数列，f(n)=f(n-1)+f(n-2),如果我们一个个递推求解，当n特别大的时候复杂度就会变的很高，对于f(n)= a*f(n-1)+b*f(n-2),在矩阵运算中我们会发现这样一组公式：

到知道这个公式后我们就采用矩阵快速幂的方法可以求解f(n)

#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
struct mat{
int at[2][2];
};
mat d;
int n,mod;
mat mul(mat a,mat b)
{
mat t;
memset(t.at,0,sizeof(t.at));
for(int i=0;i<n;++i)
{
for(int k=0;k<n;++k)
{
if(a.at[i][k])
for(int j=0;j<n;++j)
{
t.at[i][j]+=a.at[i][k]*b.at[k][j];
if(t.at[i][j]>=mod){t.at[i][j]%=mod;}
}
}
}
return t;
}
mat expo(mat p,int k)
{
if(k==1)return p;
mat e;
memset(e.at,0,sizeof(e.at));
for(int i=0;i<n;++i){e.at[i][i]=1;}
if(k==0)return e;
while(k)
{
if(k&1)e=mul(p,e);
p=mul(p,p);
k>>=1;
}
return e;
}
int main()
{
n=2;mod=10000;
d.at[1][1]=0;
d.at[0][0]=d.at[1][0]=d.at[0][1]=1;
int k;
while(~scanf("%d",&k))
{
if(k==-1)break;
mat ret=expo(d,k);
int ans=ret.at[0][1]%mod;
printf("%dn",ans);
}
return 0;
}

2、poj3233题意：给出矩阵A，求S = A + A^2 + A^3 + … + A^k 二分和

#include <iostream>
#include <cstdio>
#include <cstring>
using namespace std;
#define LL long long
int n,m,k;
int MOD;
struct mat {
int at[40][40];
};
mat d;
mat mul(mat a, mat b)
{
mat ret;
memset(ret.at,0,sizeof(ret.at));
for (int i=0;i<n;++i)
{
for (int k=0;k<n;++k)
{
if(a.at[i][k])
for (int j=0;j<n;++j)
{
ret.at[i][j]+=a.at[i][k]*b.at[k][j];
if(ret.at[i][j]>=MOD){ret.at[i][j]%=MOD;}
}
}
}
return ret;
}
mat expo(mat a, int k)
{
if(k==1)return a;
mat e;
memset(e.at,0,sizeof(e.at));
for(int i=0;i<n;++i){e.at[i][i]=1;}
if(k==0)return e;
while(k)
{
if(k&1)e=mul(a,e);
a=mul(a,a);
k>>=1;
}
return e;
}
mat add(mat a,mat b)
{
mat t;
for(int i=0;i<n;++i)
{
for(int j=0;j<n;++j)
{
t.at[i][j]=(a.at[i][j]+b.at[i][j]);
if(t.at[i][j]>=MOD){t.at[i][j]%=MOD;}
}
}
return t;
}
void print(mat ans)
{
for(int i=0;i<n;++i)
{
for(int j=0;j<n;++j)
{
if(j==0){printf("%d",ans.at[i][j]);continue;}
printf(" %d",ans.at[i][j]);
}
printf("n");
}
}
mat sum(int k)
{
if(k==1){return d;}
if(k&1)
{
return add(sum(k-1),expo(d,k));
}
else
{
mat s=sum(k>>1);
return add(s,mul(s,expo(d,k>>1)));
}
}
int main()
{
while(~scanf("%d%d%d",&n,&k,&m))
{
MOD=m;
mat ans,t;
for(int i=0;i<n;++i)
{
for(int j=0;j<n;++j)
{
scanf("%d",&d.at[i][j]);
if(d.at[i][j]>=m)
{
d.at[i][j]%=m;
}
}
}
ans=sum(k);
print(ans);
}
return 0;
}

3、poj3735

   题意：有n只猫咪，开始时每只猫咪有花生0颗，现有一组操作，由下面三个中的k个操作组成：

   1. g i 给i只猫咪一颗花生米

   2. e i 让第i只猫咪吃掉它拥有的所有花生米

   3. s i j 将猫咪i与猫咪j的拥有的花生米交换

   现将上述一组操作做m次后，问每只猫咪有多少颗花生？

分析：刚开始每只猫都没有花生，所以我们要在单位矩阵上构建矩阵。给第i只猫一个花生米，那么++met[0][i]，让第i只猫吃掉所有的花生米，就令第i列清空，喵咪i与猫咪j交换花生米，就令第i列和第j列互换。矩阵就这样构造完毕，操作m次，我们就可以矩阵快速幂计算了。

#include <iostream>
#include <cstring>
#include <cstdio>
#define LL long long
using namespace std;
struct met{
LL at[105][105];
};
met ret,d;
LL n,m,k;
met mul(met a,met b)
{
memset(ret.at,0,sizeof(ret.at));
for(int i=0;i<=n;++i)
{
for(int k=0;k<=n;++k)
{
if(a.at[i][k])
{
for(int j=0;j<=n;++j)
{
ret.at[i][j]+=a.at[i][k]*b.at[k][j];
}
}
}
}
return ret;
}
met expo(met a,LL k)
{
if(k==1) return a;
met e;
memset(e.at,0,sizeof(e.at));
for(int i=0;i<=n;++i){e.at[i][i]=1;}
if(k==0)return e;
while(k)
{
if(k&1)e=mul(e,a);
k>>=1;
a=mul(a,a);
}
return e;
}
int main()
{
while(~scanf("%lld%lld%lld",&n,&m,&k))
{
LL a,b;
char ch[5];
if(!n&&!k&&!m)break;
memset(d.at,0,sizeof(d.at));
for(int i=0;i<=n;++i)
{d.at[i][i]=1;}
while(k--)
{
scanf("%s",ch);
if(ch[0]=='g')
{
scanf("%lld",&a);
d.at[0][a]++;
}
else if(ch[0]=='e')
{
scanf("%lld",&a);
for(int i=0;i<=n;++i)
{
d.at[i][a]=0;
}
}
else {
scanf("%lld%lld",&a,&b);
for(int i=0;i<=n;++i)
{
LL t=d.at[i][a];
d.at[i][a]=d.at[i][b];
d.at[i][b]=t;
}
}
}
met ans=expo(d,m);
printf("%lld",ans.at[0][1]);
for(int i=2;i<=n;++i)
{
printf(" %lld",ans.at[0][i]);
}
printf("n");
}
return 0;
}

4、poj3150题目大意：给定n（1<=n<=500）个数字和一个数字m，这n个数字组成一个环(a0,a1.....an-1)。如果对ai进行一次d-step操作，那么ai的值变为与ai的距离小于d的所有数字之和模m。求对此环进行K次d-step（K<=10000000）后这个环的数字会变为多少。

分析：首先我们要构造矩阵，我们会得到一个500*500的矩阵，那么代码的复杂度就会变成O(log(k)*n^3)，很明显这么高的复杂度会超时的。但是我们发现这个矩阵是一个循环矩阵， 第i行都是第i-1行，右移一位得到的，即a[i][j]=a[i-1][j-1]。很容易我们就可以发现循环矩阵a和循环矩阵b的乘积矩阵c，c[i][j]=sum(a[i][k]*b[k][j])=sum(a[i-1][k-1]*b[j-1][k-1])=c[i-1][j-1]。那么矩阵c也是一个循环矩阵，在做矩阵乘法的时候我们只需要算出第一行的值，其余行直接右移就可以得到，那么算法的复杂度就会变为O(log(k)*n^2)。还需注意的是对于数据范围会超int，要用long long，还有由于矩阵太大了，在函数中申请不了那么大得空间，所以采用指针的方法去写函数。

#include <iostream>
#include <cstdio>
#include <cstring>
#define LL long long
using namespace std;
const int maxn=502;
int n,m,d,k;
LL tmp[maxn][maxn],e[maxn][maxn],c[maxn][maxn];
void mul(LL a[][maxn],LL b[][maxn])
{
memset(c,0,sizeof(c));
for(int k=0;k<n;++k)
{
if(a[0][k])
for(int j=0;j<n;++j)
{
c[0][j]+=a[0][k]*b[k][j];
if(c[0][j]>=m){c[0][j]%=m;}
}
}
for(int i=1;i<n;++i)
{
for(int j=0;j<n;++j)
{
c[i][j]=c[i-1][(j-1+n)%n];
}
}
for(int i=0;i<n;++i)
{
for(int j=0;j<n;++j)
{
b[i][j]=c[i][j];
}
}
}
void expo(LL a[][maxn],int k)
{
if(k==1){
for(int i=0;i<n;++i)
{
for(int j=0;j<n;++j)
{
e[i][j]=a[i][j];
}
}
return;
}
memset(e,0,sizeof(e));
for(int i=0;i<n;++i){e[i][i]=1;}
while(k)
{
if(k&1){mul(a,e);}
mul(a,a);
k>>=1;
}
}
int main()
{
LL dat[maxn];
scanf("%d%d%d%d",&n,&m,&d,&k);
for(int i=0;i<n;++i)
{
scanf("%lld",&dat[i]);
tmp[0][i]=0;
}
tmp[0][0]=1;
for(int i=1;i<=d;++i)
{
tmp[0][i]=tmp[0][n-i]=1;
}
for(int i=1;i<n;++i)
{
for(int j=0;j<n;++j)
{
tmp[i][j]=tmp[i-1][(j-1+n)%n];
}
}
expo(tmp,k);
LL ans[maxn];
memset(ans,0,sizeof(ans));
for(int i=0;i<n;++i)
{
for(int j=0;j<n;++j)
{
ans[i]+=e[i][j]*dat[j];
if(ans[i]>=m){ans[i]%=m;}
}
}
printf("%lld",ans[0]);
for(int i=1;i<n;++i)
{
printf(" %lld",ans[i]);
}
printf("n");
return 0;
}

对于这道题，网上还有一段神代码，在这里同样学习一下

#include <iostream>
#include <cstdio>
#include <cstring>
#define LL long long
using namespace std;
int n,m,d,k;
void mul(LL a[],LL b[])
{
int i,j;
LL c[501];
for(i=0;i<n;++i)for(c[i]=j=0;j<n;++j)c[i]+=a[j]*b[i>=j?(i-j):(n+i-j)];
for(i=0;i<n;b[i]=c[i++]%m);
}
LL init[501],tmp[501];
int main()
{
int i,j;
scanf("%d%d%d%d",&n,&m,&d,&k);
for(i=0;i<n;++i)scanf("%lld",&init[i]);
for(tmp[0]=i=1;i<=d;++i)tmp[i]=tmp[n-i]=1;
while(k)
{
if(k&1)mul(tmp,init);
mul(tmp,tmp);
k>>=1;
}
for(i=0;i<n;++i)if(i)printf(" %lld",init[i]);else printf("%lld",init[i]);
printf("n");
return 0;
}