十个利用矩阵乘法解决的经典题目(个人小结)

看了Matrix67大神的博客做的一个小结。也让我对矩阵有了一个全新的认识。

传送阵

经典题目1 给定n个点，m个操作，构造O(m+n)的算法输出m个操作后各点的位置。操作有平移、缩放、翻转和旋转

    这里的操作是对所有点同时进行的。其中翻转是以坐标轴为对称轴进行翻转（两种情况），旋转则以原点为中心。如果对每个点分别进行模拟，那么m个操作总共耗时O(mn)。利用矩阵乘法可以在O(m)的时间里把所有操作合并为一个矩阵，然后每个点与该矩阵相乘即可直接得出最终该点的位置，总共耗时O(m+n)。假设初始时某个点的坐标为x和y，下面5个矩阵可以分别对其进行平移、旋转、翻转和旋转操作。预先把所有m个操作所对应的矩阵全部乘起来，再乘以(x,y,1)，即可一步得出最终点的位置。

     


这些就是矩阵的简单操作，所以记住就行了。

经典题目2 给定矩阵A，请快速计算出A^n（n个A相乘）的结果，输出的每个数都mod p。

    由于矩阵乘法具有结合律，因此A^4 = A * A * A * A = (A*A) * (A*A) = A^2 * A^2。我们可以得到这样的结论：当n为偶数时，A^n = A^(n/2) * A^(n/2)；当n为奇数时，A^n = A^(n/2) * A^(n/2) * A （其中n/2取整）。这就告诉我们，计算A^n也可以使用二分快速求幂的方法。例如，为了算出A^25的值，我们只需要递归地计算出A^12、A^6、A^3的值即可。根据这里的一些结果，我们可以在计算过程中不断取模，避免高精度运算。

此种题目包含于下题中，所以看下题例题。

经典题目3 POJ3233 

    题目大意：给定矩阵A，求A + A^2 + A^3 + ... + A^k的结果（两个矩阵相加就是对应位置分别相加）。输出的数据mod m。k<=10^9。这道题两次二分，相当经典。首先我们知道，A^i可以二分求出。然后我们需要对整个题目的数据规模k进行二分。比如，当k=6时，有：A + A^2 + A^3 + A^4 + A^5 + A^6 =(A + A^2 + A^3) + A^3*(A + A^2 + A^3)

    应用这个式子后，规模k减小了一半。我们二分求出A^3后再递归地计算A + A^2 + A^3，即可得到原问题的答案。

附上：个人代码

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<string>
#include<cmath>
#include<set>
#include<vector>
#include<stack>
#define mem(a,b) memset(a,b,sizeof(a))
#define FOR(a,b,i) for(i=a;i<=b;++i)
#define For(a,b,i) for(i=a;i<b;++i)
#define N 1000000007
using namespace std;
inline void RD(int &ret)
{
char c;
do
{
c=getchar();
}
while(c<'0'||c>'9');
ret=c-'0';
while((c=getchar())>='0'&&c<='9')
{
ret=ret*10+(c-'0');
}
}
inline void OT(int a)
{
if(a>=10)
{
OT(a/10);
}
putchar(a%10+'0');
}
__int64 m;
int n;
__int64 tx[33][33],ty[66][66],tz[66][66];
void mat(__int64 a[66][66],__int64 b[66][66])//矩阵乘
{
__int64 d[66][66];
int i,j,l;
mem(d,0);
For(0,2*n,i)
{
For(0,2*n,j)
{
For(0,2*n,l)
{
d[i][j]=(d[i][j]+a[i][l]*b[l][j]%m)%m;
}
}
}
For(0,2*n,i)
{
For(0,2*n,j)
{
a[i][j]=d[i][j];
}
}
}
void g(__int64 x)//矩阵快速幂，用的都是自己用的模板
{
mem(ty,0);
mem(tz,0);
int i,j;
For(0,n,i)
{
For(0,n,j)
{
tz[i+n][j+n]=ty[i][j+n]=tx[i][j];
}
}
For(0,2,j)
{
For(0,n,i)
{
tz[i+j*n][i]=1;
}
}
while(x>0)
{
if(x%2==1)
{
mat(ty,tz);
}
mat(tz,tz);
x/=2;
}
}
int main()
{
__int64 k;
int i,j;
RD(n);
scanf("%I64d%I64d",&k,&m);
For(0,n,i)
{
For(0,n,j)
{
scanf("%I64d",&tx[i][j]);
tx[i][j]%=m;
}
}
g(k);
For(0,n,i)
{
printf("%d",ty[i][0]);
For(1,n,j)
{
printf(" %d",ty[i][j]);
}
printf("\n");
}
return 0;
}

