首页 > 代码库 > 矩阵快速幂
矩阵快速幂
poj 3070 Fibonacci http://poj.org/problem?id=3070
模板题，矩阵都给你写好了。
1 #include<cstdio> 2 #include<cstring> 3 #define mt(a,b) memset(a,b,sizeof(a)) 4 class Matrix { ///矩阵 5 typedef int typev;///权值类型 6 static const int MV=2;///矩阵长度 7 static const int mod=10000;///%mod 8 friend Matrix operator *(const Matrix &a,const Matrix &b) { ///矩阵乘法，必须a的列等于b的行 9 Matrix ret;10 ret.n=a.n;11 ret.m=b.m;12 ret.zero();13 if(a.m==b.n) {14 for(int k=0; k<a.m; k++) {15 for(int i=0; i<a.n; i++) {16 if(a.val[i][k]) {17 for(int j=0; j<b.m; j++) {18 ret.val[i][j]+=a.val[i][k]*b.val[k][j];19 ret.val[i][j]%=mod;///看具体需要20 }21 }22 }23 }24 }25 return ret;26 }27 friend Matrix operator ^ (Matrix &a,int b) {///必须是n*n方阵才能快速幂28 Matrix ret;29 ret.n=ret.m=a.n;30 ret.unit();31 for(; b; a=a*a,b>>=1) {32 if(b&1) {33 ret=ret*a;34 }35 }36 return ret;37 }38 public:39 int n,m;///n行m列40 typev val[MV][MV];41 void zero() {42 mt(val,0);43 }44 void unit() {45 zero();46 for(int i=0; i<MV; i++)47 val[i][i]=1;48 }49 } A;50 int main() {51 int n;52 while(~scanf("%d",&n),n>-1) {53 A.n=A.m=2;54 A.val[0][0]=1;55 A.val[0][1]=1;56 A.val[1][0]=1;57 A.val[1][1]=0;58 printf("%d\n",(A^n).val[0][1]);59 }60 return 0;61 }
hdu 1575 Tr A http://acm.hdu.edu.cn/showproblem.php?pid=1575
模板题
1 #include<cstdio> 2 #include<cstring> 3 #define mt(a,b) memset(a,b,sizeof(a)) 4 const int MOD=9973; 5 class Matrix { ///矩阵 下标0开始 6 typedef int typev;///权值类型 7 static const int MV=16;///矩阵长度 8 static const int mod=MOD;///%mod 9 friend Matrix operator *(const Matrix &a,const Matrix &b) { ///矩阵乘法，必须a的列等于b的行10 Matrix ret;11 ret.n=a.n;12 ret.m=b.m;13 ret.zero();14 if(a.m==b.n) {15 for(int k=0; k<a.m; k++) {16 for(int i=0; i<a.n; i++) {17 if(a.val[i][k]) {18 for(int j=0; j<b.m; j++) {19 ret.val[i][j]+=a.val[i][k]*b.val[k][j];20 ret.val[i][j]%=mod;///看具体需要21 }22 }23 }24 }25 }26 return ret;27 }28 friend Matrix operator ^ (Matrix &a,int b) {///必须是n*n方阵才能快速幂29 Matrix ret;30 ret.n=ret.m=a.n;31 ret.unit();32 for(; b; a=a*a,b>>=1) {33 if(b&1) {34 ret=ret*a;35 }36 }37 return ret;38 }39 public:40 int n,m;///n行m列41 typev val[MV][MV];42 void zero() {43 mt(val,0);44 }45 void unit() {46 zero();47 for(int i=0; i<MV; i++)48 val[i][i]=1;49 }50 } A;51 int main() {52 int t,n,k;53 while(~scanf("%d",&t)){54 while(t--){55 scanf("%d%d",&n,&k);56 A.n=A.m=n;57 for(int i=0;i<n;i++){58 for(int j=0;j<n;j++){59 scanf("%d",&A.val[i][j]);60 }61 }62 A=A^k;63 int ans=0;64 for(int i=0;i<n;i++){65 ans+=A.val[i][i];66 ans%=MOD;67 }68 printf("%d\n",ans);69 }70 }71 return 0;72 }
Kiki & Little Kiki 2 http://acm.hdu.edu.cn/showproblem.php?pid=2276
矩阵先不说了。
1 #include<cstdio> 2 #include<cstring> 3 #define mt(a,b) memset(a,b,sizeof(a)) 4 class Matrix { ///矩阵 下标0开始 5 typedef int typev;///权值类型 6 static const int MV=128;///矩阵长度 7 static const int mod=2;///%mod 8 friend Matrix operator *(const Matrix &a,const Matrix &b) { ///矩阵乘法，必须a的列等于b的行 9 Matrix ret;10 ret.n=a.n;11 ret.m=b.m;12 ret.zero();13 if(a.m==b.n) {14 for(int k=0; k<a.m; k++) {15 for(int i=0; i<a.n; i++) {16 if(a.val[i][k]) {17 for(int j=0; j<b.m; j++) {18 ret.val[i][j]+=a.val[i][k]*b.val[k][j];19 ret.val[i][j]%=mod;///看具体需要20 }21 }22 }23 }24 }25 return ret;26 }27 friend Matrix operator ^ (Matrix &a,int b) {///必须是n*n方阵才能快速幂28 Matrix ret;29 ret.n=ret.m=a.n;30 ret.unit();31 for(; b; a=a*a,b>>=1) {32 if(b&1) {33 ret=ret*a;34 }35 }36 return ret;37 }38 public:39 int n,m;///n行m列40 typev val[MV][MV];41 void zero() {42 mt(val,0);43 }44 void unit() {45 zero();46 for(int i=0; i<MV; i++)47 val[i][i]=1;48 }49 } A,T;50 int main(){51 int n;52 char s[128];53 while(~scanf("%d%s",&n,s)){54 int ls=strlen(s);55 A.n=A.m=T.n=T.m=ls;56 A.zero();57 T.unit();58 for(int i=0;i<ls;i++){59 if(s[i]==‘1‘){60 A.val[0][i]=1;61 }62 T.val[i][i+1]=1;63 }64 T.val[ls-1][0]=1;65 A=A*(T^n);66 for(int i=0;i<ls;i++){67 printf("%d",A.val[0][i]);68 }69 puts("");70 }71 return 0;72 }
Matrix Power Series http://poj.org/problem?id=3233
测模板，具体构造再说吧。
1 #include<cstdio> 2 #include<cstring> 3 #define mt(a,b) memset(a,b,sizeof(a)) 4 int mod; 5 class Matrix { ///矩阵 下标0开始 6 typedef int typev;///权值类型 7 static const int MV=64;///矩阵长度 8 friend Matrix operator *(const Matrix &a,const Matrix &b) { ///矩阵乘法，必须a的列等于b的行 9 Matrix ret;10 ret.n=a.n;11 ret.m=b.m;12 ret.zero();13 if(a.m==b.n) {14 for(int k=0; k<a.m; k++) {15 for(int i=0; i<a.n; i++) {16 if(a.val[i][k]) {17 for(int j=0; j<b.m; j++) {18 ret.val[i][j]+=a.val[i][k]*b.val[k][j];19 ret.val[i][j]%=mod;///看具体需要20 }21 }22 }23 }24 }25 return ret;26 }27 friend Matrix operator ^ (Matrix &a,int b) {///必须是n*n方阵才能快速幂28 Matrix ret;29 ret.n=ret.m=a.n;30 ret.unit();31 for(; b; a=a*a,b>>=1) {32 if(b&1) {33 ret=ret*a;34 }35 }36 return ret;37 }38 public:39 int n,m;///n行m列40 typev val[MV][MV];41 void zero() {42 mt(val,0);43 }44 void unit() {45 zero();46 for(int i=0; i<MV; i++)47 val[i][i]=1;48 }49 } A;50 int main(){51 int n,k;52 scanf("%d%d%d",&n,&k,&mod);53 A.n=A.m=n<<1;54 for(int i=0;i<n;i++){55 for(int j=0;j<n;j++){56 scanf("%d",&A.val[i][j]);57 A.val[i][j+n]=A.val[i][j];58 A.val[i+n][j]=0;59 A.val[i+n][j+n]=0;60 if(i+n==j+n)61 A.val[i+n][j+n]=1;62 }63 }64 A=A^k;65 for(int i=0;i<n;i++){66 for(int j=n;j<A.n;j++)67 printf("%d ",A.val[i][j]);68 puts("");69 }70 return 0;71 }
Training little cats http://poj.org/problem?id=3735
1 #include<cstdio> 2 #include<cstring> 3 #include<algorithm> 4 #define mt(a,b) memset(a,b,sizeof(a)) 5 using namespace std; 6 typedef __int64 LL; 7 class Matrix { ///矩阵 下标0开始 8 typedef LL typev;///权值类型 9 static const int MV=128;///矩阵长度10 static const int mod=10000;///%mod11 friend Matrix operator *(const Matrix &a,const Matrix &b) { ///矩阵乘法，必须a的列等于b的行12 Matrix ret;13 ret.n=a.n;14 ret.m=b.m;15 ret.zero();16 if(a.m==b.n) {17 for(int k=0; k<a.m; k++) {18 for(int i=0; i<a.n; i++) {19 if(a.val[i][k]) {20 for(int j=0; j<b.m; j++) {21 ret.val[i][j]+=a.val[i][k]*b.val[k][j];22 // ret.val[i][j]%=mod;///看具体需要23 }24 }25 }26 }27 }28 return ret;29 }30 friend Matrix operator ^ (Matrix &a,int b) {///必须是n*n方阵才能快速幂31 Matrix ret;32 ret.n=ret.m=a.n;33 ret.unit();34 for(; b; a=a*a,b>>=1) {35 if(b&1) {36 ret=ret*a;37 }38 }39 return ret;40 }41 public:42 int n,m;///n行m列43 typev val[MV][MV];44 void zero() {45 mt(val,0);46 }47 void unit() {48 zero();49 for(int i=0; i<MV; i++)50 val[i][i]=1;51 }52 } A;53 int main(){54 int x,y,n,m,k;55 char op[4];56 while(~scanf("%d%d%d",&n,&m,&k),n+m+k){57 A.unit();58 A.n=A.m=n+1;59 while(k--){60 scanf("%s%d",op,&x);61 switch(op[0]){62 case ‘g‘: A.val[0][x]++; break;63 case ‘e‘: for(int i=0;i<=n;i++) A.val[i][x]=0; break;64 default: scanf("%d",&y); for(int i=0;i<=n;i++) swap(A.val[i][x],A.val[i][y]);65 }66 }67 A=A^m;68 for(int i=1;i<=n;i++) printf("%I64d ",A.val[0][i]);69 puts("");70 }71 return 0;72 }
Fast Matrix Calculation http://acm.hdu.edu.cn/showproblem.php?pid=4965
1 #include<cstdio> 2 #include<cstring> 3 #define mt(a,b) memset(a,b,sizeof(a)) 4 typedef __int64 LL; 5 class Matrix { ///矩阵 下标0开始 6 typedef int typev;///权值类型 7 static const int MV=8;///矩阵长度 8 static const int mod=6;///%mod 9 friend Matrix operator *(const Matrix &a,const Matrix &b) { ///矩阵乘法，必须a的列等于b的行 10 Matrix ret; 11 ret.n=a.n; 12 ret.m=b.m; 13 ret.zero(); 14 if(a.m==b.n) { 15 for(int k=0; k<a.m; k++) { 16 for(int i=0; i<a.n; i++) { 17 if(a.val[i][k]) { 18 for(int j=0; j<b.m; j++) { 19 ret.val[i][j]+=a.val[i][k]*b.val[k][j]; 20 ret.val[i][j]%=mod;///看具体需要 21 } 22 } 23 } 24 } 25 } 26 return ret; 27 } 28 friend Matrix operator ^ (Matrix &a,int b) {///必须是n*n方阵才能快速幂 29 Matrix ret; 30 ret.n=ret.m=a.n; 31 ret.unit(); 32 for(; b; a=a*a,b>>=1) { 33 if(b&1) { 34 ret=ret*a; 35 } 36 } 37 return ret; 38 } 39 public: 40 int n,m;///n行m列 41 typev val[MV][MV]; 42 void zero() { 43 mt(val,0); 44 } 45 void unit() { 46 zero(); 47 for(int i=0; i<MV; i++) 48 val[i][i]=1; 49 } 50 } D; 51 struct mat { 52 int n,m; 53 int data[1024][1024]; 54 }A,B,C,E; 55 class MatrixOp { 56 public: 57 int mul(mat& c,const mat& a,const mat& b) { //c=a*b 58 int i,j,k; 59 if (a.m!=b.n) return 0; 60 c.n=a.n,c.m=b.m; 61 for (i=0; i<c.n; i++) 62 for (j=0; j<c.m; j++) 63 for (c.data[i][j]=k=0; k<a.m; k++) 64 c.data[i][j]+=a.data[i][k]*b.data[k][j]; 65 return 1; 66 } 67 } gx; 68 int main(){ 69 int N,K; 70 while(~scanf("%d%d",&N,&K),N|K){ 71 A.n=N; 72 A.m=K; 73 B.n=K; 74 B.m=N; 75 for(int i=0;i<N;i++){ 76 for(int j=0;j<K;j++){ 77 scanf("%d",&A.data[i][j]); 78 } 79 } 80 for(int i=0;i<K;i++){ 81 for(int j=0;j<N;j++){ 82 scanf("%d",&B.data[i][j]); 83 } 84 } 85 gx.mul(C,B,A); 86 D.n=D.m=C.n; 87 for(int i=0;i<C.n;i++){ 88 for(int j=0;j<C.m;j++){ 89 D.val[i][j]=C.data[i][j]; 90 } 91 } 92 D=D^(N*N-1); 93 for(int i=0;i<C.n;i++){ 94 for(int j=0;j<C.m;j++){ 95 C.data[i][j]=D.val[i][j]; 96 } 97 } 98 gx.mul(E,A,C); 99 gx.mul(C,E,B);100 LL ans=0;101 for(int i=0;i<C.n;i++){102 for(int j=0;j<C.m;j++){103 ans+=C.data[i][j]%6;104 }105 }106 printf("%I64d\n",ans);107 }108 return 0;109 }
ZOJ Problem Set - 2974 Just Pour the Water http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemId=1973
nenucontest3,一个序列经过相同变化，变化10^9次，明显是构造矩阵，然后快速幂。
1 #include<cstdio> 2 #include<cstring> 3 #define mt(a,b) memset(a,b,sizeof(a)) 4 class Matrix { ///矩阵 下标0开始 5 typedef double typev;///权值类型 6 static const int MV=128;///矩阵长度 7 // static const int mod=10000;///%mod 8 friend Matrix operator *(const Matrix &a,const Matrix &b) { ///矩阵乘法，必须a的列等于b的行 9 Matrix ret;10 ret.n=a.n;11 ret.m=b.m;12 ret.zero();13 if(a.m==b.n) {14 for(int k=0; k<a.m; k++) {15 for(int i=0; i<a.n; i++) {16 if(a.val[i][k]) {17 for(int j=0; j<b.m; j++) {18 ret.val[i][j]+=a.val[i][k]*b.val[k][j];19 // ret.val[i][j]%=mod;///看具体需要20 }21 }22 }23 }24 }25 return ret;26 }27 friend Matrix operator ^ (Matrix &a,int b) {///必须是n*n方阵才能快速幂28 Matrix ret;29 ret.n=ret.m=a.n;30 ret.unit();31 for(; b; a=a*a,b>>=1) {32 if(b&1) {33 ret=ret*a;34 }35 }36 return ret;37 }38 public:39 int n,m;///n行m列40 typev val[MV][MV];41 void zero() {42 mt(val,0);43 }44 void unit() {45 zero();46 for(int i=0; i<MV; i++)47 val[i][i]=1;48 }49 } A;50 double a[128];51 int main(){52 int t,n;53 scanf("%d",&t);54 while(t--){55 scanf("%d",&n);56 for(int i=0;i<n;i++){57 scanf("%lf",&a[i]);58 }59 A.unit();60 int k,id,m;61 for(int i=0;i<n;i++){62 scanf("%d",&k);63 if(k) A.val[i][i]=0;64 for(int j=0;j<k;j++){65 scanf("%d",&id);66 A.val[id-1][i]=1.0/k;67 }68 }69 A.n=A.m=n;70 scanf("%d",&m);71 A=A^m;72 for(int i=0;i<n;i++){73 if(i) printf(" ");74 double sum=0;75 for(int j=0;j<n;j++){76 sum+=A.val[i][j]*a[j];77 }78 printf("%.2f",sum);79 }80 puts("");81 }82 return 0;83 }
end
矩阵快速幂
声明：以上内容来自用户投稿及互联网公开渠道收集整理发布，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任，若内容有误或涉及侵权可进行投诉： 投诉/举报 工作人员会在5个工作日内联系你，一经查实，本站将立刻删除涉嫌侵权内容。