[数学系列#1] 扩展GCD/扩展BSGS
EXGCD
来放一下个人认为最精简的写法 ww:
void exgcd(ll a, ll b, ll &x, ll &y) {
if(b==0) x=1, y=0;
else exgcd(b, a%b, y, x), y-=a/b*x;
}
BSGS/EXBSGS
考虑原方程: $y^x \equiv z \pmod {p}$
由欧拉定理,在 $y$ 与 $p$ 互质时,$y^{x} \equiv y^{x\% \phi(p)} \pmod {p}$。
因此只要 $y$ 与 $p$ 互质,我们就可以大胆放心的只需讨论 $x$ 从 $1$ 到 $p-1$ 的解情况。若您比较闲,也可以线筛或公式求一发 $\phi(p)$。
然后令 $x=im-j$ ,则有: $y^{im}\equiv y^jz \pmod {p}$
令 $m=\left \lceil \sqrt{p} \right \rceil$,则预处理 $y^jz \% p$ 塞到哈希表就行。
但是 $y$ 与 $p$ 不互质的话,就只能这么推:
$$g=gcd(y, p)$$
$$\Rightarrow y^{x-1}\equiv (\frac{y}{g})^{-1} \frac {z}{g} \pmod {\frac{p}{g}}$$
由于本质上是一个同余方程,若 $g\nmid z$ ,则方程无解
然后递归处理就可以辣
两个复杂度基本上都是 $O(\sqrt{m})$ ~
下面是代码:
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
void exgcd(ll a, ll b, ll &x, ll &y) {
if(b==0) x=1, y=0;
else exgcd(b, a%b, y, x), y-=a/b*x;
}
ll inv(ll x, ll mo) {
ll a, b;
exgcd(x, mo, a, b);
return (a%mo+mo)%mo;
}
const int N=1e5, M=1e6+10;
int h[N], nexp[M], p=1;
ll w[M], v[M];
inline void ins(int a, int b, int c) { nexp[p]=h[a], h[a]=p, w[p]=b, v[p]=c, p++; }
void clr() { memset(h, 0, sizeof(h)), p=1; }
int exbsgs(ll y, ll z, ll p) {
if(z==1 || p==1) return 0;
ll g=__gcd(y, p);
if(g>1) {
if(z%g) return -0x3f3f3f3f;
z/=g, p/=g, z=z*inv(y/g, p)%p;
return exbsgs(y, z, p)+1;
}
int unit=ceil(sqrt(p));
clr();
ll yz=1;
for(int i=0;i<unit;i++, yz=yz*y%p) ins(yz*z%p%N, yz*z%p, i);
ll ny=yz;
for(int i=1;i<=unit;i++, ny=ny*yz%p)
for(int u=h[ny%N];u;u=nexp[u])
if(ny==w[u]) return i*unit-v[u];
return -0x3f3f3f3f;
}
int main() {
ios::sync_with_stdio(false), cin.tie(0);
ll a, p, b;
while(cin>>a>>p>>b, a || p || b) {
ll ans=exbsgs(a, b, p);
if(ans<0) cout<<"No Solution"<<endl;
else cout<<ans<<endl;
}
return 0;
}
本文链接:https://pst.iorinn.moe/archives/exgcd-exbsgs.html
许可: https://pst.iorinn.moe/license.html若无特别说明,博客内的文章默认将采用 CC BY 4.0 许可协议 进行许可☆