【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所有的质因子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 }