bzoj 2219: 数论之神 数论

       首先将2*k+1分解质因数，然后用天朝剩余定理合并答案。具体就是每个方程解的个数的乘积。

       然后看方程x^a≡b(mod p^q)的解的个数，当b=0时：

       x^a≡0(mod p^q)，因此显然x必须是p^[q/a]的倍数，这里[x]表示x向上取整。那么就相当于统计0~p^q-1中p^[q/a]的倍数的个数；

       x^a≡b(mod p^q)，其中b=c*p^r，那么显然有a|r，否则无解；那么x应该为a'*p^(r/a)的形式，两边除以p^r后取指标变为a' lndx≡lndc(mod p^(q-r))，用BSGS得到indc进而得到答案。注意此时的答案并不是该方程的答案，因为x,x+p^(q-r),x+2*p^(q-r)都是原方程的答案，因此需要乘上p^(r-r/a)。

AC代码如下：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#define ll long long
#define N 100005
using namespace std;

int c[N];
int ksm(int x,int y){
	int t=1; for (; y; y>>=1,x*=x) if (y&1) t*=x; return t;
}
int ksm_p(int x,int y,int p){
	int t=1; for (; y; y>>=1,x=(ll)x*x%p) if (y&1) t=(ll)t*x%p; return t;
}
int gcd(int x,int y){ return (y)?gcd(y,x%y):x; }
int p_rt(int p){
	int i,j,cnt=0,q=p-1;
	for (i=2; i*i<=q; i++) if (!(q%i)){
		c[++cnt]=i; if (i*i!=q) c[++cnt]=q/i;
	}
	for (i=2; i<p; i++){
		for (j=1; j<=cnt; j++)
			if (ksm_p(i,c[j],p)==1) break;
		if (j>cnt) return i;
	}
}
struct hash_cnt{
	int tot,fst[100004],pnt[N],edg[N],nxt[N];
	void clr(){
		tot=0; memset(fst,0,sizeof(fst));
	}
	void ins(int x,int y){
		int t=x%100003+1,p;
		for (p=fst[t]; p; p=nxt[p]) if (pnt[p]==x) return;
		pnt[++tot]=x; edg[tot]=y; nxt[tot]=fst[t]; fst[t]=tot;
	}
	int find(int x){
		int t=x%100003+1,p;
		for (p=fst[t]; p; p=nxt[p]) if (pnt[p]==x) return edg[p];
		return 0;
	}
}hsh;
int calc(int x,int y,int p,int q){
	int mod=ksm(p,q); y%=mod;
	if (!y) return ksm(p,q-(q-1)/x-1);
	int cnt=0; for (; !(y%p); y/=p) cnt++;
	if (cnt%x) return 0;
	q-=cnt; cnt/=x; int phi=mod-mod/p,g=p_rt(p);
	int m=(int)sqrt(mod)+1,now=1,i;
	hsh.clr();
	for (i=1; i<=m; i++){
		now=(ll)now*g%mod; hsh.ins(now,i);
	}
	int w=ksm_p(now,phi-1,mod),k; now=y;
	for (i=0; i<m; i++,now=(ll)now*w%mod)
		if (k=hsh.find(now)) break;
	int ind=i*m+k,d=gcd(x,phi);
	return (ind%d)?0:d*ksm(p,(x-1)*cnt);
}
int main(){
	int x,y,p,cas; scanf("%d",&cas);
	while (cas--){
		scanf("%d%d%d",&x,&y,&p); p=p<<1|1; int i,ans=1;
		for (i=2; i*i<=p; i++) if (!(p%i)){
			int cnt=0; for (; !(p%i); p/=i) cnt++;
			if (ans) ans*=calc(x,y,i,cnt);
		}
		if (ans && p!=1) ans*=calc(x,y,p,1);
		printf("%d\n",ans);
	}
	return 0;
}

by lych

2016.4.18