bzoj 2219: 数论之神 （BSGS+GCD+数论）

题目描述

传送门

题目大意：
求$X^A \equiv B (\space mod \space 2*K+1)$,在$[0,2k]$范围内解得个数。

题解

这道题用到的数论知识好多啊，还是一点点从头顺吧。

首先将$2*k+1$分解质因数，得到质因子个同余方程

$$X^A\equiv B (\space mod \space p_i^{d_i})$$
根据中国剩余定理的推论，模 $2*k+1$意义下的解等于模 $p_i^{d_i}$意义下解的乘积。

然后我们对于$X^A\equiv B (\space mod \space p_i^{d_i})$ 中$A,B, p_i^{d_i}$中的关系进行分类讨论。
(1)$B\equiv 0 (\space mod \space p^d)$
$X^A$是$p^d$的倍数，那么X中一定含有因子p。设p的指数是a，那么$a*A \ge d$.
我们只要找到最小的a，那么$p^a$在$[0,p^d)$范围内的所有倍数都是合法解。

$$a=\lfloor {d-1\over A}\rfloor +1$$
所以解的个数就是 ${p^d\over p^a}=p^{d-a}$

(2)$gcd(B,p^d)=1$
所有奇素数以及他们的幂都有原根。
我们求出$p^d$的原根，那么式子其实可以转化成

$$(g^x)^A \equiv g^{ind(B)} (\space mod \space p^d),其中ind(B)表示的是g^{ind(B)}\equiv B (\space mod \space p^d)$$
那么如何求原根呢？
原根：如果g是P的原根，那么g的 $1...p-1$次幂 $mod \space p$的结果一定互不相同.
根据欧拉定理
$a^{\phi(p)}\equiv 1 (\space mod \space p) ,其中gcd(a,p)=1$
$a^x\equiv 1 (\space mod \space p)$,当最小的x等于 $\phi(p)$时，a是p的原根。
那么我们是否需要将 $2..\phi(p)$中的数都判断一下呢，实际上只需要判断 $\phi(p)$的因子即可。如果因子中没有满足上式的，那么a就是p的原根。
一般原根的大小不会太大，所有我直接暴力枚举判断即可。
求出原根后，我们就得到了方程
$$Ax\equiv inv(B) (\space mod \space \phi(p^d))$$
上面方程的解的个数就是原式的解的个数。
形如 $$Ax\equiv C (\space mod \space m)$$ 的方程，根据扩展欧几里得的推论可知，
令 $t=gcd(A,m)$，如果t|C，那么方程解的个数为t，否则无解。
如何证明？扩展欧几里得定理，要求 $gcd(A,m)=1$，如果 $gcd(A,m)=g,g\neq1$,那么需要同时除去g，所以要求 $g|C$.
设 $x1,x2$为方程的两个解。
那么 $$A*x1\equiv A*x2 (\space mod \space m)$$
$$A(x1-x2) \equiv 0 (\space mod \space m)$$
$(m/g)| (A/g)*(x1-x2)$，因为 $m/g,A/g$互质，所以 $(m/g) |(x1-x2)$
那么x其实可以表示成 $x=x0+k*(m/g),其中k\in[0,g-1]$，所以解的个数为g。

(3)$gcd(B,p^d)\neq 1$
虽然B是$p^d$的倍数，但是B可能是p的倍数。
令$B=p^a*q,其中gcd(p,q)=1$
那么方程可以写成

