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 }
View Code

题目链接: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),用扩展欧几里得算法即可实现。

posted @ 2016-08-15 16:23  oyzx~  阅读(245)  评论(0编辑  收藏  举报