BSGS (北上广深) 算法，用来解决一类 a ^ x\equiv b \pmod p (即离散对数) 的问题。

---

1. 基础 BSGS

我们考虑 a 与 p 互质的情况。既然 a 与 p 互质，我们不难想到欧拉定理 a^{\varphi{}(p)} \equiv 1\pmod p，所以 x 的 取值范围为 [1,\varphi{}(p)]，直接枚举时间复杂度 O(p)

BSGS 算法的灵魂在于，对整个枚举的过程进行了优化。

我们考虑把原式变成这样的形式 a^{Ax - B} \equiv b\pmod p，由此可得 B 的取值范围为 [0,x)。这里我们显然可以把 a^B 除过去，得到我们想要的目标式子：

a^{Ax} \equiv b \times a ^ B\pmod p

注意到了 a,b 均是定值，所以我们可以考虑固定三个变量中的两个，然后枚举。所以我们接下来考虑 A,B 之间取值范围的关系。

我们有 Ax-B\in [1,\varphi{}(p)]，目标是让 A+B 最小化。根据基本不等式知识易得当 x = \lceil \sqrt{\varphi{}(p)} \rceil 时 A+B 最小。当然我们这里为了避免计算 \varphi，所以我们可以取 A = \lceil \sqrt{p} \rceil,显然这个 A 取的范围是覆盖原本的计算范围的。

接下来就是整个算法最精华的部分了，我们显然要对左右两边进行枚举，然后找有没有相等的值。我们先对左边的部分进行枚举，并存进一个 hash 表内，每个值对应的是 A 的大小(非 OI 赛场上完全可用 unordered_map 代替 hash 表)

unordered_map <ll,ll> hash;
ll tmp = b * a % p,t = sqrt(p) + 1;
for (ll i = 1;i <= t;++i) {
    hash[tmp] = i; 
    //printf("%lld %lld\n",tmp,i);
    tmp = (tmp * a) % p;
}

然后，我们类似的对右边进行枚举，并且在 hash 表内查询，最后算出答案。

ll at = ffpow(a,t,p),now = at;
for (ll i = 1;i <= t;++i) {
    if (hash[now]) {
        return i * t - hash[now];
    }
    now = (now * at) % p;
}

时间复杂度 O(\sqrt{p} \log p)

Code:

ll ffpow(ll a,ll b,ll p) {
    ll ans = 1;
    for (;b;b >>= 1) {
        if (b & 1) ans = ans * a % p;
        a = a * a % p;
    }
    return ans % p;
}

ll BSGS(ll a,ll b,ll p) {
    unordered_map <ll,ll> hash;
    ll tmp = b * a % p,t = sqrt(p) + 1;
    for (ll i = 1;i <= t;++i) {
        hash[tmp] = i; 
        //printf("%lld %lld\n",tmp,i);
        tmp = (tmp * a) % p;
    }
    ll at = ffpow(a,t,p),now = at;
    for (ll i = 1;i <= t;++i) {
        if (hash[now]) {
            return i * t - hash[now];
        }
        now = (now * at) % p;
    }
    return -1;
}

2. 拓展 BSGS

接下来我们考虑 a 和 p 不互质的情况。首先我们需要先确定是否有解的问题，我们可以把原本的同余方程转化成

a ^ x + p\times y = b \\ a \times a ^{x-1} + p \times y = b

根据贝祖定理，方程有解当且仅当 b | gcd(a,p)。

很自然想到把不互质的情况转化成互质的情况。我们设 d = gcd(a,p)，将上面方程中的系数不断除以 d，假设我们这个过程进行了 cnt 次，我们的式子化成这样

\frac{a}{d^{cnt}} \times a ^{x-cnt-1} + \frac{p}{d^{cnt}} \times y = \frac{b}{d^{cnt}}

再将其转回我们熟悉的同余方程形式就是

a ^{x-cnt-1} \equiv \frac{b}{d^{cnt}}\pmod{ \frac{p}{d^{cnt}}}

