我们现在想解决的是这样一类问题：目标是在群 $\mathbb{F}^*_p$ 里面求离散对数 $\log_gh$. 众所周知，这是一个非常困难的问题；但如果有一些限定条件，问题就会简单很多。假设这里的 $p$ 选取不当，导致 $g$ 的阶，也就是 $p-1$ 有很多比较小的质因子：$$N = p-1 = q_1 ^ {e_1} \cdot q_2 ^ {e_2} \cdots q_t ^ {e_t},\qquad q_i\text{ is  small  primes}$$则会产生相应的漏洞。接下来，让我们一步步发明 Pohlig-Hellman 算法。

大步小步算法

　　大步小步算法(Shanks’s Babystep–Giantstep Algorithm)是通用的求解离散对数的方法，没有什么前置要求。它的思想是基于 meet-in-the-middle. 今假设在 $\mathbb{F} ^ *_p$ 有 $$g^x = h$$那我们选一个和 $\sqrt{\text{ord}(g)}$ 接近的数 $T$,则 $x$ 可以唯一地写成 $aT+b$，其中 $ a, b\leq T$. 则有：$$\begin{aligned}g ^ {aT + b} &= h \\ g ^ {aT} &= h \cdot g ^ {-b}\end{aligned}$$ 注意到等式左边的 $g^ {aT}$ 一共只有 $T$ 个，右边的 $g ^ {-b}$ 亦然。于是我们可以枚举 $a, b$，将所有可能的 $g ^ {aT}$ 和 $h\cdot g ^ {-b}$ 写进表，寻找碰撞。产生碰撞的那组 $a, b$ 为我们提供了答案 $aT + b$.

　　不难看出，大步小步算法的时间复杂度是 $O\left (\sqrt{\text{ord}(g)} \right )$. SageMath里面实现的离散对数就是基于大步小步算法和 Pohlig-Hellman 的，以下代码用于求 $\mathbb{F} ^ *_p$ 中 $\log_{627} 390$ 的结果。

R = GF(941)
h = R(390)
g = R(627)
x = discrete_log(h, g)  # 347
assert g**x == h

# 另一种方法
x = h.log(g)  # 347
assert g**x == h

　　值得指出，大步小步算法可以求解任意群内的离散对数，不仅限于 $GF(p)$.

分离与合并

　　回到最初的问题。给定 $(p, g)$, 我们是要在已知$$\text{ord}(g) = N = q_1 ^ {e_1} \cdot q_2 ^ {e_2} \cdots q_t ^ {e_t}$$的情况下求解离散对数 $g ^ x = h$. 如果直接采用大步小步算法，则复杂度是 $O(\sqrt n)$ 的。有没有更好的方法呢？假如我们有高效算法求出 $$x\bmod q_1 ^ {e_1},\quad x\bmod q_2 ^ {e_2},\quad \cdots,\quad x\bmod q_t ^ {e_t}$$ 那我们就可以利用CRT将它们合并起来！CRT是非常快的，所以只要我们能快速求出以上这些余数，就能很快地得到 $x\bmod N$，也就是原问题的结果。

　　现在设 $r_i = x \bmod q_i ^ {e_i}$，简记 $Q_i = q_i ^ {e_i}$，则有$$x = k_i Q_i + r_i$$现在来考虑如何才能找到 $r_i$. 由于 $g$ 的阶是 $N$，故有$g^N = 1$，从而$$\begin{aligned}h ^ {\frac{N}{Q_i}} &= (g ^ x) ^ {\frac{N}{Q_i}} \\ &= g ^ {k_i N} \cdot g ^ {r_i \frac{N}{Q_i}} \\ &=  \big (g ^ \frac{N}{Q_i}\big)^{r_i }\end{aligned}$$ 左边的$h ^ {\frac{N}{Q_i}}$ 是可以直接算出来的，右边的 $\big (g ^ \frac{N}{Q_i}\big)$ 是知道的，所以问题转化成了一个新的离散对数问题。只要解出这个问题，就可以得到 $r_i$. 但是这个算法的高明之处在哪里呢？它高明之处在于：$\big (g ^ \frac{N}{Q_i}\big)$ 的阶并不是 $p-1$，而是 $Q_i$，也就是 $q_i ^ {e_i}$，这比起原来的 $N$，降低了很多很多——从 $N$ 降到了质数的几次幂。

