[数学系列#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://acxblog.site/archives/exgcd-exbsgs.html
文章采用知识共享署名-非商业性使用 4.0 国际许可协议进行许可。