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

题目描述

传送门

题目大意：
求XA≡B( mod 2∗K+1),在[0,2k]范围内解得个数。

题解

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

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

XA≡B( mod pdii)

根据中国剩余定理的推论，模2∗k+1意义下的解等于模pdii意义下解的乘积。

然后我们对于XA≡B( mod pdii) 中A,B,pdii中的关系进行分类讨论。
(1)B≡0( mod pd)
XA是pd的倍数，那么X中一定含有因子p。设p的指数是a，那么a∗A≥d.
我们只要找到最小的a，那么pa在[0,pd)范围内的所有倍数都是合法解。

a=⌊d−1A⌋+1

所以解的个数就是pdpa=pd−a

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

(gx)A≡gind(B)( mod pd),其中ind(B)表示的是gind(B)≡B( mod pd)

那么如何求原根呢？
原根：如果g是P的原根，那么g的1...p−1次幂mod p的结果一定互不相同.
根据欧拉定理
aϕ(p)≡1( mod p),其中gcd(a,p)=1
ax≡1( mod p),当最小的x等于ϕ(p)时，a是p的原根。
那么我们是否需要将2..ϕ(p)中的数都判断一下呢，实际上只需要判断ϕ(p)的因子即可。如果因子中没有满足上式的，那么a就是p的原根。
一般原根的大小不会太大，所有我直接暴力枚举判断即可。
求出原根后，我们就得到了方程

Ax≡inv(B)( mod ϕ(pd))

上面方程的解的个数就是原式的解的个数。
形如

Ax≡C( mod m)

的方程，根据扩展欧几里得的推论可知，
令t=gcd(A,m)，如果t|C，那么方程解的个数为t，否则无解。
如何证明？扩展欧几里得定理，要求gcd(A,m)=1，如果gcd(A,m)=g,g≠1,那么需要同时除去g，所以要求g|C.
设x1,x2为方程的两个解。
那么

A∗x1≡A∗x2( mod m)

A(x1−x2)≡0( mod m)

(m/g)|(A/g)∗(x1−x2)，因为m/g,A/g互质，所以(m/g)|(x1−x2)
那么x其实可以表示成x=x0+k∗(m/g),其中k∈[0,g−1]，所以解的个数为g。

(3)gcd(B,pd)≠1
虽然B是pd的倍数，但是B可能是p的倍数。
令B=pa∗q,其中gcd(p,q)=1
那么方程可以写成

XA≡pa∗q( mod pd)

其中A|a,否则无解。因为如果A不是a的因数，那么XA中约数p的个数不可能和B中的相等。
所以可以设a=k∗A,那么x实际上可以写成pk∗y的形式。
那么方程就变成了

pa∗yA≡pa∗q( mod pd)

那么方程解的数量是否等于方程yA≡q(modpd−a)的数量呢？
其实是不对的，因为我们在化简的时候定义域的范围改变了，x的取值范围是[0,pd)，因为x=y∗pa,所以y的取值范围是[0,pd−k),而在第二个方程中y的取值范围变成了[0,pd−a)，所以原始方程的解的数量应该是yA≡q( mod pd−a)数量的pa−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);    }}