BSGS算法lyh's blog

BSGS算法（Baby-Step-Giant-Step，大步小步法，简称BSGS），~~也简称北上广深算法~~。是一类用于求解高次同余方程的算法。

引入：现在需要解这么一个同余方程：

$A^x\equiv B(mod\ C)$

（这里C是质数）要求使得上式成立的最小x。一种方法是暴力测试，复杂度线性，但是当实际答案特别大时效率低下。BSGS算法可以将其压到根号级别。
这里需要A与C互质。假设$m=\lceil \sqrt C\rceil$（注意这里是向上取整），然后以m为除数，根据带余除法原理可以将x写成$x=im+j$的形式，这样同余式就是：

$A^{im+j}\equiv B(mod\ C)$

既然A与C互质，故A的幂次与C也应该互质，那么A的幂次关于C的逆元存在，用它的i次方去乘同余式两侧，得到：

$A^j\equiv B(A^{-m})^i(mod\ C)$

        考虑到j的范围是[0, m)中的整数，可以先枚举这个区间中的值，求得数对$(A^j\%C,j)$，把这些数对用hash表存起来。注意对于多个j满足$A^j\%C$为同一个值时，记录最小的j，这是为了保证最后求出的x最小。这一步为算法中的“小步”，时间复杂度$O(\sqrt C)$。
        然后枚举右边的i，对于每一个i，求出$B(A^{-m})^i$，然后到hash表中找对应的j，答案即为$x=im+j$，这一步为算法中的“大步”。
        i的范围如何确定？这里i的取值范围定为$[0,m)$中的整数。取这个范围的原因是如果i在这个范围内都无解，那么$i≥m$时必然也无解。根据j的范围和i的范围，可知x可能的区间是$[0,m(m-1)+m-1]$。这个区间的右端点$m^2-1$一定不比C-1小。根据欧拉定理有$A^{\phi(C)}\equiv 1(mod C)$，可知x不可能大于$\phi(C)$，而当C为质数时，$\phi(C)$不会超过C-1，因此所有可能的范围我们都考虑到了。
        伪模板题。
        就是求：

$\frac {10^x-1} {9}\equiv k(mod\ m)$

m不是3的时候，这个式子等价于：

$10^x-1\equiv 9k(mod\ m)$

那么就是：

$10^x\equiv 9k+1(mod\ m)$

当m不是2和5的时候（此时m与10不互质），用BSGS算法解决即可。对于m=2、3、5的情况特判。在乘过程中可能会爆long long，需要用快速乘。

#include<iostream>
#include<cmath>
#include <map>

using namespace std;
typedef long long ll;
ll k, m;
map<ll, ll> hs;

ll qMutiple(ll a, ll b) {//快速乘
    ll ans = 0, s = a;
    while (b) {
        if (b & 1)ans = (ans + s) % m;
        b >>= 1, s = (s << 1) % m;
    }
    return ans;
}

ll qPow(ll a, ll b) {//快速幂
    ll ans = 1ll, s = a;
    while (b) {
        if (b & 1)ans = qMutiple(ans, s);
        b >>= 1, s = qMutiple(s, s);
    }
    return ans;
}

ll BSGS(ll x, ll p) {
    ll s = static_cast<ll>(sqrt(m * 1.0) + 2), ps = 1ll, tp = 1ll;
    for (ll i = 0; i < s; i++, ps = qMutiple(ps, x))if (!hs.count(ps))hs[ps] = i;
    ps = qPow(x, m - s - 1);//Fermat小定理:逆元
    for (ll i = 0; i < s; i++, tp = qMutiple(tp, ps)) {
        if (hs.count(qMutiple(tp, p)))return i * s + hs[qMutiple(tp, p)];
    }
    return -1;//无解-1
}

int main() {
    cin >> k >> m;
    if (m == 2) {//特判
        k %= 2;
        if (k == 0)cout << -1;
        else cout << 1;
        return 0;
    } else if (m == 3) {
        k %= 3;
        if (k == 0)cout << 3;
        else if (k == 1)cout << 1;
        else cout << 2;
        return 0;
    } else if (m == 5) {
        k %= 5;
        if (k != 1)cout << -1;
        else cout << 1;
        return 0;
    }
    cout << BSGS(10, (9 * k + 1) % m);
    return 0;
}

exBSGS算法

当A与C不互质（或者说C为非质数时）BSGS算法是不适用的，此时就需要使用exBSGS算法。观察下面的同余式：

$A^x\equiv B(mod\ C)$

这里有一个结论：若$d=gcd(A,C)$，则当$d\not|C$且$B≠1$时，方程无自然数解，$B=1$时有解0。
有解时，将同余式两侧同时去除d，得到：

$A^{x-1}\frac {A} {d}\equiv \frac {B} {d}(mod\ \frac {C} {d})$

然后就有：

$A^{x-1}\equiv \frac {B} {d} (\frac {A} {d})^{-1}(mod\ \frac {C} {d})$

反复使用这种变换手段，直到满足模数与$A$互质为止，此时用BSGS算法解决。模板题：戳这里，注意这题不要用快速乘（~~一点也不快~~）。

#include<bits/stdc++.h>

using namespace std;
typedef long long ll;

ll gcd(ll x, ll y) {
    if (y == 0)return x;
    return gcd(y, x % y);
}

ll exGcd(ll a, ll b, ll &x, ll &y) {//扩展欧几里得
    if (b == 0) {
        x = 1, y = 0;
        return a;
    }
    ll t = exGcd(b, a % b, y, x);
    y -= a / b * x;
    return t;
}

unordered_map<ll, ll> hs;//比map快

ll BSGS(ll x, ll p, ll m) {
    hs.clear();
    ll s = static_cast<ll>(sqrt(m * 1.0) + 2), ps = 1ll, tp = 1ll, a, b;
    for (ll i = 0; i < s; i++, ps = ps * x % m)if (!hs.count(ps))hs[ps] = i;
    exGcd(ps, m, a, b), ps = (a % m + m) % m;//a为逆元，注意正数化
    for (ll i = 0; i < s; i++, tp = tp * ps % m) {
        if (hs.count(tp * p % m))return i * s + hs[tp * p % m];
    }
    return -1;//无解-1
}

ll exBSGS(ll a, ll b, ll p) {
    ll d, k = 0, x, y, ans;
    while ((d = gcd(a, p)) != 1) {
        if (b % d)return b == 1 ? 0 : -1;
        b /= d, p /= d, ++k, exGcd(a / d, p, x, y), x = (x % p + p) % p, b = b * x % p;
    }
    return (ans = BSGS(a, b, p)) == -1 ? -1 : ans + k;
}

int main() {
    int a, p, b;
    while (cin >> a >> p >> b) {
        if (a == 0 && b == 0 && p == 0)return 0;
        ll x = exBSGS(a, b, p);
        if (x == -1)cout << "No Solution" << endl;
        else cout << x << endl;
    }
}