[SPOJ3105][洛谷P4195]Mod(exBSGS模板)

题面

https://www.luogu.com.cn/problem/P4195

题解

需要用到exBSGS（扩展BSGS）算法。

BSGS算法见 https://www.cnblogs.com/xh092113/p/12255049.html

exBSGS可以快速解决一般的\(a^x{\equiv}r{\mod p}\)的指数同余方程。

如果方程中的\(gcd(a,p)=1\)，那么可以用BSGS解决；

否则设\(g=gcd(a,p)\)，可以如下处理：

\[a^x{\equiv}r{\mod p} \]

\[{a\over g}*a^{x-1}{\equiv}{r \over g}{\mod {\frac}{p}{g}} \]

\[a^{x-1}{\equiv}{\frac{r}{g}}*{(\frac{a}{g})^{-1}}{\mod {\frac{p}{g}}} \]

• 这里\(({\frac{a}{g}})^{-1}\)指的是\({\frac{a}{g}}\)关于\({\frac{p}{g}}\)的逆元

接下来就可以递归地做了。因为\(p\)变成了\({\frac{p}{g}}\)，\(p\)在不断变小，所以最后总会使得\(gcd(a,p)=1\)，从而使用BSGS解决。

时间复杂度方面，递归的层数最多是\(O(\log{p})\)，每层求gcd和逆元是\(O(\log{p})\)，最后一层BSGS是\(O(\sqrt{p})\)，故总时间复杂度是\(O(\log^2{p}+\sqrt{p})\)。

代码

#include<bits/stdc++.h>

using namespace std;

#define ll long long
#define rg register
#define inf 0x3f3f3f3f3f3f3f3f

namespace ModCalc{
	inline void Inc(ll &x,ll y,ll mod){
		x += y;if(x >= mod)x -= mod;
	}
	
	inline void Dec(ll &x,ll y,ll mod){
		x -= y;if(x < 0)x += mod;
	}
	
	inline void Adjust(ll &x,ll mod){
		x = (x % mod + mod) % mod;
	} 
	
	inline ll Add(ll x,ll y,ll mod){
		Inc(x,y,mod);return x;
	}
	
	inline ll Sub(ll x,ll y,ll mod){
		Dec(x,y,mod);return x;
	}
	
	inline ll check(ll x,ll mod){
		Adjust(x,mod);return x;
	}
}
using namespace ModCalc;

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);
}

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

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

inline void exgcd(ll a,ll &x,ll b,ll &y){
	if(!b){
		x = 1,y = 0;
		return;
	}
	exgcd(b,y,a%b,x);
	y -= a / b * x;
}

unordered_map<ll,ll>H;

inline ll BSGS(ll a,ll r,ll mod){
	a %= mod,r %= mod;
	ll T = (ll)round(sqrt(mod));
	ll a_T = power(a,T,mod);
	H.clear();
	for(rg ll i = 1,cur = r * a % mod;i <= T;i++,cur = cur * a % mod)
		H[cur] = i;
	for(rg ll i = 1,cur = a_T;i * T - T < mod;i++,cur = cur * a_T % mod)
	if(H.count(cur))return i * T - H[cur];
	return -inf;
}

inline ll exBSGS(ll a,ll r,ll mod){
	a %= mod,r %= mod;
	ll g = gcd(mod,a);
	if(r % g != 0){
		if(r == 1)return 0;
		else return -inf;
	}
	if(g == 1)return BSGS(a,r,mod);
	else{
		ll iv,y;
		exgcd(a/g,iv,mod/g,y);
		return exBSGS(a,r/g*check(iv,mod/g),mod/g) + 1;
	}
}

int main(){
	ll a = read(),mod = read(),r = read();
	while(a || r || mod){
		ll x = exBSGS(a,r,mod);
		if(x < 0)puts("No Solution");
		else write(x),putchar('\n');
		a = read(),mod = read(),r = read();
	}
	return 0;
}

• 第\(88\)行的\(mod\)操作要注意，如果此题中\(a\)可以为\(0\)，那么需要加一手\(a=0\)的特判。原因是：如果方程形如\(a^x{\equiv}1{\mod p}\)，而又有\(p|a\)成立时，此时若\(a=0\)则无解，若\(a\not=0\)则答案取\(x=0\)，需要进行分类讨论；如果先进行\(a\%=mod\)，那么这两种情况就会被混淆。