高次同余方程之——BSGS算法

BSGS算法

BSGS算法，即大小步算法（Baby-Step-Giant-Step）俗称“北上广深、拔山盖世“算法。

用途

可以在\(O(\sqrt{P})\)的时间复杂度求出形如 \(A^x\equiv B(mod\;P),\,(gcd(P,A)=1)\) 的方程的最小非负整数解。

实现

根据费马小定理，我们不难证明若最小非负整数 \(x\) 存在，那么一定有 \(x<P\)，我们不妨假设

\[x=kn-m,\,A^{kn-m}\equiv B(mod\,P) \]

其中 \(n\) 为 \(\lceil \sqrt{P}\,\rceil\) ，\(k、m\) 为不超过 \(n\) 的正整数。

这样做有什么好处呢？我们发现 \(k、m\) 都是取值范围小于 \(\lceil \sqrt{P}\,\rceil\) 的正整数，我们想要利用好这条性质，就要把式子变一下形:

\[A^{kn−m}\equiv B(mod\;P) \]

\[A^{kn}\equiv A^m⋅B(mod\;P) \]

我们令 \(G=A^n( mod\,P )\)，则

\[G^k\equiv A^m⋅B(mod\;P) \]

这样，我们就可以先处理每一个$ A^{m}\cdot B $，存进哈希表中，再求出 \(G=A^n\)，并依次枚举\(\,G^k\,\)在哈希表中是否存在，优先取 \(k\) 小的即可，最终就有 \(x=kn-m\)。

扩展BSGS

不难看出，以上过程依托费马小定理限定了 \(x\) 的取值范围，所以限制是\((gcd(P,A)=1)\)，但如果不是这样呢？

我们不妨假设 \((gcd(P,A)=d),A=a\times d，B=b\times d，P=p\times d\)（如果 \(B\) 不是 \(d\) 的倍数则方程在整数范围内无解）则：

\[(a⋅d)^x\equiv b⋅d(mod\;p⋅d) \]

根据等式的性质，一定有

\[a⋅(a⋅d)^{x−1}\equiv b(mod\;p) \]

\[(a⋅d)^{x−1}\equiv b⋅a^{−1}(mod\;p) \]

依照上文求解即可。

例子：SDOI2015计算器

题目描述

你被要求设计一个计算器完成以下三项任务：

1、给定 \(y,z,p\), 计算 \(y^z\ mod\ p\) 的值；

2、给定 \(y,z,p\)，计算满足 \(xy≡ z ( mod\ p)\) 的最小非负整数；

3、给定 \(y,z,p\)，计算满足 \(y^x ≡ z ( mod\ p)\) 的最小非负整数。

输入格式

输入包含多组数据。
第一行包含两个正整数 \(T,K\) 分别表示数据组数和询问类型（对于一个测试点内的所有数据，询问类型相同）。
以下行每行包含三个正整数 \(y,z,p\)，描述一个询问。

输出格式

对于每个询问，输出一行答案。对于询问类型 2 和 3，如果不存在满足条件的，则输出 “Orz, I cannot find x!”，注意逗号与 “I” 之间有一个空格。

样例

样例输入1

3 1
2 1 3
2 2 3
2 3 3

样例输出1

2
1
2

样例输入2

3 2
2 1 3
2 2 3
2 3 3

样例输出2

2
1
0

数据范围与提示

对于 100% 的数据，\(1\le y,z,p\le 10^9\)，为质数，\(1\le T\le 10\)。

1.快速幂取模。

2.拓展欧几里得：

　　\(xy≡ z ( mod\ p )\)

　　\(xy-bp=z;\)

　　设 \(t=gcd(y,p)\)，若 \(z\%t\ne 0\),则无解;

　　令 \(y/t，p/t，z/t\)，此时 \(gcd(y,p)=1\),

　　用 \(exgcd\) 求解 \(xy-bp=1\)

　得到 \(x\) 之后 \(\ast z\) 即为所求。

3.扩展BSGS算法，注意特殊情况处理

#include <iostream>
#include <cmath>
#include <unordered_map>
typedef long long LL;
int qpow(LL a, int b, const int &M) {
    int ans = 1;
    while (b) {
        if (b & 1) ans = ans * a % M;
        a = a * a % M, b >>= 1;
    }
    return ans;
}
int exgcd(int a, int b, int &x, int &y) {
    if (b == 0) {
        x = 1, y = 0;
        return a;
    }
    int gcd = exgcd(b, a%b, y, x);
    y -= a / b * x;
    return gcd;
}
void inv(int a, int b, int M) {
    int x, y, gcd = exgcd(a, M, x, y);
    if (b % gcd) { puts("Orz, I cannot find x!"); return; }
    b /= gcd, M /= gcd;
    b = (LL)b * x % M + M;
    printf("%d\n", b % M);
}
void bsgs(int a, int b, int M) {
    if (!(a%M) && b%M) { puts("Orz, I cannot find x!"); return; }
    if (!(a%M) && !(b%M)) { puts("1"); return; }
    std::unordered_map<LL, int> hash;
    LL t = ceil(sqrt(M)), at = 1;
    for (int i = 1; i <= t; ++i)
        at = at * a % M, hash[b*at%M] = i;
    a = at;
    for (int i = 1; i <= t; ++i) {
        LL v = at;
        if (hash[v]) {
            printf("%lld\n", ((i*t-hash[v])%M+M) % M);
            return;
        }
        at = at * a % M;
    }
    puts("Orz, I cannot find x!");
}
int main() {
    int T, K;
    scanf("%d%d", &T, &K);
    while (T--) {
        int y, z, p;
        scanf("%d%d%d", &y, &z, &p);
        if (K == 1) printf("%d\n", qpow(y, z, p));
        else if (K == 2) inv(y, z, p);
        else bsgs(y, z, p);
    }
    return 0;
}

【模板】扩展 BSGS/exBSGS

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <unordered_map>
using namespace std;
typedef long long LL;
const int INF = 0x3f3f3f3f;
int a, b, p;
unordered_map<int, int> hs;
int exgcd(int a, int b, int &x, int &y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }
    int d = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}
int BSGS(int a, int b, int p) {
    if (1 % p == b % p) return 0;
    int k = sqrt(p) + 1;
    hs.clear();
    for (int y = 0, r = b % p; y < k; y++) {
        hs[r] = y;
        r = (LL)r * a % p;
    }
    int ak = 1;
    for (int i = 1; i <= k; i++) ak = (LL)ak * a % p;
    for (int x = 1, l = ak; x <= k; x++) {
        if (hs.count(l)) return k * x - hs[l];
        l = (LL)l * ak % p;
    }
    return -INF;
}
int exBSGS(int a, int b, int p) {
    b = (b % p + p) % p;
    if (1 % p == b % p) return 0;
    int x, y;
    int d = exgcd(a, p, x, y);
    if (d > 1) {
        if (b % d) return -INF;
        exgcd(a / d, p / d, x, y);
        return exBSGS(a, (LL)b / d * x % (p / d), p / d) + 1;
    }
    return BSGS(a, b, p);
}
int main() {
    while (~scanf("%d%d%d", &a, &p, &b), a || b || p) {
        int res = exBSGS(a, b, p);
        if (res < 0) puts("No Solution");
        else printf("%d\n", res);
    }
    return 0;
}