$$X^A\equiv p^a*q (\space mod \space p^d)$$
其中A|a,否则无解。因为如果A不是a的因数，那么 $X^A$中约数p的个数不可能和B中的相等。
所以可以设 $a=k*A$,那么x实际上可以写成 $p^k*y$的形式。
那么方程就变成了
$$p^a*y^A\equiv p^a*q (\space mod \space p^d)$$
那么方程解的数量是否等于方程 $y^A\equiv q (mod p^{d-a})$的数量呢？
其实是不对的，因为我们在化简的时候定义域的范围改变了，x的取值范围是 $[0,p^d)$，因为 $x=y*p^a$,所以y的取值范围是 $[0,p^{d-k})$,而在第二个方程中y的取值范围变成了 $[0,p^{d-a})$，所以原始方程的解的数量应该是 $y^A\equiv q (\space mod \space p^{d-a})$数量的 $p^{a-k}$倍。
第二个方程的求解方式同情况(2)

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<map>
#define LL long long
#define N 100003 
using namespace std;
map<LL,int> mp;
int pd[N],prime[N],n,f[N];
LL m[N],cnt[N];
void init(int n)
{
    for (int i=2;i<=n;i++) {
        if (!pd[i]) prime[++prime[0]]=i;
        for (int j=1;j<=prime[0];j++){
            if (prime[j]*i>n) break;
            pd[prime[j]*i]=1;
            if (i%prime[j]==0) break;
        }
    }
}
LL gcd(LL x,LL y)
{
    LL r;
    while (y) {
        r=x%y;
        x=y; y=r;
    }
    return x;
}
LL get_phi(LL x)
{
    LL ans=x;
    for (int i=2;i*i<=x;i++)
     if (x%i==0) {
        ans=(LL)ans/i*(i-1);
        while (x%i==0) x/=i;
     }
    if (x>1) ans=ans/x*(x-1);
    return ans;
}
LL quickpow(LL num,int x)
{
    LL base=num; LL ans=1;
    while (x) {
        if (x&1) ans=ans*base;
        x>>=1;
        base=base*base;
    }
    return ans;
}
LL quickpow(LL num,int x,LL p){
    LL base=num%p; LL ans=1;
    while (x) {
        if (x&1) ans=ans*base%p;
        x>>=1;
        base=base*base%p;
    }
    return ans;
}
LL calc(LL p,LL phi)
{
    int c=0;
    for (int i=2;i*i<=phi;i++)
        if (phi%i==0)
            f[++c]=i,f[++c]=phi/i;
    for (int g=2;;g++)
    {
        int j;
        for (j=1;j<=c;j++)
            if (quickpow(g,f[j],p)==1) break;
        if (j==c+1) return g;
    }
    return 0;
}
LL bsgs(LL a,LL b,LL p)  
{  
    LL m=ceil(sqrt(p)); LL ans=b; LL d=1;
    mp.clear();  
    LL tmp=quickpow(a,m,p);  
    mp[ans]=0;  
    for (LL i=1;i<=m;i++) ans=ans*a%p,mp[ans]=i;  
    for (int LL i=1;i<=m+1;i++) {  
        d=d*tmp%p;  
        if (mp[d]) {  
            return (LL)i*m-mp[d];  
        }  
    }  
    return -1;  
}  
LL solve(LL A,LL B,LL P,LL k)
{
    LL p=quickpow(P,k);
    if (B%p==0) {
        LL a=(k-1)/A+1;
        return quickpow(P,k-a);
    }
    LL t=0; LL opt=1; LL X=B;
    while (B%P==0) B/=P,t++;
    if (t%A) return 0;
    if (gcd(X,p)!=1) {
        p=quickpow(P,k-t);
        opt=quickpow(P,t-t/A);
    }
    LL phi=get_phi(p); LL g=calc(p,phi);
    LL C=bsgs(g,B,p);
    if (C==-1) return 0;
    t=gcd(A,phi);
    if (C%t==0) return (LL)t*opt;
    else return 0;
}
int main()
{
    freopen("a.in","r",stdin);
    freopen("my.out","w",stdout);
    int T; scanf("%d",&T);
    init(100000);
    while (T--) {
       LL A,B,K;
       scanf("%lld%lld%lld",&A,&B,&K);
       LL p=2*K+1;
       n=0;
       for (int i=1;prime[i]*prime[i]<=p;i++)
        if (p%prime[i]==0) {
            m[++n]=prime[i]; cnt[n]=0;
            while (p%prime[i]==0) p/=prime[i],cnt[n]++;
        }
       if (p>1) m[++n]=p,cnt[n]=1;
       LL ans=1;
       for (int i=1;i<=n;i++)
        ans=ans*solve(A,B,m[i],cnt[i]);
       printf("%lld\n",ans);
    }
}