矩阵快速幂 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("%d\n",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 longint 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 longusing 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 longusing 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 longusing 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;}