Baby Step/Giant Step Algorithm\text{Baby Step/Giant Step Algorithm} 用于解决这样一类同余方程的方法

ax≡b(modc)a^x \equiv b\pmod c

结合原根的知识还能解决模数有原根的NN次剩余问题

xn≡b(modc)x^n \equiv b\pmod c

基础问题 ax≡b(modc)a^x\equiv b\pmod c

具体来说，BSGS的思路就是分块预处理+hash。

下面分情况讨论：

Case gcd(a,c) == 1

由欧拉定理可知xx的取值范围不会大于φ(c)\varphi(c)，我们直接取x<cx<c即可。

a,ca,c互素时，存在aa在模cc意义下的逆元。

令分块的大小为⌈c⌉\lceil\sqrt c\rceil，即我们将解的形式限定为x=A×⌈c⌉−B,A,B∈[0,⌈c⌉]x=A\times\lceil\sqrt c\rceil - B,A,B\in[0,\lceil\sqrt c\rceil]，直接将aBa^B的部分放到同余式右侧。

现在需要解决的问题变成：

aA×⌈c⌉≡b×aB(modc)a^{A\times\lceil\sqrt c\rceil}\equiv b\times a^{B}\pmod c

  • 首先预处理，枚举b×aBb\times a^B的取值，并将其存下来(hash)
  • 枚举AA的取值，并查看是否存在对应的BB是得两侧相等。
  • 得到结果，或者无解。

Case gcd(a,c) != 1

一个比较好的想法是，我们将其转换成互素的情况。

由于ax≡b(modc),g=gcd⁡(a,b,c)ax\equiv b\pmod c, g=\gcd(a, b, c)可以转化成 agx≡bg(modcg)\frac{a}{g}x\equiv \frac{b}{g}\pmod{\frac{c}{g}}

同样的，我们也可以不断去消去cc中与aa不互素的部分，从而达到目的

具体步骤为：

  • a×ax−1≡b(modc)a\times a^{x-1}\equiv b\pmod c
  • g=(a,c)g=(a,c)，若g∤bg\nmid b则无解
  • ag×ax−1≡bg(modcg)\frac{a}{g}\times a^{x-1}\equiv \frac{b}{g}\pmod{\frac{c}{g}}
  • 循环操作直到(a,c)=1(a,c)=1或akG=bG\frac{a^k}{G}=\frac{b}{G}

代码模板

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

ll qpow(ll x, ll n, ll mod) {
  ll res = 1LL;
  for (x %= mod; n > 0LL; n >>= 1, x = x * x % mod) {
    if (n & 1LL) res = res * x % mod;
  }
  return (res + mod) % mod;
}

// a^x EQUIV n (MOD mod), and gcd(a, mod) = 1
ll BSGS(ll a, ll n, ll mod) {
  a %= mod, n %= mod;
  if (n == 1LL || mod == 1LL) return 0LL;

  unordered_map<ll, ll> bs;
  ll S = sqrt(mod) + 1;

  ll base = n;
  for (ll k = 0, val = base; k <= S; k++) {
    bs[val] = k, val = val * a % mod;
  }
  base = qpow(a, S, mod);
  for (ll x = 1, val = base; x <= S; x++) {
    if (bs.count(val)) return x * S - bs[val];
    val = val * base % mod;
  }
  return -1;  // No solution
}

pair<ll, ll> exgcd(ll a, ll b) {
  bool neg_a = (a < 0), neg_b = (b < 0);
  ll x = 1, y = 0, r = 0, s = 1;
  while (b != 0LL) {
    ll t = a / b;
    r ^= x ^= r ^= x -= t * r;
    s ^= y ^= s ^= y -= t * s;
    b ^= a ^= b ^= a %= b;
  }
  return {neg_a ? -x : x, neg_b ? -y : y};
}

ll inv(ll a, ll mod) {
  auto [res, _] = exgcd(a, mod);
  return (res % mod + mod) % mod;
}

