矩阵乘法

    矩阵乘法是一种高效的算法可以把一些一维递推优化到log( n ),还可以求路径方案等,所以更是是一种应用性极强的算法。矩阵,是线性代数中的基本概念之一。一个m×n的矩阵就是m×n个数排成m行n列的一个数阵。由于它把许多数据紧凑的集中到了一起,所以有时候可以简便地表示一些复杂的模型。矩阵乘法看起来很奇怪,但实际上非常有用,应用也十分广泛。

基本定义 

  它是这样定义的,只有当矩阵A的列数与矩阵B的行数相等时A×B才有意义。一个m×n的矩阵am,n)左乘一个n×p的矩阵bn,p),会得到一个m×p的矩阵cm,p),满足

  矩阵乘法满足结合率,但不满足交换率

  一般的矩乘要结合快速幂才有效果``

  -------------------------------------------------------------------------------------------------------------

  矩阵乘法在编程语言中的运用例题(十个利用矩阵乘法解决的经典题目 )

  -------------------------------------------------------------------------------------------------------------

经典题目(转载自Matrix67原创) 

这里我们不介绍其它有关矩阵的知识,只介绍矩阵乘法和相关性质。

  不要以为数学中的矩阵也是黑色屏幕上不断变化的绿色字符。在数学中,一个矩阵说穿了就是一个二维数组。一个n行m列的矩阵可以乘以一个m行p列的矩阵,得到的结果是一个n行p列的矩阵,其中的第i行第j列位置上的数等于前一个矩阵第i行上的m个数与后一个矩阵第j列上的m个数对应相乘后所有m个乘积的和。比如,下面的算式表示一个2行2列的矩阵乘以2行3列的矩阵,其结果是一个2行3列的矩阵。其中,结果的那个4等于2*2+0*1:

  下面的算式则是一个1 x 3的矩阵乘以3 x 2的矩阵,得到一个1 x 2的矩阵:

  矩阵乘法的两个重要性质:一,矩阵乘法不满足交换律;二,矩阵乘法满足结合律

经典题目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^nnA相乘)的结果,输出的每个数都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,即可得到原问题的答案。

k=7时,有:A+A^2+A^3+A^4+A^5+A^6+A^7 =A+A^2+A^3+A^3*(A+A^2+A^3+A^3*A^1)

注意取模。
代码
#include<stdio.h>
#include
<string.h>
int n,kk,m;
int ans[31][31],a[31][31],ori[31][31];
int tmp[31][31],tp[31][31];
void calSqu()
{
int i,j,k;
memset(tmp,
0,sizeof(tmp));
for(i=0;i<n;i++)
for(j=0;j<n;j++)
for(k=0;k<n;k++)
tmp[i][j]
+=(a[i][k]*a[k][j])%m;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
a[i][j]
=tmp[i][j]%m;
}
void cal(int t)
{
int i,j,k;
if(t==1)
return;
if(t%2==0)
{
cal(t
/2);
memset(tmp,
0,sizeof(tmp));
for(i=0;i<n;i++)
for(j=0;j<n;j++)
for(k=0;k<n;k++)
tmp[i][j]
+=(a[i][k]*ans[k][j])%m;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
ans[i][j]
=(ans[i][j]+tmp[i][j])%m;
calSqu();
}
else{
cal(t
/2);
memset(tp,
0,sizeof(tp));
for(i=0;i<n;i++)
for(j=0;j<n;j++)
tp[i][j]
=ans[i][j];
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
for(k=0;k<n;k++)
tp[i][j]
+=(a[i][k]*ori[k][j])%m;
tp[i][j]
%=m;
}
memset(tmp,
0,sizeof(tmp));
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
for(k=0;k<n;k++)
tmp[i][j]
+=(a[i][k]*tp[k][j])%m;
tmp[i][j]
%=m;
}
for(i=0;i<n;i++)
for(j=0;j<n;j++)
ans[i][j]
=(ans[i][j]+tmp[i][j])%m;
calSqu();
memset(tmp,
0,sizeof(tmp));
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
for(k=0;k<n;k++)
tmp[i][j]
+=(a[i][k]*ori[k][j])%m;
tmp[i][j]
%=m;
}
for(i=0;i<n;i++)
for(j=0;j<n;j++)
a[i][j]
=tmp[i][j];
}
}
int main()
{
int i,j;
scanf(
"%d%d%d",&n,&kk,&m);
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
scanf(
"%d",&a[i][j]);
ans[i][j]
=ori[i][j]=a[i][j];
}
cal(kk);
for(i=0;i<n;i++)
{
for(j=0;j<n-1;j++)
printf(
"%d ",ans[i][j]%m);
printf(
"%d\n",ans[i][j]%m);
}
return 0;
}

