bzoj2219: 数论之神

在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)

输出一行，表示答案

3
213 46290770 80175784
3 46290770 80175784
3333 46290770 80175784

27
27
297

新加数组一组--2015.02.27

数论鸣谢 AekdyCoin

不愧是数论之神。。。

设2*k+1=p1^a1*p2^a2*...*pn^an（不会打数学符号。。。求神犇在评论中告诉我），ans(%（2*k+1）)=ans(%p1^a1)*ans(p2^a2)*...*ans(pn^an)。

所以把2*k+1拆分后分情况计算就好了。

写了好久。。。又把bsgs里两个变量弄混了，定义了两个y。。。讲道理为什么不会报ce啊。。。

  1 #include<cstdio>
  2 #include<cstring>
  3 #include<iostream>
  4 #include<cmath>
  5 using namespace std;
  6 int T,A,B,k;
  7 long long ans;
  8 long long power(int x,int y){
  9     if(y==0)return 1;
 10     if(y==1)return x;
 11     long long tmp=power(x,y/2);tmp=tmp*tmp;
 12     if(y&1)tmp=tmp*x;
 13     return tmp;
 14 }
 15 long long work1(int A,int B,int p,int cnt){return power(p,cnt-((cnt-1)/A)-1);}
 16 #define base 1000000
 17 int h[base+5],pre[50005],val[50005];
 18 void insert(int x,int y){val[x]=y;y%=base;pre[x]=h[y];h[y]=x;}
 19 int get(int x){
 20     int y=-1;
 21     for(int i=h[x%base];i!=-1;i=pre[i]){
 22         if(val[i]==x&&(y==-1||i<y))y=i;
 23     }
 24     return y;
 25 }
 26 int gcd(int a,int b){return b==0?a:gcd(b,a%b);}
 27 int exgcd(int a,int b,int &x,int &y){
 28     if(b==0){x=1;y=0;return a;}
 29     int tt=exgcd(b,a%b,x,y),xx=x;
 30     x=y;y=xx-a/b*y;return tt;
 31 }
 32 int bsgs(int a,int b,int p){
 33     int m=ceil(sqrt(p*1.0)),D=1,fuckpps=1;
 34     memset(h,-1,sizeof(h));
 35     for(int i=0;i<m;i++)insert(i,fuckpps),fuckpps=1LL*fuckpps*a%p;
 36     for(int i=0,j=-1;i<m;i++){
 37         int x,y,d=exgcd(D,p,x,y);
 38         x=(1LL*x*(b/d)%(p/d)+(p/d))%(p/d);j=get(x);
 39         if(j!=-1)return i*m+j;
 40         D=1LL*D*fuckpps%p;
 41     }
 42     return 0;
 43 }
 44 long long power(int x,int y,int p){
 45     if(y==0)return 1;
 46     if(y==1)return x;
 47     long long tmp=power(x,y/2,p);tmp=tmp*tmp%p;
 48     if(y&1)tmp=tmp*x%p;
 49     return tmp;
 50 }
 51 int tot,f[1000005];
 52 int findg(int pp){
 53     int pps=pp-1;tot=0;
 54     for(int i=2;i*i<=pp-1;i++){
 55         if(pps%i==0){
 56             f[++tot]=i;
 57             while(pps%i==0)pps/=i;
 58         }
 59     }
 60     if(pps!=1)f[++tot]=pps;
 61     for(int i=1;;i++){
 62         int j=1;
 63         for(;j<=tot;j++){
 64             if(power(i,(pp-1)/f[j],pp)==1)break;
 65         }
 66         if(j==tot+1)return i;
 67     }
 68     return 0;
 69 }
 70 long long work2(int A,int B,int p,int cnt){
 71     int pp=power(p,cnt),g=findg(p);
 72     int indB=bsgs(g,B%pp,pp);
 73     int a=A,b=pp-pp/p,c=indB,d=gcd(a,b);
 74     if(c%d)return 0;
 75     return d;
 76 }
 77 long long work3(int A,int B,int p,int cnt){
 78     int a=0;
 79     while(B%p==0)B/=p,a++;
 80     if(a%A)return 0;
 81     return work2(A,B,p,cnt-a)*power(p,a-a/A);
 82 }
 83 int main(){
 84     scanf("%d",&T);
 85     while(T--){
 86         scanf("%d%d%d",&A,&B,&k);k=k*2+1;ans=1;
 87         for(int i=2;i*i<=k;i++){
 88             if(k%i==0){
 89                 int cnt=0,y=1;
 90                 while(k%i==0)k/=i,cnt++,y*=i;
 91                 int d=gcd(B,y);
 92                 if(d==y)ans*=work1(A,B,i,cnt);
 93                 else if(d==1)ans*=work2(A,B,i,cnt);
 94                 else ans*=work3(A,B,i,cnt);
 95             }
 96         }
 97         if(k!=1){
 98             int d=gcd(B,k);
 99             if(d==k)ans*=work1(A,B,k,1);
100             else if(d==1)ans*=work2(A,B,k,1);
101         }
102         printf("%lld
",ans);
103     }
104     return 0;
105 }

View Code

bzoj2219: 数论之神

2219: 数论之神

Description

Input

Output

Sample Input

Sample Output

HINT

Source