[BZOJ2219]数论之神(BSGS)

题面

http://darkbzoj.tk/problem/2219

题解

前置知识

• BSGS算法： https://www.cnblogs.com/xh092113/p/12255049.html
• 求原根： https://blog.csdn.net/zhelong3205/article/details/78366476

我们要求方程(x^a = r({mod M}))（其中M是质数）在模M意义下的解x的数量。

首先将M的素因数分解式写出：(M = {prod_{i=1}^{k}}{p_i^{s_i}})，其中(p_1<p_2<…<p_k)为奇素数，(s_i{in}Z^*)。

那么我们可以将原方程分解为多个方程：

[egin{cases} x^a {equiv} r {mod {p_1}^{s_1}} \… \x^a {equiv} r {mod {p_k}^{s_k}} end{cases} ]

由中国剩余定理和乘法原理，这第一个方程在模({p_1}^{s_1})意义下的解数(A_1)、……第k个方程在模({p_k}^{s_k})意义下的解数(A_k)的乘积就是原方程的解数。

因此，现在只需要解决一个形如(x^a{equiv}r{mod p^s})的方程在模(p^s)意义下的解数的问题。

不妨设(r {in} [0,p^s))。可以分成三种情况讨论：

• 1、(r=0)：这代表着只要(p^{{lceil}{frac{s}{a}}{ ceil}} | x)就满足方程了，解数为(p^{s-{lceil}{frac{s}{a}}{ ceil}})。

• 2、(r { eq} 0)，设(r=p^t*u)，其中((u,p)=1)。

  • 2.1、(t=0)：那么x没有约数p。

    [x^a{equiv}u{mod p^s} ]

    设g为(p^s)的原根，再设(m,n{in}[0,phi(p^s)))使得(x=g^m,u=g^n)。（由于u,x均与p互质，所以满足条件的m,n一定存在）n可以通过BSGS算法求出，而我们想要求的答案转化为

    [am{equiv}n{mod {phi(p^s)}} ]

    在模({phi(p^s)})意义下的解数。那么取(d=(a,{phi(p^s)}))，若(d{ mid}n)则无解，否则由裴蜀定理，恰有d组解。

  • 2.2、(t{ eq}0)：那么此时x中含p的次数已经被确定。如果(a{ mid}t)则无解；否则，确定了有(v_p(x)={frac{t}{a}})。那么设(x=p^{frac{t}{a}}*v ((v,p)=1))，有

    [v^a {equiv} u {mod p^{s-t}} ]

    此时，方程转为了2.1的形式，可以用2.1的方法，计算出此处v在模(p^{s-t})意义下的解数N。

    由于(x=p^{frac{t}{a}}*v)，那么x在模(p^s)意义下的解数即为v在模(p^{s-{frac{t}{a}}})意义下的解数，也就是(N*p^{t-{frac{t}{a}}})。

代码

#include<bits/stdc++.h>

using namespace std;

#define ll long long
#define rg register 

inline ll read(){
	ll s = 0,ww = 1;
	char ch = getchar();
	while(ch < '0' || ch > '9'){if(ch == '-')ww = -1;ch = getchar();}
	while('0' <= ch && ch <= '9'){s = 10 * s + ch - '0';ch = getchar();}
	return s * ww;
} 

inline void write(ll x){
	if(x < 0)putchar('-'),x = -x;
	if(x > 9)write(x / 10);
	putchar('0' + x % 10);
}

ll pn;
bool isp[32000+5];
ll pri[32000+5];

inline void Eular(){
	pn = 0;
	for(rg ll i = 2;i <= 32000;i++)isp[i] = 1;
	for(rg ll i = 2;i <= 32000;i++){
		if(isp[i])pri[++pn] = i;
		for(rg ll j = 1;i * pri[j] <= 32000;j++){
			isp[i * pri[j]] = 0;
			if(i % pri[j] == 0)break;
		}
	}
}

inline ll power(ll a,ll n,ll mod){
	ll s = 1,x = a;
	if(!mod){
		while(n){
			if(n & 1)s *= x;
			x *= x;
			n >>= 1;
		}
		return s;
	}
	while(n){
		if(n & 1)s = s * x % mod;
		x = x * x % mod;
		n >>= 1;
	}
	return s;
}

ll k;
ll p[30+5],s[30+5],P[30+5];

inline void divide(ll M,ll &k,ll *p,ll *s,ll *P){ //素因数分解 
	k = 0;
	ll T = (ll)round(sqrt(M));
	for(rg ll i = 1;i <= pn && pri[i] <= T;i++){
		if(M % pri[i] == 0){
			k++;
			p[k] = pri[i];
			for(P[k] = 1,s[k] = 0;M % p[k] == 0;M /= p[k])P[k] *= p[k],s[k]++;
		}
	}
	if(M != 1){
		k++;
		p[k] = P[k] = M;
		s[k] = 1;
	}
}

inline ll gcd(ll a,ll b){
	return b ? gcd(b,a % b) : a;
}

ll _k;
ll _p[30+5],_s[30+5],_P[30+5];

inline ll calcg(ll p,ll P){ //计算P的原根 
	ll phi = P / p * (p - 1);
	divide(phi,_k,_p,_s,_P);
	for(rg ll g = 2;;g++){
		bool b = 1;
		for(rg ll i = 1;i <= _k;i++)if(power(g,phi/_p[i],P) == 1){
			b = 0;
			break;
		}
		if(b)return g;
	}
}

unordered_map<ll,ll>Hash;

inline ll BSGS(ll g,ll u,ll phi,ll P){ //计算方程g^x=u(mod P)在模phi(P)意义下的解x 
	Hash.clear();
	ll T = (ll)ceil(sqrt(P));
	ll g_T = power(g,T,P);
	for(rg ll cur = u * g % P,i = 1;i <= T;i++,cur = cur * g % P)Hash[cur] = i;
	for(rg ll cur = g_T,i = 1;i <= T;i++,cur = cur * g_T % P)
	if(Hash.count(cur))return i * T - Hash[cur];	
}

inline ll solve(ll a,ll u,ll p,ll P){ //计算方程x^a=u(mod P)在模P意义下解x的个数 
	ll g = calcg(p,P);
	ll phi = P / p * (p - 1);
	ll n = BSGS(g,u,phi,P); 
	ll d = gcd(a,phi);
	if(n % d)return 0;
	else return d;
}

inline ll calc(ll a,ll r,ll p,ll s,ll P){ //计算各A[i]。P=p^s 
	r %= P;
	if(!r)return power(p,s - (ll)ceil(1.0*s/a),0);
	ll t = 0,u;
	while(r % p == 0)r /= p,P /= p,t++; 
	u = r;
	if(!t)return solve(a,u,p,P);
	else{
		if(t % a)return 0;
		else return power(p,t - t / a,0) * solve(a,u,p,P); //此处的P已经变为p^(s-t)	
	}
}

int main(){
	Eular();
	ll T = read();
	while(T--){
		ll a = read(),r = read(),M = read() * 2 + 1;
		divide(M,k,p,s,P);
		ll ans = 1;
		for(rg ll i = 1;i <= k;i++)ans *= calc(a,r,p[i],s[i],P[i]);	
		write(ans),putchar('
');
	}
	return 0;
}