exBSGS
已知数(a,p,b)，求满足(a^x≡b (mod p))的最小自然数(x)。
(100\%)的数据，(a,p,b≤10^9)。
_皎月半洒花的题解
其实本质上，当(p)不为素数时，我们无法进行朴素 BSGS 的原因是我们的欧拉定理(a^{varphi(p)} equiv b(mod p)) 只能处理((a,p)=1)的情况。那么我们知道，朴素的 BSGS 的关键在于，可以保证最小解是有界的——(x)一定在([1,varphi(p)])中。所以最后 BSGS 的复杂度才会是(Theta(sqrt{varphi(p)})) 的——比如说比较常见的(p)是素数的情况下，时间复杂度为(Theta(p))。
那么也就是说，我们只需要进行一些操作，保证((a,p)=1)即可。
我们思考，对于同余式(a^xequiv b~(mod p))而言，我们先假定((a,p)>1)。而此时如果有(((a,p), b)=1)，那么说明此式只有可能在(x=0,b=1)的时候有解——这个结论是平凡的。因为假设我们把它展开成(acdot a^{x-1} +kp=b)的形式，必须要有((a,p) ~|~ b)的情况下，才能保证(a^{x-1})和(k)都是整数。
那么对于((a,p)>1)且((a,p)~|~b)，我们令原式变成
[ a^{x-1}cdot frac{a}{(a,p)} equiv frac{b}{(a,p)} (mod frac{p}{(a,p)}) $$ 的样子，如果此时$(a^{x-1},frac{p}{(a,p)})=1$ 的话，我们就直接解
$$ a^{x-1}equiv frac{frac{b}{(a,p)}} {frac{a}{(a,p)} }(mod frac{p}{(a,p)}) $$ 这个方程即可。否则我们继续分解直至$(p',a)=1$。
那么此时有个问题需要注意，就是如果们在解这个方程时，出现了
$$ (a^{x-1}, frac{p}{(a,p)})
mid frac{frac{b}{(a,p)}} {frac{a}{(a,p)} } $$ 的情况，那我们需要特判并`return -1` ；另一种情况，如果我们出现了
$$ a^{x-1}equiv frac{frac{b}{(a,p)}} {frac{a}{(a,p)} } equiv1(mod frac{p}{(a,p)}) $$ 的情况，也需要特判并输出此$k$（此时同余式左边是$a^{x-k}$，因为$a^{x-k}equiv1~(mod p)$所以直接输出$k$），不过也有可能不需要，完全看你写的 BSGS 能不能判断$x=0$的情况……一般情况下不能。
算法流程大体就是这些。有些需要注意的是，BSGS 有两种写法。我的这种写法是小步为$a^j$，大步为$a^{-im}b$。此时由于$p$不再是素数，所以不能用费马小定理，需要我们用exgcd的方法求逆元，包括但不限于$frac{b}{(a,p)}$的逆元和$a^{-im}$。
```cpp
#include<bits/stdc++.h>
#define co const
#define il inline
template<class T> T read(){
T x=0,w=1;char c=getchar();
for(;!isdigit(c);c=getchar())if(c=='-') w=-w;
for(;isdigit(c);c=getchar()) x=x*10+c-'0';
return x*w;
}
template<class T>il T read(T&x){
return x=read<T>();
}
using namespace std;
typedef long long LL;
unordered_map<LL,int> H;
int N,M,P,ans;
LL gcd(LL a,LL b){
if(!b) return a;
return gcd(b,a%b);
}
LL expow(LL a,LL b,LL mod){
LL res=1;
for(;b;b>>=1,a=a*a%mod)
if(b&1) res=res*a%mod;
return res;
}
LL exgcd(LL&x,LL&y,LL a,LL b){
if(!b){
x=1,y=0;
return a;
}
LL t=exgcd(y,x,b,a%b);
y-=x*(a/b);
return t;
}
LL BSGS(LL a,LL b,LL mod,LL d){
H.clear();
LL Q,p=ceil(sqrt(mod)),x,y;
exgcd(x,y,d,mod),b=(b*x%mod+mod)%mod;
Q=expow(a,p,mod),exgcd(x,y,Q,mod),Q=(x%mod+mod)%mod;
for(LL i=1,j=0;j<=p;++j,i=i*a%mod)
if(!H.count(i)) H[i]=j;
for(LL i=b,j=0;j<=p;++j,i=i*Q%mod)
if(H.count(i)) return j*p+H[i];
return -1;
}
LL exBSGS(){
if(M==1) return 0;
LL d=1,k=0;
for(LL nd;(nd=gcd(N,P))>1;){
if(M%nd) return -1;
++k,M/=nd,P/=nd,d=d*(N/nd)%P;
if(d==M) return k;
}
LL ans;
return (ans=BSGS(N,M,P,d))==-1?-1:ans+k;
}
int main(){
while(read(N)|read(P)|read(M)){
N%=P,M%=P,ans=exBSGS();
if(ans<0) puts("No Solution");
else printf("%d
",ans);
}
return 0;
}
```]