今天开始学习矩阵方面的知识,主要参照大牛的博客十个利用矩阵乘法解决的经典题目
经典题目一:
给定n个点,m个操作,构造O(m+n)的算法输出m个操作后各点的位置。操作有平移、缩放、翻转和旋转
这里的操作是对所有点同时进行的。其中翻转是以坐标轴为对称轴进行翻转(两种情况),旋转则以原点为中心。如果对每个点分别进行模拟,那么m个操作总共耗时O(mn)。利用矩阵乘法可以在O(m)的时间里把所有操作合并为一个矩阵,然后每个点与该矩阵相乘即可直接得出最终该点的位置,总共耗时O(m+n)。假设初始时某个点的坐标为x和y,下面5个矩阵可以分别对其进行平移、旋转、翻转和旋转操作。预先把所有m个操作所对应的矩阵全部乘起来,再乘以(x,y,1),即可一步得出最终点的位置。
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#define pi acos(-1.0)
double p[1000005][3];
double y[3][3]={-1,0,0,0,1,0,0,0,1};
double x[3][3]={1,0,0,0,-1,0,0,0,1};
double k[3][3]={1,0,0,0,1,0,0,0,1};
double ans[3][3]={1,0,0,0,1,0,0,0,1};
void multi(double a[3][3])
{
int i,j,l;
double p;
for(i=0;i<3;i++)
for(j=0;j<3;j++)
k[i][j]=ans[i][j];
for(l=0;l<3;l++)
{
for(i=0;i<3;i++)
{
p=0;
for(j=0;j<3;j++)
{
p+=a[l][j]*k[j][i];
}
ans[l][i]=p;
}
}
}
int main()
{
int n,m,i,j,l,v;
char s[2];
double a,b,c;
scanf("%d%d",&n,&m);
for(i=0;i<n;i++)
{
scanf("%lf%lf",&p[i][0],&p[i][1]);
p[i][2]=1;
}
for(i=0;i<m;i++)
{
double w[3][3]={1,0,0,0,1,0,0,0,1};
getchar();
scanf("%s",s);
if(s[0]=='X')
{
multi(x);
}
if(s[0]=='Y')
{
multi(y);
}
if(s[0]=='M')
{
scanf("%lf%lf",&w[0][2],&w[1][2]);
multi(w);
}
if(s[0]=='S')
{
scanf("%lf",&a);
w[0][0]=a;
w[1][1]=a;
multi(w);
}
if(s[0]=='R')
{
scanf("%lf",&a);
a=(a/180.0)*pi;
w[0][0]=cos(a);
w[0][1]=-sin(a);
w[0][2]=0;
w[1][0]=sin(a);
w[1][1]=cos(a);
w[1][2]=0;
w[2][0]=0;
w[2][1]=0;
w[2][2]=1;
multi(w);
}
}
for(i=0;i<n;i++)
{
for(j=0;j<3;j++)
{c=0;
for(l=0;l<3;l++)
{ c+=ans[j][l]*p[i][l];}
if(j<2)
printf("%.1lf ",c);
}
printf("\n");
}
return 0;
}
经典题目二:
给定矩阵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的值即可。我们可以在计算过程中不断取模,避免高精度运算。
简单题,根据矩阵的乘法求fib的值,这一道题的传授的求fib的方法挺重要的
#include<stdio.h>
void fun(int a1[][2],int a2[][2])
{
int c[2][2],i,j,k;
for(i=0;i<2;i++)
for(j=0;j<2;j++)
{
c[i][j]=0;
for(k=0;k<2;k++)
c[i][j]=(c[i][j]+a1[i][k]*a2[k][j])%10000;
}
for(i=0;i<2;i++)
for(j=0;j<2;j++)
a1[i][j]=c[i][j];
}
int main()
{
int n;
while(~scanf("%d",&n))
{
int a[2][2]={{1,1},{1,0}};
int b[2][2]={{1,0},{0,1}};//初始化为单位矩阵,用来保存矩阵乘积的结果
if(n==-1) break;
while(n)
{
if(n&1) //n的值最后一定会变为1,最后一定会执行这一句将结果储存在b中
fun(b,a);
fun(a,a);
n>>=1;
}
printf("%d\n",b[1][0]);
}
}
hdu 1575 Tr A
和上一题思路差别不大,最后求的值为主对角线上值之和
#include<stdio.h>
#define M 9973
int n;
void fun(int a[][15],int b[][15])
{
int i,j,k,c[15][15];
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
c[i][j]=0;
for(k=0;k<n;k++)
c[i][j]=(c[i][j]+a[i][k]*b[k][j]%M)%M;
}
for(i=0;i<n;i++)
for(j=0;j<n;j++)
a[i][j]=c[i][j];
}
int main()
{
int a[15][15],b[15][15];//b保存矩阵乘积最后的结果
int t,k,i,j,count;
scanf("%d",&t);
while(t--)
{
scanf("%d%d",&n,&k);
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
scanf("%d",&a[i][j]);
if(i==j)
b[i][j]=1;
else
b[i][j]=0;
}
while(k)
{
if(k&1)
fun(b,a);
fun(a,a);
k>>=1;
}
count=0;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
if(i==j)
{
count+=b[i][j];
count%=M;
}
printf("%d\n",count);
}
return 0;
}
经典题目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,即可得到原问题的答案。
S=A + A^2 + A^3 + ... + A^m
1.求A^m
m为偶数:m&1==0,m=2k,A^2k=A^k * A^k
m为奇数:m&1==1,m=2k+1 A^(2k+1)=A^k * A^k *A
2.求S
m为偶数:m&1==0,m=2k,S=A+A^2....A^k + A^k(A+A^2....A^k );
m为奇数:m&1==1,m=2k+1,S=A+A^2....A^k +A^(k+1)+ A^(k+1)(A+A^2....A^k )
#include<stdio.h>
#include<string.h>
struct M
{
int p[35][35];
};
int n,m,k;
M add(M a,M b) //求两个矩阵的和
{
M res;
int i,j;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
res.p[i][j]=(a.p[i][j]+b.p[i][j])%m;
return res;
}
M mult(M a,M b) //求两个矩阵的乘积
{
M c;
int i,j,k;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
c.p[i][j]=0;
for(k=0;k<n;k++)
c.p[i][j]=(c.p[i][j]+a.p[i][k]*b.p[k][j]%m)%m;
}
return c;
}
M pow(M a,int k) //求矩阵a的k次方
{
M b;
int i;
memset(b.p,0,sizeof(b.p));
for(i=0;i<n;i++)
b.p[i][i]=1;
while(k)
{
if(k&1)
b=mult(b,a);
a=mult(a,a);
k>>=1;
}
return b;
}
M sum(M b,int k) //求b^1+……+b^k
{
M res,a=b;
if(k==1)
{
res=a;
return res;
}
else
{
M c=sum(a,k/2);
if(k&1)
{
M t=pow(a,k/2+1);
return add(add(c,t),mult(t,c));
}
else
{
M t=pow(a,k/2);
return add(c,mult(c,t));
}
}
}
int main()
{
M a;
int i,j;
while(~scanf("%d%d%d",&n,&k,&m))
{
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
scanf("%d",&a.p[i][j]);
a.p[i][j]%=m;
}
M b=sum(a,k);
for(i=0;i<n;i++)
{
for(j=0;j<n-1;j++)
printf("%d ",b.p[i][j]);
printf("%d\n",b.p[i][j]);
}
}
return 0;
}
求下标为g(x)=k*i+b,(i=0,1,…,n-1)的Fibonacci数列的和,比上一个有难度,难在怎么转化为上一题的模型
题解:sum=A^(k*0+b) + A^(k*1+b) + … + A^(k*(n-1)+b)
=A^b*[A^(k*0) + A^(k*1) + … + A^(k*(n-1))]
=A^b*[(A^k)^0 + (A^k)^1 + … + (A^k)^(n-1)]
我们可以二分分别求出A^b 和 A^k。再把A^k看成整体,中括号里
的等比数列再进行二分求解。
#include<stdio.h>
#include<string.h>
int m;
struct M
{
__int64 p[2][2];
};
M add(M a,M b)
{
M c;
int i,j;
for(i=0;i<2;i++)
for(j=0;j<2;j++)
c.p[i][j]=(a.p[i][j]+b.p[i][j])%m;
return c;
}
M mult(M a,M b)
{
int i,j,k;
M c;
for(i=0;i<2;i++)
for(j=0;j<2;j++)
{
c.p[i][j]=0;
for(k=0;k<2;k++)
c.p[i][j]=(c.p[i][j]+a.p[i][k]*b.p[k][j]%m)%m;
}
return c;
}
M pow(M a,int k)
{
M b={1,0,0,1};
while(k)
{
if(k&1)
b=mult(b,a);
a=mult(a,a);
k>>=1;
}
return b;
}
M sum(M a,int k)
{
M res;
if(k==1)
{
res=a;
return res;
}
else
{
M b=sum(a,k/2);
if(k&1)
{
M c=pow(a,k/2+1);
return add(add(b,c),mult(c,b));
}
else
{
M c=pow(a,k/2);
return add(b,mult(b,c));
}
}
}
int main()
{
int k,b,n,i,j;
while(~scanf("%d%d%d%d",&k,&b,&n,&m))
{
M a,s,ans,q={1,1,1,0},x={1,0,0,1};
s=pow(q,k);
a=pow(q,b);
ans=add(x,sum(s,n-1));
ans=mult(ans,a);
printf("%I64d\n",ans.p[1][0]);
}
return 0;
}
hdu 3306 Another kind of Fibonacci
有一定的难度,构造举证的方式比较难想到,如下
#include<stdio.h>
#include<string.h>
const int N=10007;
struct M
{
int p[4][4];
};
M multi(M a,M b)
{
int i,j,k;
M c;
for(i=0;i<4;i++)
for(j=0;j<4;j++)
{
c.p[i][j]=0;
for(k=0;k<4;k++)
{
c.p[i][j]=(c.p[i][j]+a.p[i][k]*b.p[k][j]%N)%N;
}
}
return c;
}
M pow(M a,int k)
{
M b;
memset(b.p,0,sizeof(b.p));
for(int i=0;i<4;i++)
b.p[i][0]=1;
while(k)
{
if(k&1)
b=multi(a,b);
a=multi(a,a);
k>>=1;
}
return b;
}
int main()
{
int n,x,y;
while(~scanf("%d%d%d",&n,&x,&y))
{
x%=N;
y%=N;
M a;
memset(a.p,0,sizeof(a.p));
a.p[0][0]=a.p[0][1]=a.p[2][1]=1;
a.p[3][1]=x;a.p[3][3]=y;
a.p[1][1]=x*x%N;
a.p[1][2]=y*y%N;
a.p[1][3]=2*x*y%N;
M b=pow(a,n);
printf("%d\n",b.p[0][0]%N);
}
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次方,再乘以初始序列即可。做出来了别忙着高兴,得意之时就是你灭亡之日,别忘了最后可能还有几个置换需要模拟。
简单的矩阵运用题,就是把一个字符串按照一定的操作,执行一定的次数后,输出字符串
这题目描述有点错误,不要被误导了
这个地方
hello" -> "elhol" -> "lhelo" -> "helol".
正确的方式是hello->elhlo——>helol->lhelo->elhol
#include<stdio.h>
#include<string.h>
int n;
struct M
{
int p[85][85];
};
M mult(M a,M b)
{
int i,j,k;
M c;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
c.p[i][j]=0;
for(k=0;k<n;k++)
c.p[i][j]=(c.p[i][j]+a.p[i][k]*b.p[k][j]);
}
return c;
}
M pow(M a,int k)
{
M b;
int i,j;
memset(b.p,0,sizeof(b.p));
for(i=0;i<n;i++)
b.p[i][i]=1;
while(k)
{
if(k&1)
b=mult(b,a);
a=mult(a,a);
k>>=1;
}
return b;
}
int main()
{
int i,j,k,l,m;
char s[85];
while(~scanf("%d%d",&n,&m),n+m)
{
M a,b;
memset(a.p,0,sizeof(a.p));
for(i=0;i<n;i++)
{
scanf("%d",&k);
a.p[k-1][i]=1;
}
getchar();
gets(s);
b=pow(a,m);
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
if(b.p[i][j])
printf("%c",s[j]);
}
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。例如,我们可以用下面的矩阵乘法来二分计算f(n) = 4f(n-1) - 3f(n-2) + 2f(n-4)的第k项:
利用矩阵乘法求解线性递推关系的题目我能编出一卡车来。这里给出的例题是系数全为1的情况。
hdu 1575 A Simple Math Problem
矩阵求线性递推方程
构造一个10*10的矩阵A,运用矩阵二分快速幂求得矩阵B=A^(k-9),用矩阵(f9,f8,f7
* ......f2,f1,f0)乘以B,求出的B[0][0]%m就是f(k)%m的值
#include<stdio.h>
#include<string.h>
int m;
struct M
{
__int64 p[10][10];
};
M mult(M a,M b)
{
M c;
memset(c.p,0,sizeof(c.p));
int i,j,k;
for(i=0;i<10;i++)
for(j=0;j<10;j++)
for(k=0;k<10;k++)
c.p[i][j]=(c.p[i][j]+(a.p[i][k]*b.p[k][j])%m)%m;
return c;
}
M pow(M a,int k)
{
M b;
int i;
memset(b.p,0,sizeof(b.p));
for(i=0;i<10;i++)
b.p[i][i]=1;
while(k)
{
if(k&1)
b=mult(b,a);
a=mult(a,a);
k>>=1;
}
return b;
}
int main()
{
__int64 k;
int aa[10],i,j,l;
while(~scanf("%I64d%d",&k,&m))
{
for(i=0;i<10;i++)
scanf("%d",&aa[i]);
if(k<10) {printf("%d\n",k);continue;}
M a;
memset(a.p,0,sizeof(a.p));
__int64 b[1][10]={9,8,7,6,5,4,3,2,1,0};
__int64 s[1][10];
memset(s,0,sizeof(s));
for(i=0;i<10;i++)
{
if(i<9)
a.p[i][i+1]=1;
a.p[i][0]=aa[i];
}
M c=pow(a,k-9);
for(i=0;i<1;i++)
for(j=0;j<10;j++)
for(l=0;l<10;l++)
s[i][j]=(s[i][j]+(b[i][l]*c.p[l][j])%m)%m;
printf("%I64d\n",s[0][0]%m);
}
return 0;
}