　　所以我们的算法模型很清楚了：解出上述 $t$ 个新的离散对数问题，得到 $r_1, r_2,\cdots,r_t$，最后利用 CRT 将结果合并起来。复杂度从原来的 $O(N)$ 降为 $$O\left(\sum _{i=1} ^ t \sqrt{q _i ^ {e _ i}} + \log N\right)$$

进一步优化

　　上面的式子有个很严重的问题：如果 $N$ 就是一个很小质数的某次方，那么等于什么优化也没有做。以$O(\sqrt{q _i ^ {e _ i}})$去解一个离散对数问题，仍然是很耗时的过程。不过我们有个好消息——需要解的离散对数$$\big (g ^ \frac{N}{Q_i}\big)^{r_i } = h ^ {\frac{N}{Q_i}} $$ 底数的阶是 $Q_i= q_i ^ {e _ i}$，也就是一个质数的幂。

　　有个算法能在 $O(e\sqrt{q})$ 的时间复杂度内，求出离散对数问题 $g ^ x = h$的解，其中必须满足 $g$ 的阶恰是质数的幂 $q^e$. 做法如下：先将 $x$ 写成 $q$ 进制的形式$$x = x_0 + x_1 q + x_2 q ^ 2 +\cdots + x _{e-1} q ^ {e-1} \quad \text{where } 0 \leq x_ i < q$$

　　来考虑 $h ^ {q ^ {e-1}}$. $$\begin{aligned} h ^ {q ^ {e-1}} &= (g ^ x) ^ {q ^ {e-1}} \\ &= g ^ {(x _ 0 + x _ 1 q + x _ 2 q ^ 2 + \cdots) q ^ {e-1}} \\ &= g ^ {x _ 0 q ^ {e-1}} \cdot g ^ {q ^ e (x_1 + x _ 2 q + \cdots)}\\&= g ^ {x _ 0 q ^ {e-1}}\\&= \left (g ^ {q ^ {e-1}}\right)^{x _ 0}\end{aligned}$$ 左边的 $h ^ {q ^ {e-1}}$ 和右边的 $\left (g ^ {q ^ {e-1}}\right)^{x _ 0}$ 都可以求出来，于是求 $x_0$ 又转化为了一个新的离散对数问题。但这里最好的消息是：底数 $\left (g ^ {q ^ {e-1}}\right)$ 的阶仅为 $q$. 我们在 $O(\sqrt q)$ 的时间复杂度内，就可以把 $x_0$ 求出来。

　　现在有 $x_0$ 了，如何求出 $x_1$？如法炮制，考虑 $h ^ {q ^ {e-2}}$.$$\begin{aligned} h ^ {q ^ {e-2}} &= (g ^ x) ^ {q ^ {e-2}} \\ &= g ^ {(x _ 0 + x _ 1 q + x _ 2 q ^ 2 + \cdots) q ^ {e-2}} \\ &= g ^ {x _ 0 q ^ {e-2}} \cdot  g ^ {x _ 1 q ^ {e-1}} \cdot g ^ {q ^ e (x_2 + x _ 3 q + \cdots)}\\&= g ^ {x _ 0 q ^ {e-2}} \cdot  g ^ {x _ 1 q ^ {e-1}}\\&= \left (g ^ {q ^ {e-1}}\right)^{x _ 1} \cdot g ^ {x _ 0 q ^ {e-2}}\end{aligned}$$于是有$$ \left (g ^ {q ^ {e-1}}\right)^{x _ 1} = \left( h \cdot g ^ {-x _ 0} \right)^ {q ^ {e-2}} $$又是一个底数的阶为 $q$ 的离散对数问题，于是 $x_2$ 也能在 $O(\sqrt{q})$ 复杂度内解出来。依次类推，我们可以解出所有的 $x_i$，时间复杂度是 $O(e \sqrt q)$.

　　到头来，最后解决整个问题的复杂度是：

$$O\left(\sum _{i=1} ^ t e_i\sqrt{q _i} + \log N\right)$$

代码实现

　　SageMath 里面的 discrete_log 就是利用 Pohlig-Hellman 和大步小步算法实现的。只需要调用 discrete_log 即可。 

R = GF(11251)
x = R(9689).log(23)
assert R(23) ^ x == 9689   # sage特色语法，这里 `^` 是幂运算