// a^x EQUIV n (MOD mod), and gcd(a, mod) != 1
ll exBSGS(ll a, ll n, ll mod) {
  a %= mod, n %= mod;
  if (n == 1LL || mod == 1LL) return 0LL;

  ll k = 0, val = 1;
  for (ll g = __gcd(a, mod); g != 1LL; g = __gcd(a, mod)) {
    if (n % g != 0LL) return -1;  // No solution
    mod /= g, n /= g;
    val = val * (a / g) % mod, k++;
    if (val == n) return k;
  }
  ll res = BSGS(a, n * inv(val, mod) % mod, mod);
  return ~res ? res + k : res;
}

NN次剩余问题（离散对数）

要求模数要有原根，令为gg。

  • 两侧取离散对数得到Nlog⁡x≡log⁡b(modφ(c))N\log x\equiv \log b\pmod{\varphi(c)}
  • 求出右侧log⁡b\log b，即求gy≡b(modc)g^y\equiv b\pmod c
  • 现在转换为一次同余式的求解

例题代码

51nod 1038 X^A mod P

题目链接

X^A mod P = B，其中P为质数。给出P和A B，求< P的所有X。

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111

#include <bits/stdc++.h>
using namespace std;

using ll = long long;
#define all(a) begin(a), end(a)

ll qpow(ll x, ll n, ll mod) {
  ll res = 1LL;
  for (x %= mod; n > 0LL; n >>= 1, x = x * x % mod) {
    if (n & 1LL) res = res * x % mod;
  }
  return (res + mod) % mod;
}

// a^x EQUIV n (MOD mod), and gcd(a, mod) = 1
ll BSGS(ll a, ll n, ll mod) {
  a %= mod, n %= mod;
  if (n == 1LL || mod == 1LL) return 0LL;

  unordered_map<ll, ll> bs;
  ll S = sqrt(mod) + 1;

  ll base = n;
  for (ll k = 0, val = base; k <= S; k++) {
    bs[val] = k, val = val * a % mod;
  }
  base = qpow(a, S, mod);
  for (ll x = 1, val = base; x <= S; x++) {
    if (bs.count(val)) return x * S - bs[val];
    val = val * base % mod;
  }
  return -1;  // No solution
}

tuple<ll, ll> exgcd(ll a, ll b) {
  ll x = 1, y = 0, r = 0, s = 1;
  while (b != 0LL) {
    ll t = a / b;
    r ^= x ^= r ^= x -= t * r;
    s ^= y ^= s ^= y -= t * s;
    b ^= a ^= b ^= a %= b;
  }
  return {a, x};
}

void solve() {
  ll a, b, p;
  cin >> p >> a >> b;

  ll phi = p - 1;  // p is a prime number

  b %= p, a %= phi;

  vector<int> fact = [](int n) {
    vector<int> res;
    for (int i = 2; i * i <= n; i++) {
      if (n % i == 0) {
        res.push_back(i);
        while (n % i == 0) n /= i;
      }
    }
    if (n != 1) res.push_back(n);
    return res;
  }(phi);

  ll root = [&fact, &phi](int p) {
    for (ll res = 2; res < p; res++) {
      bool ok = true;
      for (const ll &x : fact) {
        if (qpow(res, phi / x, p) != 1) continue;
        ok = false;
        break;
      }
      if (ok) return res;
    }
    return -1LL;
  }(p);

  ll L = BSGS(root, b, p);
  auto [g, inv] = exgcd(a, phi);
  /**
   * root^log(b) \equiv b (mod p)
   * 
   * a * log(ans) \equiv log(b) (mod phi)
   **/
  if (L % g != 0) return cout << "No Solution\n", void();

  a /= g, L /= g;
  ll mod = phi / g;

  L = (inv * L % mod + mod) % mod;

  vector<ll> ans;
  for (; L < phi; L += mod) {
    ans.push_back(qpow(root, L, p));
  }

  sort(all(ans));
  for (const ll &x : ans) cout << x << ' ';

  cout << '\n';
}

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);
  int T;
  cin >> T;

  while (T--) solve();
  return 0;
}