Bzoj2219 数论之神
Time Limit: 3 Sec Memory Limit: 259 MB
Submit: 954 Solved: 268
Description
在ACM_DIY群中，有一位叫做“傻崽”的同学由于在数论方面造诣很高，被称为数轮之神！对于任何数论问题，他都能瞬间秒杀！一天他在群里面问了一个神题: 对于给定的3个非负整数 A,B,K 求出满足 (1) X^A = B(mod 2*K + 1) (2) X 在范围[0, 2K] 内的X的个数！自然数论之神是可以瞬间秒杀此题的，那么你呢？
Input
第一行有一个正整数T,表示接下来的数据的组数( T <= 1000) 之后对于每组数据，给出了3个整数A,B,K (1 <= A, B <= 10^9, 1 <= K <= 5 * 10^8)
Output
输出一行，表示答案
Sample Input
213 46290770 80175784
3 46290770 80175784
3333 46290770 80175784
Sample Output
27
297
HINT
新加数组一组--2015.02.27
Source
数学问题 原根 阶 指标 中国剩余定理 脑洞题
吼题
学姐的讲解很棒棒 http://blog.csdn.net/regina8023/article/details/44863519
模数$P=2*K+1$很大，在这个范围下没什么方法可以有效计算，所以需要优化数据范围。
将P分解质因数，$p_{1}^{a1}+p_{2}^{a2}+p_{3}^{a3}+...$
分别在每个模$ p_{i}^{ai} $ 的意义下计算出答案，将这些答案累乘起来就是最终的答案。（在每个模意义下选出一个解，则构成了一组同余方程，根据中国剩余定理，每组同余方程对应一个唯一解，所以可以用乘法原理统计总答案数）
假设当前我们在处理一个$p^a$
现在需要求$x^a \equiv b (\mod p^a)$的解个数
分三类情况讨论：
如果$ gcd(p^a,b)==b $，则x是p的某次幂的形式，直接出解
如果$ gcd(p^a,b)==1 $，问题转化为求$A*indx \equiv indb (\mod p^a)$的解的个数 (ind是指标)
根据扩展欧几里得的推论，若$indb \% gcd(A,phi(p^a)) == 0 $，解的个数为 gcd(A,phi(p^a))
如果$ gcd(p^a,b)>1 $，可以通过放缩范围转化成第二类情况。
求逆元时候要注意模数不一定是质数，所以要用欧拉定理来求。
没意识到这个，直接算p-2次幂，WA得飞起。
吐槽一下数据真的弱，我求原根的时候x%pri==0写成x%pri==1，居然过了样例和discuss的hack
隔壁yhx注释掉BSGS直接交发现A了，也就是说数据其实保证有解，根本不用判0……
1 #include<iostream> 2 #include<algorithm> 3 #include<cstdio> 4 #include<cmath> 5 #include<cstring> 6 #include<map> 7 #define LL long long 8 using namespace std; 9 const int mxn=100010; 10 int read(){ 11 int x=0,f=1;char ch=getchar(); 12 while(ch<'0' || ch>'9'){if(ch=='-')f=-1;ch=getchar();} 13 while(ch>='0' && ch<='9'){x=x*10+ch-'0';ch=getchar();} 14 return x*f; 15 } 16 LL ksm(LL a,LL k){ 17 LL res=1; 18 while(k){ 19 if(k&1)res=res*a; 20 a=a*a; 21 k>>=1; 22 } 23 return res; 24 } 25 LL ksm(LL a,LL k,int p){ 26 LL res=1;a%=p; 27 while(k){ 28 if(k&1)res=(LL)res*a%p; 29 a=(LL)a*a%p; 30 k>>=1; 31 } 32 return res; 33 } 34 int gcd(int a,int b){return (!b)?a:gcd(b,a%b);} 35 int exgcd(int a,int b,LL &x,LL &y){ 36 if(!b){ x=1;y=0;return a;} 37 int tmp=exgcd(b,a%b,x,y); 38 int c=x;x=y;y=c-a/b*x; 39 return tmp; 40 } 41 // 42 int P; 43 int pri[mxn],cnt=0; 44 int GetG(int p){//求原根 45 int tmp=p-1; 46 int m=sqrt(tmp+0.5); 47 cnt=0; 48 for(int i=2;i<=m;i++) 49 if(tmp%i==0){ 50 while(tmp%i==0)tmp/=i; 51 pri[++cnt]=i; 52 } 53 if(tmp>1)pri[++cnt]=tmp; 54 for(int g=2;g<p;g++){ 55 bool is=1; 56 for(int i=1;i<=cnt;i++) 57 if(ksm(g,(p-1)/pri[i],p)==1){is=0;break;} 58 if(is)return g; 59 } 60 return 0; 61 } 62 map<int,int>mp; 63 int BSGS(int a,int b,int p,int pp){ 64 mp.clear(); 65 int m=sqrt(p),i,j; 66 int tmp=1; 67 for(i=1;i<=m;i++){ 68 tmp=(LL)tmp*a%p; 69 if(tmp==b)return i; 70 if(!mp[tmp])mp[tmp]=i; 71 } 72 mp[1]=0; 73 int res=1,phi; 74 if(pp==p)phi=p-1; 75 else phi=p/pp*(pp-1); 76 for(i=0;i<=p;i+=m){ 77 LL inv=ksm(res,phi-1,p); 78 inv=inv*b%p; 79 if(mp[inv])return mp[inv]+i; 80 res=(LL)res*tmp%p; 81 } 82 return 0; 83 } 84 int G,A,B,K; 85 int solve(int p,int t){//p^t 86 int pr=ksm(p,t),b=B; 87 b%=pr; 88 if(!b)return ksm(p,t-((t-1)/A+1)); 89 int c=0; 90 while(b%p==0)b/=p,c++; 91 if(c%A)return 0;//无解 92 // 93 int tmp=c/A; 94 t-=c; 95 int phi=pr-pr/p; 96 G=GetG(p); 97 int ind=BSGS(G,b,pr,p); 98 int d=gcd(A,phi); 99 if(ind%d)return 0; 100 return d*ksm(p,(A-1)*tmp); 101 } 102 int main(){ 103 int i,j; 104 int T=read(); 105 while(T--){ 106 A=read();B=read();K=read(); 107 P=2*K+1;K=P; 108 int m=sqrt(P),c; 109 int ans=1; 110 for(i=2;i<=m;i++){ 111 if(!ans)break; 112 if(K%i==0){ 113 c=0; 114 while(K%i==0)K/=i,c++; 115 ans*=solve(i,c); 116 } 117 } 118 if(ans && K>1)ans*=solve(K,1); 119 printf("%d\n",ans); 120 } 121 return 0; 122 }
