BZOJ 2219: 数论之神

#include<cstring>
#include<iostream>
#include<cstdio>
#include<map>
#include<cmath>
#include<algorithm>
#define rep(i,l,r) for (int i=l;i<=r;i++)
#define down(i,l,r) for (int i=l;i>=r;i--)
#define clr(x,y) memset(x,y,sizeof(x))
#define maxn 50050
#define inf 2000000000
using namespace std;
int a,b,k,p,t,ans;
map<int,int> mp;
int read(){
    int x=0,f=1; char ch=getchar();
    while (!isdigit(ch)) {if (ch==‘-‘) f=-1; ch=getchar();}
    while (isdigit(ch)) {x=x*10+ch-‘0‘; ch=getchar();}
    return x*f;
}
int Pow(int x,int y){
    int ans=1;
    while (y){
        if (y&1) ans=1LL*ans*x;
        x=1LL*x*x; y>>=1;
    }
    return ans;
}
int Pow(int x,int y,int mo){
    int ans=1; while (y){
        if (y&1) ans=1LL*ans*x%mo;
        x=1LL*x*x%mo; y>>=1;
    }
    return ans;
}
int divs[maxn],dtot;
int groot(int p){//是对于p-1的每一个素因子a去检查(p-1)/a！！ 
    dtot=0;
    int q=p-1;
    for (int i=2;i*i<=q;i++){
        if (q%i==0) {
            divs[++dtot]=i;
            if (i*i!=q) divs[++dtot]=q/i;
        }
    }
    for (int i=2;i<p;i++){
        int ok=0;
        rep(j,1,dtot) if (Pow(i,divs[j],p)==1) {ok=1;break;}
        if (!ok) return i;
    }
}
int exgcd(int a,int b,int &x,int &y){
    if (b==0) {
        x=1; y=0; return a;
    }
    int d=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return d;
}
int gcd(int a,int b){
    if (b==0) return a;
    return gcd(b,a%b);
}
int NI(int a,int b){
    int x,y; 
    exgcd(a,b,x,y);
    return (x%b+b)%b;
}
int bsgs(int a,int b,int p){
    mp.clear();
    int m=sqrt(p)+1,mul,tmp;
    int now=1; rep(i,0,m-1) {if (!mp.count(now)) mp[now]=i; now=1LL*now*a%p;}
    now=1; mul=Pow(a,m,p);
    rep(i,0,m-1){
        tmp=1LL*b*NI(now,p)%p;
        if (mp.count(tmp)) return mp[tmp]+i*m;
        now=1LL*now*mul%p;
        if (now==0) break;
    } 
    return inf;
}
int solve(int a,int b,int p,int d){
    int pd=Pow(p,d);
    b=b%pd;
    if (b==0) return Pow(p,d-((d-1)/a+1));
    else {
        int bet=0;        
        while (b%p==0) bet++,b/=p;
        if (bet%a!=0) return 0;
        int tt=bet/a;
        int g=groot(p),phi=pd-pd/p;
        int ind=bsgs(g,b,pd);
        int d=gcd(a,phi);
        if (ind%d==0) return d*Pow(p,(a-1)*tt);
        else return 0;
    }    
}
int main(){
    t=read();
    while (t--){
        a=read(); b=read(); k=read();
        p=2*k+1;
        ans=1;
        for (int i=2;i*i<=p&&ans;i++) if ((p%i)==0) {
            int now=0; while ((p%i)==0) p/=i,now++;
            ans=ans*solve(a,b,i,now);
        }
        if (p!=1&&ans) ans=ans*solve(a,b,p,1); 
        printf("%d\n",ans);
    }
    return 0;
}