「刷题」数论之神

#include<iostream>
#include<cstdio>
#include<cmath>
#include<map>
using  namespace std;
int T,a,b,k,num,ans,d,dp[10000000];
map<int,int> mp;
int qw(int a,int b,int mod)
{
    int ans=1;
    for(;b;b>>=1,a=1LL*a*a%mod) if(b&1) ans=1LL*ans*a%mod;
    return ans;
}
int gcd(int a,int b)
{
    return b==0?a:gcd(b,a%b);
}
int BSGS(int y,int z,int p)
{
    int n=sqrt(p)+1;
    z%=p;
    if(y==0) return z==0?0:-1;
    mp.clear();
    for(int i=0,t=z;i<n;i++,t=1LL*t*y%p) mp[t]=i;
    int fac=qw(y,n,p);
    for(int i=0,t=1;i<=n;i++,t=1LL*t*fac%p)
    {
        int rank=mp.find(t)==mp.end()?-1:mp[t];
        if(rank>=0&&i*n-rank>=0) return i*n-rank;
    }
    return -1;
}
int primeroot(int p,int phi)
{
    d=0;
    for(int i=2;1LL*i*i<=phi;i++)
        if(phi%i==0)
            dp[++d]=i,dp[++d]=phi/i;
    for(int i=2;;i++)
    {
        int j;
        for(j=1;j<=d;j++)
            if(qw(i,dp[j],p)==1)
                break;
        if(j==d+1) return i;
    }
    return 0;
}
int solit(int a,int b,int p,int k,int pk)
{
    int phi=pk-pk/p;
    int g=primeroot(pk,phi);
    int indB=BSGS(g,b,pk);
    int ans=gcd(phi,a);
    if(indB%ans!=0) return 0;
    return ans;
}
int stcit(int a,int b,int p,int k,int pk)
{
    if(b%pk==0) return qw(p,k-(k-1)/a-1,2100000000);
    b=b%pk;
    if(gcd(pk,b)>1)
    {
        int GCD=gcd(b,pk),cnt=0;
        while(GCD%p==0) GCD/=p,b/=p,pk/=p,cnt++;
        if(a%cnt!=0) return 0;
        return qw(p,-cnt/a+cnt,2100000000)*solit(a,b,p,k,pk);
    }
    return solit(a,b,p,k,pk);
}
void work(int a,int b,int k)
{
    for(int i=2;1LL*i*i<=k;i++)
        if(k%i==0)
        {
            int p=i,pk=1,num=0;
            for(;k%p==0;k/=p,pk*=p,num++);
            ans=ans*stcit(a,b,p,num,pk);
        }
    if(k>1) ans=ans*stcit(a,b,k,1,k);
}
int main()
{
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d%d%d",&a,&b,&k);
        k=k<<1|1;ans=1;
        work(a,b,k);
        printf("%d\n",ans);
    }
}