bzoj2219: 数论之神
1 #include <iostream> 2 #include <cstdio> 3 #include <cstring> 4 #include <cmath> 5 #include <algorithm> 6 using namespace std; 7 typedef long long int64; 8 int ca,top,prim[50005]; 9 int64 A,B,K,pi,pk,ans,ti,q; 10 int64 ksm(int64 x,int64 y){ 11 if (y==0) return 1; 12 if (y==1) return x; 13 int64 d=ksm(x,y/2); 14 if (y%2==1) return d*d*x; 15 else return d*d; 16 } 17 int64 exgcd(int64 a,int64 b,int64 &x,int64 &y){ 18 if (b==0){ 19 x=1,y=0; 20 return a; 21 } 22 int64 temp=exgcd(b,a%b,x,y),tmp; 23 tmp=x,x=y,y=tmp-a/b*y; 24 return temp; 25 } 26 int64 YuanGen(int64 x){ 27 bool flag; 28 top=0; 29 int64 temp=x-1; 30 for (int i=2;i<=sqrt(x-1);i++){ 31 if (temp%i==0){ 32 prim[++top]=i; 33 while (temp%i==0) temp/=i; 34 } 35 } 36 if (temp>1) prim[++top]=temp; 37 for (int i=1;;i++){ 38 flag=1; 39 for (int j=1;j<=top;j++){ 40 if (ksm(i,(x-1)/prim[j])%x==1){ 41 flag=0; 42 break; 43 } 44 } 45 if (flag) return i; 46 } 47 } 48 #define maxn 100005 49 #define maxm 400005 50 int now[maxn],prep[maxm]; 51 int64 val[maxm]; 52 void insert(int x,int64 y){ 53 int64 pos=y%maxn; 54 prep[x]=now[pos],now[pos]=x,val[x]=y; 55 } 56 int find(int64 x){ 57 int64 pos=x%maxn; int ans=maxm*4; 58 for (int i=now[pos];i!=-1;i=prep[i]){ 59 if (val[i]==x) ans=min(ans,i); 60 } 61 if (ans==maxm*4) return -1; 62 else return ans; 63 } 64 int64 BSGS(int64 A,int64 B,int64 C){ 65 memset(now,-1,sizeof(now)); 66 int64 pos,temp=ceil(sqrt(C*1.0)),D=1,R=1; 67 for (int i=0;i<temp;i++){ 68 insert(i,D); 69 D=D*A%C; 70 } 71 int64 tmp,x,y; 72 for (int i=0;i<temp;i++){ 73 tmp=exgcd(R,C,x,y); 74 x=(x*(B/tmp)%C+C)%C; 75 pos=find(x); 76 if (pos!=-1) return i*temp+pos; 77 R=R*D%C; 78 } 79 return -1; 80 } 81 int64 calc(int64 A,int64 B,int64 C){ 82 q=YuanGen(pi); 83 int64 a,b,t,x,y,c=C-C/pi; 84 b=BSGS(q,B,C); 85 t=exgcd(A,c,x,y); 86 c/=t; 87 x=(x*(b/t)%c+c)%c; 88 int sum=0; 89 for (int i=x;i<c*t;) sum++,i+=c; 90 return sum; 91 } 92 int64 work(int64 A,int64 B){ 93 B%=pk; int64 a,b,c,x; 94 if (B==0){ 95 a=(ti-1)/A+1; 96 return pk/ksm(pi,a); 97 } 98 b=0,c=1; 99 while (B%pi==0){ 100 b++,c*=pi; 101 B/=pi; 102 } 103 if (b%A!=0) return 0; 104 x=b/A; 105 return calc(A,B,pk/c)*ksm(pi,b-x); 106 } 107 int main(){ 108 int64 temp; 109 scanf("%d",&ca); 110 while (ca--){ 111 scanf("%lld%lld%lld",&A,&B,&K); 112 K=K*2+1,temp=K; ans=1; 113 for (int i=2;i<=sqrt(K);i++){ 114 if (temp%i==0){ 115 pi=i,pk=1,ti=0; 116 while (temp%i==0){ 117 pk*=i; ti++; 118 temp/=i; 119 } 120 ans*=work(A,B); 121 } 122 } 123 if (temp>1){ 124 pi=pk=temp; ti=1; 125 ans*=work(A,B); 126 } 127 printf("%lld\n",ans); 128 } 129 return 0; 130 }
题目链接：http://www.lydsy.com/JudgeOnline/problem.php?id=2219
题目大意：在ACM_DIY群中，有一位叫做“傻崽”的同学由于在数论方面造诣很高，被称为数轮之神！对于任何数论问题，他都能瞬间秒杀！一天他在群里面问了一个神题: 对于给定的3个非负整数 A,B,K 求出满足 (1) X^A = B(mod 2*K + 1) (2) X 在范围[0, 2K] 内的X的个数！自然数论之神是可以瞬间秒杀此题的，那么你呢？
第一行有一个正整数T,表示接下来的数据的组数( T <= 1000) 之后对于每组数据，给出了3个整数A,B,K (1 <= A, B <= 10^9, 1 <= K <= 5 * 10^8)
输出一行，表示答案。
做法：这题很像之前写过的一个题，用原根+指标+bsgs即可，但是那一题B，C互质，且C一定存在原根。这题不一样了，设C=2*k+1，C不一定存在原根，但是我们可以将其标准分解，这些pi^ti一定存在原根，因为C%2==1，一定不会有2的幂，而奇素数的a次幂，a>=1,一定存在原根。我们对每个x^A=B(mod pi^ti)（0=<x<pi^ti）求出解的个数后乘起来就是原来的解，至于为什么，右转百度搜中国剩余定理推论。
怎么求x^A=B(mod pi^ti)（0=<x<pi^ti）的解的个数呢？我们将B%pi^ti，为了一些边界情况，我们要分情况讨论，当B=0时，那么x^A一定是pi^ti的倍数，我们x中存在因子p^a,且a*A>=ti，那么只要是p^a的倍数就行，可以手推式子。
当B>0时，我们将B分解为pi^b * T，我们要想办法把pi^b约去，那么x^A中也要刚好有pi^b次方，否则之后同余不成立，所以当A不整除b时无解，令k=b/A,那么x=p^k * G,G属于[0,p^(ti-k)),约去后化简为G^A = T (mod pi^(ti-b)),G属于[0,pi^(ti-b)),mod数存在原根，T与mod数互质，转化为了弱化版，可以用bsgs+原根求解，设解数为ans，最后应该ans*=p^(b-k)，因为定义域缩小了，且存在大小为p^(ti-b)的循环节，除一下就知道最后要扩大的倍数了。
那么怎么求x^A=B(mod C)呢？ C存在原根，且（B，C）=1。我们找到C的原根g，找到B在原根g下modC的指标，x也必定为原根g的若干次方%C，式子变为
g^(aA)=g^b(mod C),由于循环节phi（C），变为aA=b（mod phi(C)）,0=<a<phi(C)，用扩展欧几里得算法即可实现。