这样我们的 a 显然与 \frac{p}{d^{cnt}} 互质。

所以对之后的式子再变化一下：

a ^{Ax} \equiv \frac{b}{d^{cnt}} \times a ^ {b} \times a ^{cnt+1} \pmod {\frac{p}{d^{cnt}}}

套用一次 BSGS 算法即可，本质上是一样的。

Code:

ll exBSGS(ll a,ll b,ll p) {
    a %= p,b %= p;
    if (p == 1) return 0;//注意你当 p = 1 时是有解的
    if (b == 1) return 0;
    ll g,cnt = 0,c = 1;
    while ((g = gcd(a,p)) > 1) {
        if (b % g) return -1;
        b /= g,p /= g,c = c * (a / g) % p,++cnt;
        //这里的c实际上是\frac{a}{g} ^ cnt，之后算右半边的时候回用到
        if (b == c) return cnt;
    }
    unordered_map <ll,ll> hash;
    ll tmp = b * a % p,t = sqrt(p) + 1;
    for (ll i = 1;i <= t;++i) {
        hash[tmp] = i; 
        //printf("%lld %lld\n",tmp,i);
        tmp = (tmp * a) % p;
    }
    ll at = ffpow(a,t,p),now = c * at % p;
    for (ll i = 1;i <= t ;++i) {
        if (hash.count(now)) {
            return i * t - hash[now] + cnt;//注意加上次数
        }
        now = (now * at) % p;
    }
    return -1;
}

3.例题

相信大家一定会马上去把模板题秒了，所以这里就不放了。

P3306 [SDOI2013] 随机数生成器

给定序列 x_{i+1} \equiv a \times x_i + b\pmod p，求最小的 i 使得 x_i = t 成立

0 \leq a, b, x_1, t \lt p，2 \leq p \leq 10^9

p 为质数。

首先我们可以先把前几项推一下：

x_2 = a \times x_1 + b \\ x_3 = a^2 \times x_1 + a \times b + b \\ x_4 = a^3 \times x_1 + a ^ 2 \times b + a \times b + b

所以我们运用一点等比数列的知识，可以推出通项公式。

x_n = a ^{n-1} \times x_1 + b \times \frac{1-a^{n-1}}{1-a}

接下来就是激情推式子！

\begin{aligned} x_n = a ^{n-1} \times x_1 + b \times \frac{a^{n-1}-1}{a-1} &\equiv t \pmod p \\ a ^ {n-1} \times x_1 + \frac{b}{a-1} \times a ^ {n-1} - \frac{b}{a-1} &\equiv t \pmod p \\ \frac{ax_1+b-x_1}{a-1} \times a ^ {n-1}&\equiv t + \frac{b}{a-1} = \frac{at - t + b}{a-1}\pmod p \\ a^{n-1} &\equiv \frac{at-t+b}{ax_1+b-x_1}\pmod p\\ a^{n-1} &\equiv \frac{t-at-b}{x_1-ax_1-b}\pmod p \end{aligned}

注意到 p 是质数，所以可以直接用普通 BSGS 求即可。

代码实现上细节比较多，自行查看。

Code:

void solve() {
    ll p,a,b,x1,t;read(p,a,b,x1,t);
    if (x1 == t) {puts("1");return ;}
    if (a == 0) {
        if (b == t) {puts("2");}
        else {puts("-1");}
        return;
    }
    if (a == 1) {
        if ((t - x1) % gcd(b,p)) {
            puts("-1");
            return ;
        }
        t = (t - x1 + p) % p;
        printf("%lld\n",(t * ffpow(b,p-2,p)) % p + 1);
        return ;
    }
    ll r = ((t - a * t - b) % p * ffpow((x1 - x1 * a - b) % p,p-2,p) % p) % p;
    ll ans = BSGS(a,r,p);
    if (ans == -1) puts("-1");
    else printf("%lld\n",ans + 1);
}