离散对数

定义

前置知识：阶与原根．

离散对数的定义方式和对数类似．取有原根的正整数模数 $m$，设其一个原根为 $g$. 对满足 $(a,m)=1$ 的整数 $a$，我们知道必存在唯一的整数 $0\leq k<\varphi(m)$ 使得

我们称这个 $k$ 为以 $g$ 为底，模 $m$ 的离散对数，记作 $k=\operatorname{ind}_g a$，在不引起混淆的情况下可记作 $\operatorname{ind} a$.

显然 $\operatorname{ind}_g 1=0$，$\operatorname{ind}_g g=1$.

性质

离散对数的性质也和对数有诸多类似之处．

大步小步算法

目前离散对数问题仍不存在多项式时间经典算法（离散对数问题的输入规模是输入数据的位数）．在密码学中，基于这一点人们设计了许多非对称加密算法，如 Ed25519．

在算法竞赛中，BSGS（baby-step giant-step，大步小步算法）常用于求解离散对数问题．形式化地说，对 $a,b,m\in\mathbf{Z}^+$，该算法可以在 $O(\sqrt{m})$ 的时间内求解

其中 $a\perp m$．方程的解 $x$ 满足 $0 \le x < m$.（注意 $m$ 不一定是素数）

算法描述

令 $x = A \left \lceil \sqrt m \right \rceil - B$，其中 $0\le A,B \le \left \lceil \sqrt m \right \rceil$，则有 $a^{A\left \lceil \sqrt m \right \rceil -B} \equiv b \pmod m$，稍加变换，则有 $a^{A\left \lceil \sqrt m \right \rceil} \equiv ba^B \pmod m$.

我们已知的是 $a,b$，所以我们可以先算出等式右边的 $ba^B$ 的所有取值，枚举 $B$，用 hash/map 存下来，然后逐一计算 $a^{A\left \lceil \sqrt m \right \rceil}$，枚举 $A$，寻找是否有与之相等的 $ba^B$，从而我们可以得到所有的 $x$，$x=A \left \lceil \sqrt m \right \rceil - B$.

注意到 $A,B$ 均小于 $\left \lceil \sqrt m \right \rceil$，所以时间复杂度为 $\Theta\left (\sqrt m\right )$，用 map 则多一个 $\log$.

扩展 BSGS 算法

对 $a,b,m\in\mathbf{Z}^+$，求解

其中 $a,m$ 不一定互质．

当 $(a, m)=1$ 时，在模 $m$ 意义下 $a$ 存在逆元，因此可以使用 BSGS 算法求解．于是我们想办法让他们变得互质．

具体地，设 $d_1=(a, m)$. 如果 $d_1\nmid b$，则原方程无解．否则我们把方程同时除以 $d_1$，得到

如果 $a$ 和 $\frac{m}{d_1}$ 仍不互质就再除，设 $d_2=\left(a, \frac{m}{d_1}\right)$. 如果 $d_2\nmid \frac{b}{d_1}$，则方程无解；否则同时除以 $d_2$ 得到

同理，这样不停的判断下去，直到 $a\perp \dfrac{m}{d_1d_2\cdots d_k}$.

记 $D=\prod_{i=1}^kd_i$，于是方程就变成了这样：

由于 $a\perp\dfrac{m}{D}$，于是推出 $\dfrac{a^k}{D}\perp \dfrac{m}{D}$. 这样 $\dfrac{a^k}{D}$ 就有逆元了，于是把它丢到方程右边，这就是一个普通的 BSGS 问题了，于是求解 $x-k$ 后再加上 $k$ 就是原方程的解啦．

注意，不排除解小于等于 $k$ 的情况，所以在消因子之前做一下 $\Theta(k)$ 枚举，直接验证 $a^i\equiv b \pmod m$，这样就能避免这种情况．

基于值域预处理的快速离散对数

前文的 BSGS 算法时间复杂度为单次 $O(\sqrt m)$，在询问量级较大的时候效率较低．若每次求解的模数是一个固定的质数 $p$，我们就有一个基于值域预处理的快速算法．

我们已经知道 $\operatorname{ind}_g(ab)\equiv\operatorname{ind}_g a+\operatorname{ind}_g b\pmod{p-1}$，所以我们可以只对所有质数通过 BSGS 算法计算离散对数，合数的离散对数则可通过该式转化为若干已知的质数离散对数值之和．此时复杂度仍然不优，我们考虑只预处理一部分的离散对数，具体来说，我们预处理 $1$ 到 $L = \lfloor\sqrt p\rfloor + 1$ 的离散对数．注意此时的 BSGS 块长 $B$ 不能取 $O(\sqrt{L})$，因为 BSGS 预处理（插入哈希表）部分的复杂度是 $O(B)$，而查询一共需要 $O(\pi(L))$ 次，则总时间复杂度为 $O\left(B+\dfrac{\pi(L)p}{B}\right)$，此时取 $B=O(\sqrt{\pi(L)p})$ 才是最优．由 素数定理 $\pi(n)\sim\dfrac{n}{\log n}$，则总的预处理时间复杂度可以平衡为 $O\left(\dfrac{p^{3/4}}{\log^{1/2} p}\right)$．