经典题目4

      VOJ1049

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

  首先将这m个置换合并起来(算出这m个置换的乘积),然后接下来我们需要执行这个置换k/m次(取整,若有余数则剩下几步模拟即可)。注意任意一个置换都可以表示成矩阵的形式。例如,将1 2 3 4置换为3 1 2 4,相当于下面的矩阵乘法:

置换k/m次就相当于在前面乘以k/m个这样的矩阵。我们可以二分计算出该矩阵的k/m次方,再乘以初始序列即可。做出来了别忙着高兴,得意之时就是你灭亡之日,别忘了最后可能还有几个置换需要模拟。

 

经典题目5

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

  大家自己去看看吧,书上讲得很详细。解题方法和上一题类似,都是用矩阵来表示操作,然后二分求最终状态。

经典题目6

  给定np,求第nFibonaccimod p的值,n不超过2^31

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

 

经典题目7

  VOJ1067

我们可以用上面的方法二分求出任何一个线性递推式的第n项,其对应矩阵的构造方法为:在右上角的(n-1)*(n-1)的小矩阵中的主对角线上填1,矩阵第n行填对应的系数,其它地方都填0。例如,我们可以用下面的矩阵乘法来二分计算f(n) = 4f(n-1) - 3f(n-2) + 2f(n-4)的第k项:

利用矩阵乘法求解线性递推关系的题目我能编出一卡车来。

经典题目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列就表示从ij经过3条边的路径数。同理,如果要求经过k步的路径数,我们只需要二分求出A^k即可。

经典题目9

  1 x 2的多米诺骨牌填满M x N的矩形有多少种方案,M<=5N<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->111111->110->111111->000->111,这与用多米诺骨牌覆盖3x2矩形的方案一一对应。这样这个题目就转化为了我们前面的例题8

 
   此题的源代码(Matrix67原创)。你可以再次看到位运算的相关应用,熟悉了一下C++

的类和运算符重载。

代码
#include <cstdio>
#include
<iostream>
#define SIZE (1<<m)
#define MAX_SIZE 32
using namespace std;
class CMatrix
{
public:
long element[MAX_SIZE][MAX_SIZE];
void setSize(int);
void setModulo(int);
CMatrix
operator* (CMatrix);
CMatrix power(
int);
private:
int size;
long modulo;
};
void CMatrix::setSize(int a)
{
for (int i=0; i<a; i++)
for (int j=0; j<a; j++)
element[i][j]
=0;
size
= a;
}
void CMatrix::setModulo(int a)
{
modulo
= a;
}
CMatrix CMatrix::
operator* (CMatrix param)
{
CMatrix product;
product.setSize(size);
product.setModulo(modulo);
for (int i=0; i<size; i++)
for (int j=0; j<size; j++)
for (int k=0; k<size; k++)
{
product.element[i][j]
+=element[i][k]*param.element[k][j];
product.element[i][j]
%=modulo;
}
return product;
}
CMatrix CMatrix::power(
int exp)
{
CMatrix tmp
= (*this) * (*this);
if (exp==1) return *this;
else if (exp & 1) return tmp.power(exp/2) * (*this);
else return tmp.power(exp/2);
}
int main()
{
const int validSet[]={0,3,6,12,15,24,27,30};
long n, m, p;
CMatrix unit;
scanf(
"%d%d%d", &n, &m, &p);
unit.setSize(SIZE);
for(int i=0; i<SIZE; i++)
for(int j=0; j<SIZE; j++)
if( ((~i)&j) == ((~i)&(SIZE-1)) )
{
bool isValid=false;
for (int k=0; k<8; k++)isValid=isValid||(i&j)==validSet[k];
unit.element[i][j]
=isValid;
}
unit.setModulo(p);
printf(
"%d", unit.power(n).element[SIZE-1][SIZE-1] );
return 0;
}

经典题目10

  POJ2778

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

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

    题目中的数据规模保证前缀数不超过100,一次矩阵乘法是三方的,一共要乘log(n)次。因此这题总的复杂度是100^3 * log(n)AC了。

posted @ 2010-07-27 17:57  孟起  阅读(...)  评论(... 编辑 收藏