BZOJ 2219 数论之神 （BSGS + CRT + 原根）

 

 

简单粗暴的题意：求 X^A ≡ B(mod 2K+1) 在 X∈[0,2K]内有几个解。

首先，这个2*K+1并不一定是一个质数，因此我们可以考虑对2*K+1进行分解质因数。根据中国剩余定理，如果我们可以把2*K+1拆成，那么我们可以把这一个方程拆成n个方程，每个方程是一样的同余方程，但是模数分别对应这n个质数的幂。根据定理，如果这n个方程中有一个方程无解，那么整个方程无解。若有解，那么解的个数是各个方程的解的个数的乘积。那么我们的问题就变成了求每一个 X^A ≡ B(mod p^a) 在 X∈[0,2K]的解的个数，其中p为质数。

针对这个式子，我们考虑分情况讨论。

 

我们考虑当(B,p^a)=B时，有：

                                                         

我们设X=c*p^t，那么有：

                                                    

则有A*t>=a，那么，也即只要满足t的条件，c可以任取，对应就有个解。

 

当(B,p^a)=1时，两个互质且p^a有原根g。因为互质所以B在p^a下有逆元，对应的可以BSGS求出满足条件的y使得g^y=B(mod p^a)。同理，也可以用BSGS求出另外一边X对应的指数x，这样就有：

                                                  

当能够求出x的时候，解的个数就是gcd(phi(p^a),A)，否则就是0。

 

当(B,p^a)>1的时候，B=p^cnt*b。可以写成：

                                                  

我们把式子转换一下，可以有：

                                                  

当cnt%A!=0的时候，方程无解。而当等于0的时候，就变成了第二种情况，用求原根用BSGS解即可。但是这种情况注意到做左边的在式子里体现的取值比实际取值范围小，所以最后数量上要再乘多一个p^(cnt-cnt/A)。

关于代码，这里BSGS的时候，之前是用欧拉也即快速幂来求逆元的，但是这里由于不一定是指数，所以要么改成费马小，用phi或者用拓展欧几里德来求逆元。具体见代码：

#include <bits/stdc++.h>
#define INF 1e18
#define LL long long
#define sc(x) scanf("%d",&x)
#define scc(x,y) scanf("%lld%lld",&x,&y)
#define sccc(x,y,z) scanf("%lld%lld%lld",&x,&y,&z)
using namespace std;

const int N = 100010;

LL pri[N],t[N],p[N],tot;
map<LL, LL> mp;
LL K,n,m,ans;

LL qpow(LL x,LL n,LL p)
{
    LL res=1;
    while(n)
    {
        if (n&1) res=res*x%p;
        x=x*x%p; n>>=1;
    }
    return res;
}

LL ex_gcd(LL a,LL b,LL &x,LL &y)
{
    if(b==0){x=1; y=0;return a;}
    else
    {
        LL r=ex_gcd(b,a%b,y,x);
        y-=x*(a/b); return r;
    }
}

LL BSGS(LL A,LL B,LL C)
{
    int m=ceil(sqrt(C+0.5));
    mp.clear();
    LL now=1;
    for (int i=1;i<=m;i++)
    {
        now=now*A%C;
        if (!mp[now]) mp[now]=i;
    }
    mp[1]=0;
    A=qpow(A,m,C);
    now=1LL;
    for (int i=0;i<=m;i++)
    {
        LL x,y;
        ex_gcd(now,C,x,y);
        x=(x*B%C+C)%C;
        if (mp.count(x)) return i*m+mp[x];
        now=now*A%C;
    }
    return 0;
}

LL GetPrimitiveRoot(LL x,LL phi)
{
    int tt=0;
    for(int i=2;(LL)i*i<=phi;i++)
        if (phi%i==0) p[++tt]=i,p[++tt]=phi/i;
    for (int g=2;;g++)
    {
        int j;
        for (j=1;j<=tt;j++)
            if (qpow(g,p[j],x)==1) break;
        if (j==tt+1) return g;
    }
}

LL solve(LL a, LL b, LL P, LL pri)
{
    LL phi=P-P/pri;
    LL g=GetPrimitiveRoot(P,phi);
    LL x=BSGS(g,b%P,P);
    LL gcd=__gcd(a,phi);
    if (x%gcd) return 0;
    return gcd;
}

int main()
{
    int T; sc(T);
    while(T--)
    {
        tot=0;ans=1;
        sccc(n,m,K);
        K=K<<1|1;
        for(int i=2;(LL)i*i<=K;i++)
        {
            if (K%i) continue;
            pri[++tot]=i; t[tot]=0;
            while(K%i==0) K/=i,t[tot]++;
        }
        if (K!=1) pri[++tot]=K,t[tot]=1;
        for(int i=1;i<=tot;i++)
        {
            if (ans==0) break;
            LL P=qpow(pri[i],t[i],INF);
            if (m%P==0) ans=ans*qpow(pri[i],t[i]-(t[i]-1)/n-1,INF);
            else
            {
                LL tmp=m,cnt=0;
                while(tmp%pri[i]==0)
                    cnt++,tmp/=pri[i],t[i]--;
                P=qpow(pri[i],t[i],INF);
                if (cnt%n) {ans=0;break;}
                ans=ans*solve(n,tmp,P,pri[i])*qpow(pri[i],cnt-cnt/n,INF);
            }
        }
        printf("%lld\n",ans);
    }
    return 0;
}