经典题目4 VOJ1049

    题目大意：顺次给出m个置换，反复使用这m个置换对初始序列进行操作，问k次置换后的序列。m<=10, k<2^31。

    首先将这m个置换“合并”起来（算出这m个置换的乘积），然后接下来我们需要执行这个置换k/m次（取整，若有余数则剩下几步模拟即可）。注意任意一个置换都可以表示成矩阵的形式。置换k/m次就相当于在前面乘以k/m个这样的矩阵。我们可以二分计算出该矩阵的k/m次方，再乘以初始序列即可。做出来了别忙着高兴，得意之时就是你灭亡之日，别忘了最后可能还有几个置换需要模拟。
附上：本人代码

VOJ传送阵

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<string>
#include<cmath>
#include<set>
#include<vector>
#include<stack>
#define mem(a,b) memset(a,b,sizeof(a))
#define FOR(a,b,i) for(i=a;i<=b;++i)
#define For(a,b,i) for(i=a;i<b;++i)
//#define N 1000000007
using namespace std;
inline void RD(int &ret)
{
char c;
do
{
c=getchar();
}
while(c<'0'||c>'9');
ret=c-'0';
while((c=getchar())>='0'&&c<='9')
{
ret=ret*10+(c-'0');
}
}
inline void OT(int a)
{
if(a>=10)
{
OT(a/10);
}
putchar(a%10+'0');
}
int n;
int tx[101][101],ty[101][101],tz[101][101],tw[101][101];
void mat(int a[101][101],int b[101][101])//矩阵乘
{
int d[101][101];
int i,j,l;
mem(d,0);
FOR(1,n,i)
{
FOR(1,n,j)
{
FOR(1,n,l)
{
d[i][j]+=a[i][l]*b[l][j];
}
}
}
FOR(1,n,i)
{
FOR(1,n,j)
{
a[i][j]=d[i][j];
}
}
}
void g(int x)//矩阵快速幂，用的都是自己用的模板
{
mem(ty,0);
int i,j;
FOR(1,n,i)
{
ty[i][i]=1;
}
while(x>0)
{
if(x%2==1)
{
mat(ty,tz);
}
mat(tz,tz);
x/=2;
}
}
int main()
{
int m,k,i,j,l;
while(scanf("%d%d%d",&n,&m,&k)!=EOF)
{
mem(tx,0);
FOR(1,n,i)
{
tx[i][i]=1;
}
FOR(1,m,i)
{
FOR(1,n,j)
{
RD(tw[i][j]);
}
mem(tz,0);
FOR(1,n,j)
{
tz[j][tw[i][j]]=1;
}
mat(tz,tx);
FOR(1,n,l)
{
FOR(1,n,j)
{
tx[l][j]=tz[l][j];
}
}
}
g(k/m);
mem(tz,0);
FOR(1,n,i)
{
tz[i][1]=i;
}
mat(ty,tz);
FOR(1,k%m,i)
{
mem(tz,0);
FOR(1,n,j)
{
tz[j][tw[i][j]]=1;
}
mat(tz,ty);
FOR(1,n,l)
{
FOR(1,n,j)
{
ty[l][j]=tz[l][j];
}
}
}
FOR(1,n,i)
{
printf("%d",ty[i][1]);
if(i<n)
{
printf(" ");
}
}
printf("\n");
}
return 0;
}

经典题目5 《算法艺术与信息学竞赛》207页（2.1代数方法和模型，[例题5]细菌，版次不同可能页码有偏差）

    大家自己去看看吧，书上讲得很详细。解题方法和上一题类似，都是用矩阵来表示操作，然后二分求最终状态。这就是简单的矩阵操作和矩阵快速幂，我就不多写了，书上写得很详细。

经典题目6 给定n和p，求第n个Fibonacci数mod p的值，n不超过2^31

    根据前面的一些思路，现在我们需要构造一个2 x 2的矩阵，使得它乘以(a,b)得到的结果是(b,a+b)。每多乘一次这个矩阵，这两个数就会多迭代一次。那么，我们把这个2 x 2的矩阵自乘n次，再乘以(0,1)就可以得到第n个Fibonacci数了。不用多想，这个2 x 2的矩阵很容易构造出来：


这类题题解很多，属于简单题就不多写了。

经典题目7 VOJ1067

    我们可以用上面的方法二分求出任何一个线性递推式的第n项，其对应矩阵的构造方法为：在右上角的(n-1)*(n-1)的小矩阵中的主对角线上填1，矩阵第n行填对应的系数，其它地方都填0。

这题，有人说矩阵应该建反的，但我就是用Matrix67大神说的建图方法过的。。。

