HDU 3221 上海09 B题矩阵乘法数论

最新推荐文章于 2022-02-25 23:45:16 发布

原创最新推荐文章于 2022-02-25 23:45:16 发布 · 1k 阅读

0 ·

本内容遵循CC 4.0 BY-SA版权协议

标签

#matrix #class #c

本文介绍了一种利用矩阵快速幂解决模指数问题的方法，并通过一个具体的编程实现展示了如何高效地进行大数的模指数运算。该方法适用于密码学及算法竞赛等场景。

用到一个公式A^x%C=A^(x%phi(C)+phi(C))%C (x>=phi(C));

注意全部用__int64,用int有可能溢出

#include<iostream> using namespace std; #define MAXN 20 __int64 Modal; bool flag; __int64 eular(__int64 n)//取出1到n中跟n互质的数的个数 { __int64 ret = 1,i; for (i = 2;i * i <= n;i++) if (n % i == 0) { n /= i; ret *= (i - 1); while (n % i == 0) { n /= i; ret *= i; } } if (n > 1) ret *= (n - 1); return ret; } class Matrix { public: __int64 M[MAXN][MAXN]; __int64 n,m;//n*m的矩阵 friend Matrix operator *(const Matrix &M1,const Matrix &M2)//矩阵乘法MatrixMultiply(MM?- -!) { __int64 i,j,k; Matrix M3; memset(M3.M,0,sizeof(M3.M)); if (M1.m!=M2.n) cout<<"These two Matrix can't be Multiplied!"<<endl; for(i=1;i<=M1.n;i++) for(j=1;j<=M2.m;j++) for(k=1;k<=M1.m;k++) { if (M3.M[i][j]+M1.M[i][k]*M2.M[k][j]>=Modal) { flag=true; } M3.M[i][j]=(M3.M[i][j]+M1.M[i][k]*M2.M[k][j]%Modal)%Modal;//对Modal求模 } M3.n=M1.n;M3.m=M2.m; return M3; } }; Matrix QuickMM(Matrix M,__int64 n) { if (n==1) return M; if (n==2) return M*M; if (n%2==0) return QuickMM(QuickMM(M,n/2),2); else return QuickMM(QuickMM(M,n/2),2)*M; } Matrix M1,M2,ans; int main() { __int64 n,p,a,b,x,v,t_case,i; memset(M2.M,0,sizeof(M2.M)); M2.n=2;M2.m=2; M2.M[1][1]=0;M2.M[1][2]=1; M2.M[2][1]=1;M2.M[2][2]=1; memset(M1.M,0,sizeof(M1.M)); M1.n=1;M1.m=2; M1.M[1][1]=1;M1.M[1][2]=0; __int64 t; scanf("%I64d",&t); t_case=0; while(++t_case<=t) { scanf("%I64d %I64d %I64d %I64d",&a,&b,&p,&n); Modal=eular(p); M1.M[1][1]=1;M1.M[1][2]=0; flag=false; if (n>=3) { ans=M1*QuickMM(M2,n-2); x=ans.M[1][2]; } else x=M1.M[1][n]; if (flag) x+=Modal; v=1; for(i=1;i<=x;i++) v=v*a%p; flag=false; M1.M[1][1]=0;M1.M[1][2]=1; if (n>=3) { ans=M1*QuickMM(M2,n-2); x=ans.M[1][2]; } else x=M1.M[1][n]; if (flag) x+=Modal; for(i=1;i<=x;i++) v=v*b%p; printf("Case #%I64d: %I64d/n",t_case,v); } }