PH算法
Pohlig Hellman算法：
1 #include<bits/stdc++.h> 2 using namespace std; 3 #define ll long long 4 #define ld long double 5 const int pr[]={2,3,5,7,11,13,31}; 6 int tot;ll fac[105],M[105],R[105];map<ll,ll>mp,h; 7 inline ll mul(ll a,ll b,ll p){return (a*b-(ll)((ld)a*b/p)*p+p)%p;} 8 inline ll ml(ll a,ll b,ll p){ll r=0;for(;b;b>>=1,a=(a+a)%p)if(b&1)r=(r+a)%p;return r;} 9 inline ll pw(ll a,ll b,ll p){ll r=1;for(;b;b>>=1,a=mul(a,a,p))if(b&1)r=mul(r,a,p);return r;} 10 ll gcd(ll a,ll b){return !b?a:gcd(b,a%b);} 11 ll exgcd(ll a,ll b,ll&x,ll&y){if(!b){x=1;y=0;return a;}int g=exgcd(b,a%b,y,x);y-=a/b*x;return g;} 12 ll crt(int n,ll*m,ll*r)//x%m_i=r_i 13 { 14 ll M=1,R=0,x,y; 15 for(int i=1;i<=n;i++)M*=m[i]; 16 for(int i=1;i<=n;i++) 17 { 18 ll t=M/m[i],d=exgcd(t,m[i],x,y); 19 R=(R+mul(mul(r[i],t,M),(x%m[i]+m[i])%m[i],M))%M; 20 } 21 return (R+M)%M; 22 } 23 inline bool Miller_Rabin(ll n) 24 { 25 if(n==2||n==3||n==5||n==7||n==11||n==13||n==31)return 1; 26 if(!(n&1))return 0;ll r=n-1; 27 int k=0;while((r&1)==0)r>>=1,k++; 28 for(int i=0;i<7;i++) 29 { 30 ll x=pw(pr[i],r,n),y=mul(x,x,n); 31 for(int i=0;i<k;i++,x=y,y=mul(x,x,n))if(y==1&&x!=1&&x!=n-1)return 0; 32 if(y!=1)return 0; 33 } 34 return 1; 35 } 36 #define Func(x) mul(x,x,n)+1 37 inline ll Pollard_Rho(ll n) 38 { 39 ll x=rand()%n+1,y=Func(x);ll k=1;int tst=0; 40 while(x!=y) 41 { 42 ll t=k;k=mul(k,abs(y-x),n); 43 if(!k){t=gcd(t,n);if(t>1)return t;break;} 44 if(tst==100){tst=0;ll t=gcd(k,n);if(t>1)return t;} 45 x=Func(x);y=Func(Func(y));tst++; 46 } 47 ll t=gcd(k,n);return t>1?t:-1; 48 } 49 void Factorization(ll n) 50 { 51 while(!(n&1)){fac[++tot]=2;mp[2]++;n/=2;} 52 if(n==1)return;else if(Miller_Rabin(n)){fac[++tot]=n;mp[n]++;return;} 53 ll p;do p=Pollard_Rho(n);while(p==-1);Factorization(p);Factorization(n/p); 54 } 55 ll Ord(ll x,ll a,ll p) 56 { 57 tot=0;mp.clear();Factorization(x); 58 for(int i=1;i<=tot;i++)if(!(x%fac[i])&&pw(a,x/fac[i],p)==1)x/=fac[i]; 59 tot=0;mp.clear();Factorization(x);return x; 60 } 61 ll qry(ll b,ll p,ll q){for(int i=0;i<q;i++)if(b==h[i])return i;return -1;} 62 ll Pohlig_Hellman(ll a,ll b,ll p) 63 { 64 ll n=Ord(p-1,a,p);int cc=0; 65 for(auto it:mp) 66 { 67 h.clear();ll q=it.first,e=it.second; 68 ll t=1,c=pw(a,n/q,p),bb=b,ans=0,x=n; 69 for(int i=0;i<q;i++,t=mul(t,c,p))h[i]=t; 70 for(ll i=1,qi=1;i<=e;i++,qi*=q) 71 { 72 x=n/qi/q;ll res=pw(bb,x,p),c; 73 c=qry(res,p,q);if(c==-1)return -1; 74 ll t=c*qi,iv=pw(pw(a,t,p),p-2,p); 75 ans+=t;bb=mul(bb,iv,p); 76 } 77 cc++;M[cc]=pw(q,e,p);R[cc]=ans; 78 } 79 return cc?crt(cc,M,R):-1; 80 } 81 int main() 82 { 83 int T;scanf("%d",&T); 84 while(T--) 85 { 86 ll p,a,b;scanf("%lld%lld%lld",&p,&a,&b); 87 printf("%lld\n",Pohlig_Hellman(a,b,p)); 88 } 89 return 0; 90 }
