【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的原根。就可以很快找出来啦。
解释一下情况1和情况2为什么范围扩大之后就直接乘：
例如情况1：t=(a-1)/A+1，[0,p^t]中有一个解，范围[0,p^a)中有p^(a-t)个这样的范围。
假设p^t就是解。那下一个小区间中的p^(2t)也是解，以此类推。
PS：找原根的方法：
1 预判n有没有原根，有原根的数为：1、2、4、P^a,2*P^a,P为任意奇素数
2
3 快速求所有原根：
4 m=phi(n)
5 找到m所有的质因子y
6 找出n最小的原根a：gcd(a,n)==1 && a^(m/y) %n都!=1
7 则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 1
128 if(B%pa==0) return quickpow(p,a-(((a-1)/A)+1),Inf);
129 //case 2
130 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 3
139 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 }
