高次同余方程

已知整数 $a,b,p$，其中 $p$ 为质数，求最小的非负整数 $x$，满足 $a^x \equiv b \pmod p$。

因为 $p$ 为质数，根据费马小定理，$a^{p-1} \equiv 1 \pmod p$。又知道 $a^0 \equiv 1 \pmod p$，所以幂有长度不超过 $p$ 的循环节。暴力枚举小于 $p$ 的数 $x$ 即可。

但是如果 $p$ 比较大怎么办呢？

BSGS Baby steps giant steps

$BSGS$ 算法的思想是 $meet \ in \ the \ middle$，时间复杂度 $\sqrt{p}$

先搜索前一半结果，把结果扔到哈希表中，然后搜索后一半，在表中检索

设 $x = im - j$，其中 $m = \left\lceil {\sqrt p } \right\rceil$，这样原式变为 $a^{im-j} \equiv b \pmod p$，进一步变形可得 $(a^m)^i \equiv b * a^j$

$j$ 的范围是 $[0 , m)$，枚举等号右边，把 $(ba^j,j)$ 加入哈希表。如果有重复的 $key$，用新的 $j$ 替代旧的

$i$ 的取值范围是 $[1,m]$，枚举等号左边，到哈希表中查找是否有相等的值。找到的第一个相等的 $i$，即可得出对应的 $im - j$ 是最小的 $x$

模版题：可爱的质数

可爱的质数

C++ 代码

#include <bits/stdc++.h>

using std::cin;
using std::cout;
using std::unordered_map;
using ll = long long;

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

void bsgs(ll a, ll b, ll p) {
    ll m = ceil(sqrt(p));
    ll t = b;
    unordered_map<ll, ll> dict;
    dict[b] = 0;
    for (int j = 1; j < m; ++j) {
        t = t * a % p;
        dict[t] = j;
    }
    ll mi = fastpow(a, m, p);
    t = 1;
    for (int i = 1; i <= m; ++i) {
        t = t * mi % p;
        if (dict.count(t) != 0) {
            cout << i * m - dict[t] << '\n';
            return;
        }
    }
    cout << "no solution\n";
}

int main() {
    ll p, b, n;
    cin >> p >> b >> n;

    bsgs(b, n, p);

    return 0;
}

exBSGS

扩展BSGS算法，针对 $p$ 不是质数的情况。

无解的判定条件是：$\gcd(a,p) \nmid b$ 且 $b \neq 1$

对于式子：

$$a^x \equiv b\pmod p, g = \gcd(a,p)$$

两边除掉 $g$

$$a^{x-1} \frac{a}{g} = \frac{b}{g} \pmod {\frac{p}{g}}$$

设 $a’ = \frac{a}{g}, b’ = \frac{b}{d’}, p’ = \frac{p}{g}$

原式可以变成

$$a’a^{x-1} = b’ \pmod p’$$

继续迭代，直到无解或者 $a,p$ 互质即可

// 求a^x=b(mod p)，其中p不一定是质数
ll exBSGS(ll a, ll b, ll p) {
    a = a % p;
    b = b % p;
    if (b == 1) return 0;
    ll cnt = 0, g, k = 1; // cnt是下面while循环迭代了几次
    while (true) {
        g = gcd(a, p);
        if (g == 1) break;
        if (b % g != 0) return -1;
        cnt++;
        k = k * (a / g) % p;
        b = b / g;
        p = p / g;
        if (k == b) return  cnt;
    } 
    ll m = ceil(sqrt(p));
    ll t = b;
    unordered_map<ll, ll> dict;
    dict[b] = 0;
    for (int j = 1; j < m; ++j) {
        t = t * a % p;
        dict[t] = j;
    }
    ll am = fastpow(a, m, p);
    t = k;
    for (int i = 1; i <= m; ++i) {
        t = t * am % p;
        if (dict.count(t) != 0) {
            return i * m - dict[t] + cnt;
        }
    }
    return -1;
} 

一些练习题

1. 计算器
2. 破解D-H协议
3. 多少个1？