bzoj2219 数论之神 欧拉降幂+BSGS+原根

题目分析

题目要求$x^A\equiv B(\mathrm{m}\mathrm{o}\mathrm{d}P)$的解的个数。

首先将$P$分解质因数，对于每个方程$x^A\equiv b_i(\mathrm{m}\mathrm{o}\mathrm{d}{p}_i^{a_i})$，求出解的个数。假设我们确定了这个方程的解，最后用中国剩余定理合并，不同方程组肯定对应不同的解，所以将每个方程的解数乘起来就是最终解的个数。

对于单个方程，分三种情况讨论。

情况1:$b_i=0$。

说明若将$x$写成$p^{t}q$的形式，$t$因子的大小至少是$\lceil \frac{a}{A}\rceil$，那么$q$的大小就可以取$[1,p^{a-(\lceil \frac{a}{A}\rceil )}]$，并且其中一些取值有$p$因子所以也会包括$t$取更大的值的情况，所以该情况下的解数是$p^{a-(\lceil \frac{a}{A}\rceil )}$

情况2:$p|b$

设$b=p^{t}d$，则$x^A\equiv p^{t}d(\mathrm{m}\mathrm{o}\mathrm{d}{p}^a)$。若$A̸|t$，则无解，否则

$$(\frac{x}{p^{\frac{t}{A}}}{)}^A\equiv d(\mathrm{m}\mathrm{o}\mathrm{d}{p}^{a-t})$$

就转化成情况3了。

但是在情况2中，本来$\frac{x}{p^{\frac{t}{A}}}$的取值范围是$[0,p^{a-\frac{t}{A}})$的，所以解数还要乘一个$p^{t-\frac{t}{A}}$

情况3:$p̸|b$

求出$p^a$的原根$g$（$p^a$不是质数啊，剩余系大小是$\varphi (p^a)$啊QAQ）

那么让$g^r\equiv b(\mathrm{m}\mathrm{o}\mathrm{d}{p}^a)$，$r$用BSGS求出，则根据欧拉降幂得出下面的式子，原$x$的解数是下面这个式子的解数：

$Ax\equiv r(\mathrm{m}\mathrm{o}\mathrm{d}{p}^a\mathrm{m}\mathrm{o}\mathrm{d}\varphi (p^a))$

后面那一串记作$p^{\prime }$吧。

那么这个式子，若$gcd(p^{\prime },A)̸|r$，无解，否则根据同余方程通解的表现形式（解出的任意$x$，$x+t\frac{p^{\prime }}{gcd(p^{\prime },A)}$），解数为$gcd(p^{\prime },A)$。

代码

#include<bits/stdc++.h>
using namespace std;
#define RI register int
int T,A,B,P,kP,ans;
int pri[100005];

int gcd(int x,int y) {return y?gcd(y,x%y):x;}
int fast_pow(int x,int y) {
	int re=1;
	for(;y;y>>=1,x=x*x) if(y&1) re=re*x;
	return re;
}
int fast_pow(int x,int y,int p) {
	int re=1;
	for(;y;y>>=1,x=1LL*x*x%p) if(y&1) re=1LL*re*x%p;
	return re;
}
map<int,int> mp;
int bsgs(int x,int y,int p) {
	if(y==1) return 0;
	mp.clear();int m=ceil(sqrt((double)p)),q=1;
	for(RI i=0;i<m;++i) mp[(int)(1LL*q*y%p)]=i+1,q=1LL*q*x%p;
	for(RI i=m,x=1;;i+=m) {
		x=1LL*x*q%p;
		if(mp[x]) return i-mp[x]+1;
		if(i>p) return -1;
	}
}
int getg(int x,int phi) {
	int js=0,kphi=phi;
	for(RI i=2;i*i<=kphi;++i)
		if(kphi%i==0) {pri[++js]=i;while(kphi%i==0) kphi/=i;}
	if(kphi!=1) pri[++js]=kphi;
	for(RI i=2;;++i) {
		int flag=1;
		if(fast_pow(i,phi,x)!=1) continue;
		for(RI j=1;j<=js;++j)
			if(fast_pow(i,phi/pri[j],x)==1) {flag=0;break;}
		if(flag) return i;
	}
}
int getans(int p,int a,int kB) {
	if(kB==0) return fast_pow(p,a-((a-1)/A+1));
	int t=0;while(kB%p==0) kB/=p,++t;
	if(t%A) return 0;
	a-=t;int pa=fast_pow(p,a);
	int phi=pa-pa/p,g=getg(pa,phi),r=bsgs(g,kB,pa),d=gcd(phi,A);
	if(r%d) return 0;
	return d*fast_pow(p,t-t/A);
}

int main()
{
	scanf("%d",&T);
	while(T--) {
		scanf("%d%d%d",&A,&B,&P);
		P=P+P+1,kP=P,ans=1;
		for(RI i=2;i*i<=kP;++i) if(kP%i==0) {
			int js=0,pa=1;
			while(kP%i==0) kP/=i,++js,pa*=i;
			ans*=getans(i,js,B%pa);
		}
		if(kP!=1) ans*=getans(kP,1,B%kP);
		printf("%d\n",ans);
	}
    return 0;
}