BZOJ2219 数论之神

题意为求$x^A=B\;(mod\;P)$解的个数。
把$P$分解质因数得$P=p_1^{q_1}*p_2^{q_2}*...*p_n^{q_n}$
之后列出n个方程：$x^A=B\;(mod\;p_i^{q_i})$
那么这n个方程解的个数的乘积即为原方程解的个数。
因为我们可以从每个解集中选一个解，根据中国剩余定理，它对应了一个原方程的解，所以原方程解的个数就是每个方程解的个数的乘积。
考虑求$x^A=B\;(mod\;p^a)$解的个数。
1.
如果$B=0$，那么$x$肯定是形如$k*p^t$的形式，其中$t$为最小的$A*t>=a$的$t$，即$\lfloor{\frac{a-1}A}\rfloor+1$，解就有$p^{a-t}$个。
2.
如果$gcd(B,p^a)=1$，我们需要引入几个概念。
【阶】
阶$\delta(a)$表示最小的$x>0$使得$a^x=1\;(mod\;p)$
【原根】
如果$\delta(g)=\varphi(p)$，则$g$是$p$的原根。
【指标】
指标$ind(a)$表示在模$p$意义下以$p$的原根$g$为底，$a$的对数，即$$g^{ind(a)}=a\;(mod\;p)$$【指标的性质】
$a=b\;(mod\;p)\Leftrightarrow ind(a)=ind(b)\;(mod\;\varphi(p))$
$ind(a*b)=ind(a)+ind(b)\;(mod\;\varphi(p))$
$ind(a^k)=k*ind(a)\;(mod\;\varphi(p))$

那么原方程相当于$$ind(x^A)=ind(B)\;(mod\;\varphi(p))$$$$A*ind(x)=ind(B)\;(mod\;\varphi(p))$$这是一个线性方程，设$G=gcd(A,\varphi(p))$，当$ind(B)\%G!=0$时无解，否则解的个数为$G$（每$\frac{\varphi(p)}G$有一个解，一共就有$\frac{\varphi(p)}{\frac{\varphi(p)}G}=G$个）
3.
如果$gcd(B,p^a)>1$，设$B$是形如$p^k*b$的形式，移项得$$p^{-k}*x^A=b\;(mod\;p^{a-k})$$那么如果$k\%A\;!=0$，则方程无解。
否则可以化成$$(x*p^{-\frac kA})^A=b\;(mod\;p^{a-k})$$此时$gcd(b,p^{a-k})=1$，可以转到情况2.
需要注意的是，此方程并不与原方程同解，因为这个方程的解的取值范围是$[0,p^{a-k+\frac kA})$，而原方程的解的取值范围为$[0,p^a)$，所以最后需要乘上$p^{k-\frac kA}$

#include <cstdio>
#include <map>
#include <cmath>
using namespace std;

typedef long long ll;
const int nf=0x3f3f3f3f;
int T,a,b,p,a1,d[1000000];
map<ll,int> mp;

int gcd(int a,int b) {return b?gcd(b,a%b):a;}
ll pw(ll a,int b,int p) {ll r=1; for(;b;b>>=1,a=a*a%p) if(b&1) r=r*a%p; return r;}

int bsgs(int a,int b,int p) {
    int m=ceil(sqrt(p));
    ll A=1,a1=1,A1=pw(a,m,p);
    mp.clear();
    for(int i=1;i<=m;i++) {
        A=A*a%p;
        if(!mp[A*b%p]) mp[A*b%p]=i;
    }
    for(int i=1;i<=m;i++) {
        a1=a1*A1%p;
        if(mp[a1]) return i*m-mp[a1];
    }
    return -1;
}

int sol(int A,int b,int p,int a,int pa) {
    b%=pa;
    if(!b) return pw(p,a-(a-1)/A-1,nf);
    if(gcd(b,pa)>1) {
        int tt=0;
        while(gcd(b,pa)^1) b/=p,tt++,pa/=p;
        if(tt%A) return 0;
        return sol(A,b,p,a-tt,pa)*pw(p,tt-tt/A,nf);
    }
    int fi=pa-pa/p,tt=0,g;
    for(int i=2;i*i<=fi;i++) if(fi%i==0) d[++tt]=i,d[++tt]=fi/i;
    for(g=2;;g++) {
        for(int i=1;i<=tt;i++) if(pw(g,d[i],pa)==1) goto nx;
        break;
        nx: ;
    }
    int B=bsgs(g,b,pa);
    if(!~B) return 0;
    if(B%gcd(A,fi)) return 0;
    return gcd(A,fi);
}

int main() {
    scanf("%d",&T);
    while(T--) {
        scanf("%d%d%d",&a,&b,&p),p=2*p+1,a1=1;
        for(int i=2;i*i<=p;i++) if(p%i==0) {
            int tt=0,pa=1;
            while(p%i==0) p/=i,tt++,pa*=i;
            a1*=sol(a,b,i,tt,pa);
        }
        if(p^1) a1*=sol(a,b,p,1,p);
        printf("%d\n",a1);
    }
    return 0;
}