关于矩阵 和矩阵快速幂
矩阵相乘,就是线性代数的知识,给两个矩阵,,求出他们相乘之后的矩阵,公式就是这样
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 using namespace std; 5 const int maxn=100; 6 int c[maxn][maxn]; 7 int n,m,q; 8 int main(){ 9 int a[maxn][maxn]; 10 int b[maxn][maxn]; 11 scanf("%d%d%d",&n,&m,&q); 12 for(int i=1;i<=n;i++){ 13 for(int j=1;j<=m;j++){ 14 scanf("%d",&a[i][j]); 15 } 16 } 17 for(int i=1;i<=m;i++){ 18 for(int j=1;j<=q;j++){ 19 scanf("%d",&b[i][j]); 20 } 21 } 22 memset(c,0,sizeof(c)); 23 for(int i=1;i<=n;i++){ 24 for(int k=1;k<=m;k++){ 25 if(a[i][k]==0){ 26 continue; 27 } 28 for(int j=1;j<=q;j++){ 29 c[i][j]+=a[i][k]*b[k][j]; 30 } 31 } 32 } 33 for(int i=1;i<=n;i++){ 34 for(int j=1;j<=q;j++){ 35 printf("%d ",c[i][j]); 36 37 } 38 printf("\n"); 39 } 40 }
讲到矩阵快速幂,就先讲快速幂,
首先,快速幂的目的就是做到快速求幂,假设我们要求a^b,按照朴素算法就是把a连乘b次,这样一来时间复杂度是O(b)也即是O(n)级别,快速幂能做到O(logn)
原理:
求a^b,那么其实b是可以拆成二进制的,该二进制数第i位的权为2^(i-1),例如当b==11时,a^11=a^(2^0+2^1+2^3)
11的二进制是1011,11 = 2³×1 + 2²×0 + 2¹×1 + 2º×1,因此,我们将a¹¹转化为算 a^(2^0)*a^(2^1)*a^(2^3) ,看出来快的多了吧原来算11次,现在算三次,由于是二进制,很自然地想到用位运算这个强大的工具: & 和 >> ,&运算通常用于二进制取位操作,例如一个数 & 1 的结果就是取二进制的最末位。还可以判断奇偶x&1==0为偶,x&1==1为奇。>>运算比较单纯,二进制去掉最后一位。base就是逐步累乘,刚好符合,2^0,2^1
2^(i-1)次方,另外要注意的地方就是指数函数是爆炸性增长的,所以要注意范围,
1 int poww(int a,int b){ 2 int ans=1,base=a; 3 while(b!=0){ 4 if(b&1!=0) 5 ans*=base; 6 base*=base; 7 b>>=1; 8 } 9 return ans; 10 }
另外有的题目会要求对较大的数取mod的,所以这个快速幂也是要取模的
1 int quick(int a,int b,int c) 2 { 3 int ans=1; //记录结果 4 a=a%c; //预处理,使得a处于c的数据范围之下 5 while(b!=0) 6 { 7 if(b&1) ans=(ans*a)%c; //如果b的二进制位不是0,那么我们的结果是要参与运算的 8 b>>=1; //二进制的移位操作,相当于每次除以2,用二进制看,就是我们不断的遍历b的二进制位 9 a=(a*a)%c; //不断的加倍 10 } 11 return ans; 12 }
矩阵快速幂,思路是一样的,
这里贴下poj 3070 的代码
题目链接:http://poj.org/problem?id=3070
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 #include<cmath> 6 #include<stack> 7 #include<map> 8 #include<vector> 9 #include<queue> 10 #include<set> 11 #include<iomanip> 12 #include<cctype> 13 using namespace std; 14 const int MAXN=2e5+5; 15 const int INF=1<<30; 16 //const long long mod=1e9+7; 17 #define ll long long 18 #define edl putchar('\n') 19 #define sscc ios::sync_with_stdio(false);cin.tie(0);cout.tie(0); 20 #define FOR(i,a,b) for(int i=a;i<=b;i++) 21 #define ROF(i,a,b) for(int i=a;i>=b;i--) 22 #define FORLL(i,a,b) for(ll i=a;i<=b;i++) 23 #define ROFLL(i,a,b) for(ll i=a;i>=b;i--) 24 #define mst(a) memset(a,0,ssizeof(a)) 25 #define mstn(a,n) memset(a,n,ssizeof(a)) 26 #define zero(x)(((x)>0?(x):-(x))<eps) 27 ll mod=10000; 28 const int ssize=2; 29 struct Matrix { 30 ll a[ssize][ssize]; 31 Matrix() { 32 memset(a,0,sizeof(a)); 33 } 34 void init() { 35 for(int i=0; i<ssize; i++) 36 for(int j=0; j<ssize; j++) 37 a[i][j]=(i==j); 38 } 39 Matrix operator + (const Matrix &B)const { 40 Matrix C; 41 for(int i=0; i<ssize; i++) 42 for(int j=0; j<ssize; j++) 43 C.a[i][j]=(a[i][j]+B.a[i][j])%mod; 44 return C; 45 } 46 Matrix operator * (const Matrix &B)const { 47 Matrix C; 48 for(int i=0; i<ssize; i++) 49 for(int k=0; k<ssize; k++) 50 for(int j=0; j<ssize; j++) 51 C.a[i][j]=(C.a[i][j]+1LL*a[i][k]*B.a[k][j])%mod; 52 return C; 53 } 54 Matrix operator ^ (const ll &t)const { 55 Matrix A=(*this),res; 56 res.init(); 57 ll p=t; 58 while(p) { 59 if(p&1)res=res*A; 60 A=A*A; 61 p>>=1; 62 } 63 return res; 64 } 65 }; 66 67 int main() { 68 Matrix a,b; 69 int n; 70 while(scanf("%d",&n)==1&&n!=-1) { 71 a.a[0][0]=1; 72 a.a[0][1]=1; 73 a.a[1][0]=1; 74 a.a[1][1]=0; 75 76 b.a[0][0]=1; 77 b.a[1][0]=0; 78 if(n<=0){ 79 cout<<0<<endl; 80 } 81 else{ 82 a=a^(n-1); 83 b=a*b; 84 cout<<b.a[0][0]<<endl; 85 } 86 87 } 88 }