接下来是如何求答案．假设当前要求的是 $\operatorname{ind}_g y$，若 $y\le L$ 则直接返回，否则设 $p=vy+r$，则 $v=\left\lfloor\dfrac{p}{y}\right\rfloor<L$，$r=p\bmod y$，$y=\dfrac{p-r}{v}$，从而

注意到 $\operatorname{ind}_g (p-1)=(p-1)/2$，因此只需递归计算 $r$ 的离散对数即可．

我们还可以考虑 $y$ 的另一种表达方式，注意到 $p=vy+r=(v+1)y+r-y$，则 $y=\dfrac{p-r+y}{v+1}$，从而

我们有 $v+1 \le L$，因此只需要递归计算 $y-r$ 的离散对数即可．

综合这两种计算方式，我们有 $\min\{r,y-r\}\le \dfrac{y}{2}$，所以递归计算较小的一方即可达到 $O(\log p)$ 的查询复杂度．

至此我们得到了一个时间复杂度为 $O\left(\dfrac{p^{3/4}}{\log^{1/2} p}\right)-O(\log p)$ 的算法．

#include <cmath>
#include <iostream>
#include <unordered_map>
using namespace std;

const int N = 1'000'000;  // sqrt(P * sqrt(P) / ln(P))

unordered_map<int, int> M;
int Lg[N], p[N], t;
bool vis[N];
int P, g, B, sq, LP_1, inv;

int fadd(int x, int y, int P) {
  x += y;
  if (x >= P) x -= P;
  return x;
}

int fsub(int x, int y, int P) {
  x -= y;
  if (x < 0) x += P;
  return x;
}

int fmul(int x, int y, int P) { return 1ll * x * y % P; }

int qpow(int x, int y) {
  int ans = 1;
  while (y) {
    if (y & 1) ans = fmul(ans, x, P);
    x = fmul(x, x, P);
    y >>= 1;
  }
  return ans;
}

int calc(int x) {
  int s = x;
  for (int i = 0; i <= P / B; ++i) {
    if (M.find(s) != M.end()) return i * B + M[s];
    s = fmul(s, inv, P);
  }
  return -1;
}

void init() {
  int s = 1;
  for (int i = 0; i < B; ++i) {
    if (M.find(s) != M.end()) break;
    M[s] = i, s = fmul(s, g, P);
  }
  inv = qpow(qpow(g, B), P - 2);
  for (int i = 2; i <= sq; ++i) {
    if (!vis[i]) {
      p[++t] = i;
      Lg[i] = calc(i);
    }
    for (int j = 1; j <= t && p[j] * i <= sq; ++j) {
      vis[p[j] * i] = true;
      Lg[p[j] * i] = fadd(Lg[p[j]], Lg[i], P - 1);
      if (i % p[j] == 0) break;
    }
  }
}

int solve(int y) {
  if (y <= sq) return Lg[y];
  int v = P / y, r = P % y;
  if (r < y - r) return fadd(fsub(solve(r), Lg[v], P - 1), LP_1, P - 1);
  return fsub(solve(y - r), Lg[v + 1], P - 1);
}

int main() {
  ios::sync_with_stdio(false);
  cin.tie(nullptr);
  cin >> P >> g;
  sq = sqrt(P) + 1;
  B = sqrt(1ll * P * sqrt(P) / log(P));
  init();
  LP_1 = (P - 1) / 2;  // g ^ LP_1 = P - 1 (mod P)
  int T;
  cin >> T;
  while (T--) {
    int x;
    cin >> x;
    cout << solve(x) << '\n';
  }
  return 0;
}

习题

• SPOJ MOD 模板
• SDOI2013 随机数生成器
• SGU261 Discrete Roots 模板
• SDOI2011 计算器 模板
• Luogu4195【模板】exBSGS/Spoj3105 Mod 模板
• Codeforces - Lunar New Year and a Recursive Sequence
• LOJ6542 离散对数 index calculus 方法，非模板

本页面部分内容以及代码译自博文 Дискретное извлечение корня 与其英文翻译版 Discrete Root．其中俄文版版权协议为 Public Domain + Leave a Link；英文版版权协议为 CC-BY-SA 4.0．

参考资料

1. Discrete logarithm - Wikipedia
2. 潘承洞，潘承彪．初等数论．
3. 冯克勤．初等数论及其应用．