附上：本人代码

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<string>
#include<cmath>
#include<set>
#include<vector>
#include<stack>
#define mem(a,b) memset(a,b,sizeof(a))
#define FOR(a,b,i) for(i=a;i<=b;++i)
#define For(a,b,i) for(i=a;i<b;++i)
#define N 7777777
using namespace std;
inline void RD(int &ret)
{
char c;
do
{
c=getchar();
}
while(c<'0'||c>'9');
ret=c-'0';
while((c=getchar())>='0'&&c<='9')
{
ret=ret*10+(c-'0');
}
}
inline void OT(int a)
{
if(a>=10)
{
OT(a/10);
}
putchar(a%10+'0');
}
int n,k;
__int64 ty[11][11],tz[11][11],tw[11][11],f[11],s;
void mat(__int64 a[11][11],__int64 b[11][11]) //矩阵乘
{
__int64 d[11][11];
int i,j,l;
mem(d,0);
FOR(1,k,i)
{
FOR(1,k,j)
{
FOR(1,k,l)
{
d[i][j]=(d[i][j]+a[i][l]*b[l][j]%N)%N;
}
}
}
FOR(1,k,i)
{
FOR(1,k,j)
{
a[i][j]=d[i][j];
}
}
}
void g(int x)//矩阵快速幂，用的都是自己用的模板
{
int i,j;
FOR(1,k,i)
{
FOR(1,k,j)
{
ty[i][j]=tz[i][j];
}
}
while(x>0)
{
if(x%2==1)
{
mat(ty,tz);
}
mat(tz,tz);
x/=2;
}
}
int main()
{
int i,j;
RD(k);
RD(n);
mem(f,0);
f[1]=1;
FOR(2,k,i)
{
f[i]=1;
FOR(1,i-1,j)
{
f[i]=f[i]+f[j];
}
}

if(k==1)
{
printf("1\n");
}
else if(n<=k)
{
printf("%I64d\n",f
);
}
else
{
mem(tz,0);
FOR(2,k,i)
{
tz[i][i-1]=1;
}
FOR(1,k,i)
{
tz[i][k]=1;
}
g(n-k-1);
s=0;
FOR(1,k,i)
{
s=(s+(f[i]%N)*(ty[i][k]%N))%N;
}
printf("%I64d\n",s);
}
return 0;
}

经典题目8 给定一个有向图，问从A点恰好走k步（允许重复经过边）到达B点的方案数mod p的值

    把给定的图转为邻接矩阵，即A(i,j)=1当且仅当存在一条边i->j。令C=A*A，那么C(i,j)=ΣA(i,k)*A(k,j)，实际上就等于从点i到点j恰好经过2条边的路径数（枚举k为中转点）。类似地，C*A的第i行第j列就表示从i到j经过3条边的路径数。同理，如果要求经过k步的路径数，我们只需要二分求出A^k即可。
我感觉这类题确实不错，把矩阵和图论巧妙地结合在了一起，我之前做的POJ2888(传送阵)就是要用的这方面的知识，收获很大。

经典题目9 用1 x 2的多米诺骨牌填满M x N的矩形有多少种方案，M<=5，N<2^31，输出答案mod p的结果



    我们以M=3为例进行讲解。假设我们把这个矩形横着放在电脑屏幕上，从右往左一列一列地进行填充。其中前n-2列已经填满了，第n-1列参差不齐。现在我们要做的事情是把第n-1列也填满，将状态转移到第n列上去。由于第n-1列的状态不一样（有8种不同的状态），因此我们需要分情况进行讨论。在图中，我把转移前8种不同的状态放在左边，转移后8种不同的状态放在右边，左边的某种状态可以转移到右边的某种状态就在它们之间连一根线。注意为了保证方案不重复，状态转移时我们不允许在第n-1列竖着放一个多米诺骨牌（例如左边第2种状态不能转移到右边第4种状态），否则这将与另一种转移前的状态重复。把这8种状态的转移关系画成一个有向图，那么问题就变成了这样：从状态111出发，恰好经过n步回到这个状态有多少种方案。比如，n=2时有3种方案，111->011->111、111->110->111和111->000->111，这与用多米诺骨牌覆盖3x2矩形的方案一一对应。这样这个题目就转化为了我们前面的例题8。

8和9能转化，而且Matrix67大神在blog里有写这题代码，我就不献丑了~

