首页 > 代码库 > 【bzoj2219-数论之神】求解x^a==b(%n)-crt推论-原根-指标-BSGS
【bzoj2219-数论之神】求解x^a==b(%n)-crt推论-原根-指标-BSGS
http://www.lydsy.com/JudgeOnline/problem.php?id=2219
弄了一个晚上加一个午休再加下午一个钟。。终于ac。。TAT
数论渣渣求轻虐！！
题意：求解 x^A=B(mod n) 在0~n内解的个数。其中1 <= A, B <= 10^9, 1 <= K <= 5 * 10^8 （n=2*K+1）
首先先说这一题的弱化版：bzoj1319 http://www.lydsy.com/JudgeOnline/problem.php?id=1319
bzoj1319这题是保证了P为质数。
找到p的一个原根g,因为g^x构成一个缩系,g^x可以表示0~p-1所有数。
g^j=a(%p), g(%p)=1, (g^i)^k=1(%p)
g^ik=g^j (%p)
ik=j(%phi(p))
用BSGS求出j，exgcd求出所有i，x就是g^i。
分析这一题：P不一定是质数。
取模数不是质数，无法利用通常的方式解方程；
但是有中国剩余定理这个东西，定理的推论告诉我们：
一个取模数互质的同余方程组(未必线性)，组合起来之后，这个同余方程解的个数为各方程解的个数的乘积；
(组合起来的方程的取模数为所有数的积；实际上这里解的范围都是属于[0 ,自己取模数) )
这点十分重要呢，它不仅证明了解的求法，而且如果有任意一个方程无解，那么整个就都是无解的；
————引用自http://blog.csdn.net/ww140142/article/details/47814003
把n分解质因数。
接下来我们只需要处理方程x^A==B（%p^a）
再次引用题解。。只有第三种情况是我自己搞的。。
引用自大牛http://blog.csdn.net/regina8023/article/details/44863519
但是第三种情况我没看懂它怎么搞。。
这个时候就可以用bzoj1319的解法了！
找到p^a的一个原根g,因为g^x构成一个缩系,g^x可以表示0~p^a-1所有数。
有一个推论（我也不知道为什么）g是p的原根，则g是p^a的原根。就可以很快找出来啦。
PS：找原根的方法：
1 预判n有没有原根，有原根的数为：1、2、4、P^a,2*P^a,P为任意奇素数2 3 快速求所有原根：4 m=phi(n)5 找到m所有的质因子y6 找出n最小的原根a：gcd(a,n)==1 && a^(m/y) %n都!=17 则a^x%n(1<=x<m,gcd(x,m)==1) 是n所有原根。依照上题，化成
g^ik=g^j (%p^a)
ik=j(%phi(p^a))
用BSGS求出j，解i的个数就是答案。
这里又有一个可爱的推论。。我还是不知道为什么。。
ax-py=gcd(a,p)
解的个数为gcd(a,p)。
然后这题就做完了。
1 #include<cstdio> 2 #include<cstdlib> 3 #include<cstring> 4 #include<cmath> 5 #include<iostream> 6 #include<algorithm> 7 using namespace std; 8 9 typedef long long LL; 10 const LL N=100100,Inf=(LL)1e12; 11 struct node{LL d,id;}num[N]; 12 LL nl,fl,pxl,px[N],r[N],f[N]; 13 14 void find_px(LL n) 15 { 16 pxl=0; 17 for(LL i=2;i*i<=n;i++) 18 { 19 if(n%i==0) px[++pxl]=i,r[pxl]=0; 20 while(n>1 && n%i==0) n/=i,r[pxl]++; 21 if(n==1) break; 22 } 23 if(n>1) px[++pxl]=n,r[pxl]=1;//debug 24 } 25 26 LL gcd(LL a,LL b) 27 { 28 if(b==0) return a; 29 return gcd(b,a%b); 30 } 31 32 LL quickpow(LL x,LL y,LL n) 33 { 34 LL ans=1%n; 35 while(y) 36 { 37 if(y&1) ans=(ans*x)%n; 38 x=(x*x)%n;y>>=1; 39 } 40 return ans; 41 } 42 43 bool cmp(node x,node y){ 44 if(x.d!=y.d) return x.d<y.d; 45 return x.id<y.id; 46 } 47 48 LL find_j(LL t) 49 { 50 LL l=0,r=nl,mid; 51 while(l<=r) 52 { 53 mid=(l+r)>>1; 54 if(num[mid].d==t) return num[mid].id; 55 if(num[mid].d<t) l=mid+1; 56 if(num[mid].d>t) r=mid-1; 57 } 58 return -1; 59 } 60 61 LL BSGS(LL a,LL b,LL n,LL phi)//a^x==b(%n) 62 { 63 LL m,x,am,now,t; 64 m=(LL)ceil(sqrt((double)n)); 65 x=1%n; 66 nl=0;num[++nl].d=1,num[nl].id=0; 67 for(int i=1;i<=m;i++) 68 { 69 x=(x*a)%n; 70 num[++nl].d=x;num[nl].id=i; 71 } 72 am=x; 73 sort(num+1,num+1+nl,cmp); 74 now=1; 75 for(int i=2;i<=nl;i++) 76 { 77 if(num[i].d!=num[i-1].d) num[++now]=num[i]; 78 } 79 nl=now; 80 am=quickpow(am,phi-1,n); 81 t=b%n; 82 for(int i=0;i<=m;i++) 83 { 84 x=find_j(t); 85 if(x!=-1) return i*m+x; 86 t=(t*am)%n; 87 } 88 return -1; 89 } 90 91 LL find_root(LL p) 92 { 93 LL x=p-1; 94 fl=0; 95 for(int i=2;i*i<=p-1;i++) 96 { 97 if((p-1)%i==0) f[++fl]=i,f[++fl]=(p-1)/i;//debug不是找质因子啊。。 98 } 99 for(int i=2;i<p-1;i++)100 {101 bool bk=1;102 for(int j=1;j<=fl;j++)103 if(quickpow(i,(p-1)/f[j],p)==1) {bk=0;break;}104 if(bk) return i;105 }106 }107 108 LL solve_3(LL A,LL B,LL p,LL a)109 {110 LL phi,g,gc,j,pa;111 pa=quickpow(p,a,Inf);112 phi=(p-1)*quickpow(p,a-1,Inf);113 g=find_root(p);114 j=BSGS(g,B,pa,phi);115 gc=gcd(A,phi);116 // printf("phi = %lld j = %lld g = %lld pa = %lld\n",phi,j,g,pa);117 // printf("s3 %lld %lld %lld %lld = %lld\n\n",A,B,p,a,gc);118 if(j%gc) return 0;119 return gc;120 }121 122 LL solve(LL A,LL B,LL p,LL a)123 {124 LL g,pa,x,y,b,cnt;125 pa=quickpow(p,a,Inf);126 g=gcd(pa,B);127 //case 1128 if(B%pa==0) return quickpow(p,a-(((a-1)/A)+1),Inf);129 //case 2130 if(g>1) 131 {132 b=B/g;133 cnt=0;x=g;134 while(x%p==0) x/=p,cnt++;135 if(cnt%A) return 0;136 return solve_3(A,b,p,a-cnt)*quickpow(p,cnt-(cnt/A),Inf);137 }138 //case 3139 return solve_3(A,B,p,a);140 }141 142 int main()143 {144 freopen("a.in","r",stdin);145 // freopen("me.out","w",stdout);146 int T;147 scanf("%d",&T);148 LL A,B,n,ans;149 while(T--)150 {151 scanf("%lld%lld%lld",&A,&B,&n);152 n=2*n+1;153 find_px(n);154 ans=1;155 for(LL i=1;i<=pxl;i++)156 {157 ans*=solve(A,B,px[i],r[i]);158 if(ans==0) break;159 }160 printf("%lld\n",ans);161 }162 return 0;163 }
【bzoj2219-数论之神】求解x^a==b(%n)-crt推论-原根-指标-BSGS