[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} ]

[{aover 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('
');
		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 ot=0)则答案取(x=0)，需要进行分类讨论；如果先进行(a\%=mod)，那么这两种情况就会被混淆。