[学习笔记]BSGS及其扩展

BSGS

概述

BSGS算法是用于求解同余方程\(a^x \equiv b \ (mod \ p)\)(其中\(a,p\)互质)的最小自然数解(或正整数解)高效算法，复杂度\(O(\sqrt p)\)

思路

首先由欧拉定理，有\(a^{\varphi (p)} \equiv 1 \ (mod \ p)\)，于是\(a^{x} \equiv a^{x \mod \varphi (p)} \ (mod \ p)\)。

所以\(a^x \equiv b \ (mod \ p)\)如果有解，一定存在一个满足\(x < \varphi (p) < p\)的解，我们只需考虑\(x < p\)的部分。

对于\(a^x \equiv b \ (mod \ p)\), 令\(x = i \cdot m - j\), 其中\(i = ceil(p / m), j < m\)。

则原方程可化为\(a^{i \cdot m - j} \equiv b \ (mod \ p)\)，或\(a^{i \cdot m} \equiv ba^{j} \ (mod \ p)\)。

注意到\(i\)和\(j\)的取值均不超过\(m\)个，所以可以暴力枚举\(j\)，存下右边的值对应的\(j\)，再暴力枚举\(i\)，查询右边能否取到\(a^{im}\)。如果能，那么答案就是对应的\(i \cdot m - j\)，否则就继续枚举。

如果枚举完所以\(i,j\)还没找到答案，说明无解。

容易发现上面的过程中\(i \cdot m - j\) 不会取到0，所以如果求自然数解，需特判\(b = 1\)的情况。

上面的过程中枚举\(i\)的复杂度\(O(p / m)\)，枚举\(j\)的复杂度\(O(m)\)。为了使总复杂的最小，应取\(m = \sqrt{p}\)，总复杂度\(O(\sqrt{p})\)。

存等号右边的取值时，可以手写hash或者直接用map，用map复杂度会多一个\(log\)。

实践中可以map和unordered_map都试一下，万一出题人不卡呢

代码在最下方。

扩展BSGS

概述

普通的BSGS只能解决\(a,p\)互质的情况，那么\(a,p\)不互质的情况怎么办呢？

思路

\(a^x \equiv b \ (mod \ p)\)可以表示成\(a^x + k \cdot p = b, k为整数\)，假设\(gcd(a, p) = d\)，上面的式子有解的条件是\(d | b\)，并且可化为\(a^{x - 1} \frac{a}{d} + k \cdot \frac{p}{d} = \frac{b}{d}\)，即解方程\(a^{x - 1} \frac{a}{d} \equiv \frac{b}{d} \ (mod \ \frac{p}{d})\)。

于是我们可以不断的进行上面的转化，直到\(gcd(a, \frac{p}{\prod_{i}{d_i}}) = 1\)，这时式子变为\(a^{x - k} \frac{a^k}{\prod_{i = 1}^{k}{d_i}} \equiv \frac{b}{\prod_{i = 1}^{k}{d_i}} \ (mod \ \frac{p}{\prod_{i}{d_i}})\)，把\(\frac{a^k}{\prod_{i=1}^{k}{d_i}}\)除到右边去(乘逆元)，就变成了BSGS能解的问题，原问题的答案就是BSGS的结果加上\(k\)，BSGS无解则原方程无解。

如果中间某一步发现\(\frac{a^k}{\prod_{i = 1}^{k}{d_i}} \equiv \frac{b}{\prod_{i = 1}^{k}{d_i}} \ (mod \ \frac{p}{\prod_{i}{d_i}})\)了，答案就是\(k\)。

代码：

//qpow(a, b, p)=a的b次幂模p
//inverse(a, p)=a在模p意义下的逆元
struct BSGS {
    std::map<LL, LL> hash;
    LL solve(LL a, LL b, LL p) {
        a %= p, b %= p;
        if (b == 1) return 0;
        hash.clear();
        LL m = sqrt(p) + 1;
        for (LL i = 0, j = b; i < m; ++i, j = j * a % p)
            hash[j] = i;
        LL A = qpow(a, m, p);
        for (LL i = 1, j = A; i <= m; ++i, j = j * A % p)
            if (hash.find(j) != hash.end()) return i * m - hash[j];
        return -1;
    }
    LL solve_ex(LL a, LL b, LL p) {
        a %= p, b %= p;
        if (b == 1 || p == 1) return 0;
        LL A = 1, k = 0;
        for (LL d = gcd(a, p); d != 1; d = gcd(a, p)) {
            if (b % d) return -1;
            p /= d, b /= d, ++k;
            A = a / d * A % p;
            if (A == b) return k;
        }
        LL res = solve(a, b * inverse(A, p) % p, p);
        return (res == -1 ? -1 : res + k);
    }
} bsgs;