经典题目10 POJ2778

    题目大意是，检测所有可能的n位DNA串有多少个DNA串中不含有指定的病毒片段。合法的DNA只能由ACTG四个字符构成。题目将给出10个以内的病毒片段，每个片段长度不超过10。数据规模n<=2 000 000 000。

    下面的讲解中我们以ATC,AAA,GGC,CT这四个病毒片段为例，说明怎样像上面的题一样通过构图将问题转化为例题8。我们找出所有病毒片段的前缀，把n位DNA分为以下7类：以AT结尾、以AA结尾、以GG结尾、以?A结尾、以?G结尾、以?C结尾和以??结尾。其中问号表示“其它情况”，它可以是任一字母，只要这个字母不会让它所在的串成为某个病毒的前缀。显然，这些分类是全集的一个划分（交集为空，并集为全集）。现在，假如我们已经知道了长度为n-1的各类DNA中符合要求的DNA个数，我们需要求出长度为n时各类DNA的个数。我们可以根据各类型间的转移构造一个边上带权的有向图。例如，从AT不能转移到AA，从AT转移到??有4种方法（后面加任一字母），从?A转移到AA有1种方案（后面加个A），从?A转移到??有2种方案（后面加G或C），从GG到??有2种方案（后面加C将构成病毒片段，不合法，只能加A和T）等等。这个图的构造过程类似于用有限状态自动机做串匹配。然后，我们就把这个图转化成矩阵，让这个矩阵自乘n次即可。最后输出的是从??状态到所有其它状态的路径数总和。

这题不要看大神说得这么轻松，但对我这种不太明白AC自动机，只会敲模板的渣渣来说，写起来还是很有挑战性的。建议深刻理解题意。

附上：本人代码

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<string>
#include<cmath>
#include<set>
#include<vector>
#include<stack>
#define mem(a,b) memset(a,b,sizeof(a))
#define FOR(a,b,i) for(i=a;i<=b;++i)
#define For(a,b,i) for(i=a;i<b;++i)
#define N 100000
using namespace std;
inline void RD(int &ret)
{
char c;
do
{
c=getchar();
}
while(c<'0'||c>'9');
ret=c-'0';
while((c=getchar())>='0'&&c<='9')
{
ret=ret*10+(c-'0');
}
}
inline void OT(int a)
{
if(a>=10)
{
OT(a/10);
}
putchar(a%10+'0');
}
int n,k;
int trie[111][5],fail[111],top;
int h[128];
void init()
{
mem(trie[0],0);
top=1;
h['A']=0;
h['C']=1;
h['T']=2;
h['G']=3;
}
void maketrie(char *s)//建trie
{
int rt,nxt;
for(rt=0; *s; rt=nxt,s++)
{
if(!(nxt=trie[rt][h[*s]]))
{
mem(trie[top],0);
trie[rt][h[*s]]=nxt=top++;
}
}
trie[rt][4]=1;
}
void makeacauto()//ac自动机模板
{
int q
,bg,ed,u,v,i;
fail[0]=bg=ed=0;
FOR(0,3,i)
{
if(v=trie[0][i])
{
fail[q[ed++]=v]=0;
}
}
while(bg<ed)
{
u=q[bg++];
FOR(0,3,i)
{
if(v=trie[u][i])
{
q[ed++]=v;
fail[v]=trie[fail[u]][i];
trie[v][4]|=trie[fail[v]][4];
}
else
{
trie[u][i]|=trie[fail[u]][i];
}
}
}
}
__int64 ty[111][111],tz[111][111],tw[111][111],s;
void mat(__int64 a[111][111],__int64 b[111][111]) //矩阵乘
{
__int64 d[111][111];
int i,j,l;
mem(d,0);
For(0,top,i)
{
For(0,top,j)
{
For(0,top,l)
{
d[i][j]=(d[i][j]+a[i][l]*b[l][j])%N;
}
}
}
For(0,top,i)
{
For(0,top,j)
{
a[i][j]=d[i][j];
}
}
}
void g(int x)//矩阵快速幂，用的都是自己用的模板
{
int i,j;
mem(ty,0);
For(0,top,i)
{
ty[i][i]=1;
}
while(x>0)
{
if(x%2==1)
{
mat(ty,tw);
}
mat(tw,tw);
x/=2;
}
}
int main()
{
int i,j,m,n;
char ss[11];
RD(m);
RD(n);
init();
while(m--)
{
scanf("%s",ss);
maketrie(ss);
}
makeacauto();
mem(tw,0);
For(0,top,i)
{
FOR(0,3,j)
{
if(trie[i][4]==0)
{
tw[i][trie[i][j]]++;
}
}
}
g(n);
s=0;
For(0,top,i)
{
if(trie[i][4]==0)
{
s+=ty[0][i];
}
}
printf("%I64d\n",s%N);
return 0;
}

以上内容均是个人的矩阵小结，感觉矩阵太神奇了，潜力无穷~~~