精读笔记：Cormen, Leiserson, Rivest, Stein《Introduction to Algorithms》（第 3 版，MIT Press 2009）｜负责 PDF 第 884–1093 页（第 29 章 29.2 节末尾 — 第 34 章 34.3 节中途）

说明：PDF 页码 = 原书页码 + 21（例如 PDF p.884 = 原书 p.863）。原书伪代码统一改写为 Python 风格伪代码（数组下标沿用原书的 1 起始约定，除非另行说明），并给出时间/空间复杂度。原书带 ★ 的节为研究生难度，笔记中照样标注。PDF 文本抽取的公式有缺失符号，笔记中已按上下文还原。

第 29 章线性规划（Linear Programming）（接上一块）

29.2 把问题表述为线性规划（末尾部分）（PDF p.884–885）

（接上一块）本块从多商品流问题（multicommodity-flow problem）的线性规划表述开始。

设定：有向图 $G=(V,E)$，每条边 $(u,v)$ 有非负容量 $c(u,v)\ge 0$；有 $k$ 种商品 $K_1,\dots,K_k$，商品 $i$ 用三元组 $K_i=(s_i,t_i,d_i)$ 描述：源点 $s_i$、汇点 $t_i$、需求量 $d_i$（希望从 $s_i$ 送到 $t_i$ 的流值）。
商品 $i$ 的流 $f_i$（$f_{iuv}$ 表示商品 $i$ 从 $u$ 到 $v$ 的流量）各自满足流守恒和容量约束；聚合流（aggregate flow）$f_{uv}=\sum_{i=1}^k f_{iuv}$ 不得超过边容量。
该问题没有目标函数，只问可行流是否存在，因此写成"空目标"（null objective）线性规划：

\begin{aligned} \min\quad & 0\\ \text{s.t.}\quad & \sum_{i=1}^k f_{iuv} \le c(u,v) && \forall u,v\in V\\ & \sum_{v\in V} f_{iuv}-\sum_{v\in V} f_{ivu}=0 && \forall i,\ \forall u\in V-\{s_i,t_i\}\\ & \sum_{v\in V} f_{i,s_i,v}-\sum_{v\in V} f_{i,v,s_i}=d_i && \forall i\\ & f_{iuv}\ge 0 && \forall u,v,\ \forall i \end{aligned}

要点：目前已知的唯一多项式时间算法就是把它写成线性规划，再用多项式时间 LP 算法（椭球法、内点法）求解。这说明 LP 作为"通用求解器"的价值：很多问题没有专门的组合算法，但能写成 LP。

29.2 习题概览（PDF p.884–885）：

29.2-1：把单对最短路 LP（式 29.44–29.46）化为标准型；
29.2-2：对图 24.2(a) 写出从 $s$ 到 $y$ 的最短路 LP；
29.2-3：单源最短路的 LP，使每个 $d_v$ 等于 $s$ 到 $v$ 的最短路权；
29.2-4：对图 26.1(a) 写最大流 LP；
29.2-5：把最大流 LP 改写为只用 $O(V+E)$ 个约束；
29.2-6：用 LP 求二分图最大匹配；
29.2-7：最小费用多商品流（minimum-cost multicommodity flow）：每条边另有费用 $a(u,v)$，目标是最小化 $\sum_{u,v}a(u,v)f_{uv}$，写成 LP。

29.3 单纯形算法（The simplex algorithm）（PDF p.885–900）

总体思路

单纯形法是解 LP 的经典方法。与本书大多数算法不同，它最坏情况下不是多项式时间的，但实践中通常非常快，并且对理解 LP 很有启发。
可以把它看作"不等式版的高斯消元"：高斯消元每步把方程组改写成结构更好的等价形式，直到解一目了然；单纯形法同样每轮把松弛型（slack form）改写为等价的松弛型。
每轮迭代对应一个基本解（basic solution）：令所有非基本变量（nonbasic variables，等式右边的变量）为 0，由等式算出基本变量（basic variables，等式左边的变量）的值。每次迭代后，基本可行解的目标值不减（通常严格增大）。
做法：选一个在目标函数中系数为正的非基本变量，把它从 0 往上增大，直到某个基本变量降为 0；然后交换二者角色，重写松弛型。算法并不显式维护这个解，只是不断改写 LP，直到最优解"一眼可见"。

完整例题（PDF p.886–890）

标准型 LP：

\begin{aligned} \max\quad & 3x_1+x_2+2x_3 &(29.53)\\ \text{s.t.}\quad & x_1+x_2+3x_3\le 30 &(29.54)\\ & 2x_1+2x_2+5x_3\le 24 &(29.55)\\ & 4x_1+x_2+2x_3\le 36 &(29.56)\\ & x_1,x_2,x_3\ge 0 &(29.57) \end{aligned}

引入松弛变量 $x_4,x_5,x_6$，得松弛型：

z=3x_1+x_2+2x_3,\quad x_4=30-x_1-x_2-3x_3,\quad x_5=24-2x_1-2x_2-5x_3,\quad x_6=36-4x_1-x_2-2x_3 .

术语：对某一组非基本变量取值，若某约束的基本变量恰为 0，称该约束是紧的（tight）；若会使基本变量为负，则违反该约束。松弛变量显式记录每个约束离"紧"还有多远。
3 个方程 6 个变量，有无穷多解；所有变量非负的解为可行解。基本解：右端（非基本）变量取 0，得 $(\bar x_1,\dots,\bar x_6)=(0,0,0,30,24,36)$，目标值 0。注意基本解满足 $\bar x_i=b_i\ (i\in B)$。基本解若可行，称为基本可行解（basic feasible solution）。单纯形法运行时基本解几乎总是可行的；29.5 节会看到开始几轮可能不可行。
改写不改变 LP 本身（可行解集合完全相同），只改变对应的基本解。

第 1 次转轴：增大 $x_1$。$x_4,x_5,x_6$ 分别在 $x_1>30, 12, 9$ 时变负；约束 (29.61) 最紧，故交换 $x_1$ 与 $x_6$。由 $x_6$ 的方程解出 $x_1=9-\tfrac{x_2}{4}-\tfrac{x_3}{2}-\tfrac{x_6}{4}$，代入其余方程：

\begin{aligned} z&=27+\tfrac{x_2}{4}+\tfrac{x_3}{2}-\tfrac{3x_6}{4}\\ x_1&=9-\tfrac{x_2}{4}-\tfrac{x_3}{2}-\tfrac{x_6}{4}\\ x_4&=21-\tfrac{3x_2}{4}-\tfrac{5x_3}{2}+\tfrac{x_6}{4}\\ x_5&=6-\tfrac{3x_2}{2}-4x_3+\tfrac{x_6}{2} \end{aligned}

这一操作称为转轴（pivot）：选择入基变量（entering variable）$x_e$ 和出基变量（leaving variable）$x_l$ 并交换其角色。两种操作（移项、代入）都保持等价（习题 29.3-3）。验证：原基本解 $(0,0,0,30,24,36)$ 仍满足新方程且目标值 $27+0+0-\tfrac34\cdot36=0$；新基本解 $(9,0,0,21,6,0)$，目标值 27。

第 2 次转轴：不能增大 $x_6$（会降低目标）。选 $x_3$：三个约束分别限制 $x_3\le 18$、$42/5$、$3/2$，第三个最紧，$x_3$ 入基、$x_5$ 出基：$x_3=\tfrac32-\tfrac{3x_2}{8}-\tfrac{x_5}{4}+\tfrac{x_6}{8}$，得

\begin{aligned} z&=\tfrac{111}{4}+\tfrac{x_2}{16}-\tfrac{x_5}{8}-\tfrac{11x_6}{16}\\ x_1&=\tfrac{33}{4}-\tfrac{x_2}{16}+\tfrac{x_5}{8}-\tfrac{5x_6}{16}\\ x_3&=\tfrac32-\tfrac{3x_2}{8}-\tfrac{x_5}{4}+\tfrac{x_6}{8}\\ x_4&=\tfrac{69}{4}+\tfrac{3x_2}{16}+\tfrac{5x_5}{8}-\tfrac{x_6}{16} \end{aligned}

基本解 $(33/4,0,3/2,69/4,0,0)$，目标值 $111/4$。

第 3 次转轴：只能增大 $x_2$。三个约束给出的上界为 132、4、$\infty$（$x_4$ 的方程中 $x_2$ 系数为正，增大 $x_2$ 反而使 $x_4$ 增大，不构成限制）。$x_2$ 增到 4，$x_2$ 入基、$x_3$ 出基：

\begin{aligned} z&=28-\tfrac{x_3}{6}-\tfrac{x_5}{6}-\tfrac{2x_6}{3}\\ x_1&=8+\tfrac{x_3}{6}+\tfrac{x_5}{6}-\tfrac{x_6}{3}\\ x_2&=4-\tfrac{8x_3}{3}-\tfrac{2x_5}{3}+\tfrac{x_6}{3}\\ x_4&=18-\tfrac{x_3}{2}+\tfrac{x_5}{2} \end{aligned}

目标函数所有系数都为负，此时基本解即最优解：$(8,4,0,18,0,0)$，目标值 28。回到原问题：$x_1=8,x_2=4,x_3=0$，$3\cdot8+4=28$。

松弛变量的含义：$x_4=18$ 表示约束 (29.54) 左边 $8+4+0=12$ 比右边 30 少 18；$x_5=x_6=0$ 表示 (29.55)(29.56) 取等号（紧约束）。
常见误区：即使初始系数全为整数，中间松弛型的系数和中间解都可能是分数；最终最优解也未必是整数——本例恰好是整数纯属巧合（这正是整数规划困难的根源，见 34 章习题）。

转轴过程 PIVOT（PDF p.890–891）

输入松弛型 $(N,B,A,b,c,v)$、出基下标 $l$、入基下标 $e$；输出新松弛型 $(\hat N,\hat B,\hat A,\hat b,\hat c,\hat v)$。约定：松弛型写作 $z=v+\sum_{j\in N}c_jx_j$，$x_i=b_i-\sum_{j\in N}a_{ij}x_j\ (i\in B)$，即矩阵 $A$ 的元素是松弛型中系数的相反数。

def PIVOT(N, B, A, b, c, v, l, e):
    # 1) 新基本变量 x_e 的方程：把 x_l 那一行解出 x_e
    A_hat = new_matrix(m, n)
    b_hat[e] = b[l] / A[l][e]
    for j in N - {e}:
        A_hat[e][j] = A[l][j] / A[l][e]
    A_hat[e][l] = 1 / A[l][e]
    # 2) 其余约束：把 x_e 的新表达式代入
    for i in B - {l}:
        b_hat[i] = b[i] - A[i][e] * b_hat[e]
        for j in N - {e}:
            A_hat[i][j] = A[i][j] - A[i][e] * A_hat[e][j]
        A_hat[i][l] = -A[i][e] * A_hat[e][l]
    # 3) 目标函数同样代入
    v_hat = v + c[e] * b_hat[e]
    for j in N - {e}:
        c_hat[j] = c[j] - c[e] * A_hat[e][j]
    c_hat[l] = -c[e] * A_hat[e][l]
    # 4) 更新基/非基集合
    N_hat = (N - {e}) | {l}
    B_hat = (B - {l}) | {e}
    return N_hat, B_hat, A_hat, b_hat, c_hat, v_hat

时间 $O(mn)$（双重循环），空间 $O(mn)$（新矩阵；可原地更新）。若 $a_{le}=0$ 会除零，但后面证明只在 $a_{le}\ne0$（实际是 $a_{le}>0$）时调用。

引理 29.1：在 $a_{le}\ne0$ 时调用 PIVOT，设 $\bar x$ 为调用后的基本解，则 (1) $\bar x_j=0,\ j\in\hat N$；(2) $\bar x_e=b_l/a_{le}$；(3) $\bar x_i=b_i-a_{ie}\hat b_e,\ i\in\hat B-\{e\}$。证明：基本解令非基本变量为 0，故 $\bar x_i=\hat b_i$，再分别对照 PIVOT 第 3 行和第 9 行即得。

形式化的单纯形算法（PDF p.891–893）

需要解决的问题：如何判断 LP 可行？LP 可行但初始基本解不可行怎么办？如何判断无界？如何选入基/出基变量？前两个问题交给 29.5 节的 INITIALIZE-SIMPLEX$(A,b,c)$：输入标准型（$m\times n$ 矩阵 $A$、$m$ 维 $b$、$n$ 维 $c$），若不可行则报告并终止，否则返回一个基本解可行的松弛型。

def SIMPLEX(A, b, c):
    N, B, A, b, c, v = INITIALIZE_SIMPLEX(A, b, c)
    delta = [None] * (m + n + 1)            # 每个约束允许 x_e 增加的上限
    while any(c[j] > 0 for j in N):
        e = 选取某个 j in N 且 c[j] > 0       # 按某固定确定性规则
        for i in B:
            delta[i] = b[i] / A[i][e] if A[i][e] > 0 else INF
        l = argmin_{i in B} delta[i]          # 最小比值检验（ratio test）
        if delta[l] == INF:
            return "unbounded"
        N, B, A, b, c, v = PIVOT(N, B, A, b, c, v, l, e)
    x = [b[i] if i in B else 0 for i in 1..n]
    return x

循环在目标函数全部系数 $\le 0$ 时终止（原文说"都为负"，实质是"没有正系数"）。第 5–9 行找对 $x_e$ 增大限制最严的约束，其基本变量为出基变量 $x_l$；若没有约束限制 $x_e$，返回"无界"。
每次迭代代价 $O(mn)$（PIVOT 主导），迭代次数最坏 $\binom{n+m}{m}$（可为指数级）。空间 $O(mn)$。

引理 29.2（可行性/无界性的正确性）：若 INITIALIZE-SIMPLEX 返回的松弛型基本解可行，则 SIMPLEX 若在第 17 行返回解，该解可行；若返回 unbounded，则 LP 确实无界。

证明用三部分循环不变式：每次 while 迭代开始时，(1) 松弛型与初始松弛型等价；(2) 对所有 $i\in B$，$b_i\ge0$；(3) 对应基本解可行。

初始化：显然。
保持：(1) 由习题 29.3-3；(2) 关键：$\hat b_e=b_l/a_{le}\ge0$（因 $b_l\ge0$、$a_{le}>0$）。对其余 $i$，$\hat b_i=b_i-a_{ie}(b_l/a_{le})$：若 $a_{ie}>0$，由最小比值选法 $b_l/a_{le}\le b_i/a_{ie}$，得 $\hat b_i\ge b_i-a_{ie}(b_i/a_{ie})=0$；若 $a_{ie}\le0$，则 $\hat b_i\ge b_i\ge0$。(3) 基本解 $\bar x_i=b_i\ge0$，非基本为 0，故可行。
终止：若因第 3 行条件退出，返回可行解。若返回 unbounded，则对所有 $i\in B$ 有 $a_{ie}\le0$。构造 $\bar x_e=\infty$（理解为任意大的 $t$），其余非基本变量 0，$\bar x_i=b_i-a_{ie}\bar x_e\ge0$ 可行；目标值 $v+c_e\bar x_e\to\infty$（因 $c_e>0$），故 LP 无界。

终止性、退化与循环（PDF p.895–899）

习题 29.3-2：每次迭代不会降低目标值。但可能不变，称为退化（degeneracy）。由 PIVOT 第 14 行 $\hat v=v+c_e\hat b_e$，$c_e>0$，目标值不变当且仅当 $\hat b_e=b_l/a_{le}=0$，即 $b_l=0$。
退化例：$z=x_1+x_2+x_3$，$x_4=8-x_1-x_2$，$x_5=x_2-x_3$。$x_1$ 入、$x_4$ 出后得 $z=8+x_3-x_4$，$x_1=8-x_2-x_4$，$x_5=x_2-x_3$；此时只能 $x_3$ 入、$x_5$ 出，由于 $b_5=0$，目标仍为 8：$z=8+x_2-x_4-x_5$，$x_3=x_2-x_5$。再以 $x_2$ 入、$x_1$ 出，目标升到 16。
循环（cycling）：退化可能导致两次迭代的松弛型完全相同；SIMPLEX 是确定性算法，一旦循环就永远不终止。循环是唯一可能的不终止原因。

引理 29.3（代数引理）：若对一切 $x_j$ 取值都有 $\sum_{j\in I}\alpha_jx_j=\gamma+\sum_{j\in I}\beta_jx_j$，则 $\alpha_j=\beta_j$ 且 $\gamma=0$。证明：令所有 $x_j=0$ 得 $\gamma=0$；令某个 $x_j=1$、其余为 0 得 $\alpha_j=\beta_j$。

引理 29.4：给定基本变量集 $B$，对应的松弛型唯一确定。证明：设有两个同 $B$ 的松弛型，相减得 $\sum_{j\in N}a_{ij}x_j=(b_i-b'_i)+\sum_{j\in N}a'_{ij}x_j$，对每个 $i$ 应用引理 29.3 得 $a_{ij}=a'_{ij}$、$b_i=b'_i$；目标函数同理（习题 29.3-1）。

引理 29.5：若 SIMPLEX 在 $\binom{n+m}{m}$ 次迭代内不终止，则它循环。证明：$n+m$ 个变量中选 $m$ 个作基，至多 $\binom{n+m}{m}$ 种，每种对应唯一松弛型。

防止循环：(a) 对输入做微小扰动（perturbation），使不同解的目标值不同；(b) Bland 规则（Bland's rule）：入基、出基都在并列时选下标最小的变量。证明从略。

引理 29.6：若第 4 行和第 9 行都按最小下标打破平局，SIMPLEX 必终止。

引理 29.7：若 INITIALIZE-SIMPLEX 返回基本解可行的松弛型，SIMPLEX 要么报告无界，要么在至多 $\binom{n+m}{m}$ 次迭代内以可行解终止。（由引理 29.2、29.6 和 29.5 的逆否命题。）

29.3 习题概览（PDF p.899–900）：29.3-1 补全引理 29.4（$c=c'$、$v=v'$）；29.3-2 PIVOT 不降低 $v$；29.3-3 PIVOT 前后松弛型等价；29.3-4 标准型化为松弛型后，基本解可行当且仅当所有 $b_i\ge0$；29.3-5/6/7 用 SIMPLEX 手算三个具体 LP（如 $\max 18x_1+12.5x_2$，$x_1+x_2\le20,x_1\le12,x_2\le16$；以及一个最小化问题 $\min x_1+x_2+x_3$，$2x_1+7.5x_2+3x_3\ge10000$，$20x_1+5x_2+10x_3\ge30000$）；29.3-8 举例说明基的选法可严格少于 $\binom{m+n}{n}$。

29.4 对偶性（Duality）（PDF p.900–907）

动机

前面证明了 SIMPLEX 会终止，但没证明终止时的解最优。为此引入线性规划对偶（linear-programming duality）——它是证明最优性的工具。
先例：第 26 章最大流最小割定理。给定流 $f$，若能找到容量等于 $|f|$ 的割，就证明 $f$ 是最大流。对偶的一般含义：给一个最大化问题，定义一个相关的最小化问题，使两者最优值相同。原问题称为原始问题（primal）。

对偶 LP 的定义

标准型原问题 $\max\ c^Tx$，s.t. $Ax\le b$，$x\ge0$（式 29.16–29.18）的对偶为：

\min\ \sum_{i=1}^m b_iy_i\quad\text{s.t.}\quad \sum_{i=1}^m a_{ij}y_i\ge c_j\ (j=1,\dots,n),\qquad y_i\ge0\ (i=1,\dots,m)\qquad(29.83\text{–}29.85)

构造规则：max 变 min；右端常数与目标系数互换角色；$\le$ 变 $\ge$。原问题每个约束对应一个对偶变量 $y_i$，对偶每个约束对应一个原变量 $x_j$。

例：29.3 节例题（29.53–29.57）的对偶为 $\min 30y_1+24y_2+36y_3$，s.t. $y_1+2y_2+4y_3\ge3$，$y_1+2y_2+y_3\ge1$，$3y_1+5y_2+2y_3\ge2$，$y\ge0$。

弱对偶

引理 29.8（弱对偶，weak duality）：$\bar x$ 原问题可行、$\bar y$ 对偶可行，则 $\sum_j c_j\bar x_j\le\sum_i b_i\bar y_i$。

证明（两步放缩，一次交换求和次序）：

\sum_{j}c_j\bar x_j\le\sum_j\Big(\sum_i a_{ij}\bar y_i\Big)\bar x_j=\sum_i\Big(\sum_j a_{ij}\bar x_j\Big)\bar y_i\le\sum_i b_i\bar y_i .

第一步用对偶约束且 $\bar x_j\ge0$，第二步用原约束且 $\bar y_i\ge0$。

推论 29.9：若两可行解目标值相等，则二者分别是原、对偶的最优解（谁也无法再改进）。

从单纯形最终松弛型读出对偶最优解

例题最终松弛型 $z=28-x_3/6-x_5/6-2x_6/3$，$B=\{1,2,4\}$，$N=\{3,5,6\}$。规则：对偶变量 = 最终目标函数中对应松弛变量系数的相反数：

\bar y_i=\begin{cases}-c'_{n+i}, & (n+i)\in N\\ 0, & \text{否则}\end{cases}\qquad(29.91)

例中 $\bar y_1=0$（$x_4\in B$），$\bar y_2=-c'_5=1/6$，$\bar y_3=-c'_6=2/3$；对偶目标 $30\cdot0+24\cdot\tfrac16+36\cdot\tfrac23=28$，与原问题相等，于是由引理 29.8 证明了最优值为 28。（这就是经济学里"影子价格"（shadow price）的来源：$y_i$ 表示第 $i$ 种资源多一单位可增加的最优目标值。）

强对偶

定理 29.10（线性规划对偶定理）：SIMPLEX 对原问题 $(A,b,c)$ 返回 $\bar x$，$N,B$ 为最终松弛型的非基/基变量集合，$c'$ 为最终目标系数，$\bar y$ 由 (29.91) 定义。则 $\bar x$ 原问题最优、$\bar y$ 对偶最优，且 $\sum_j c_j\bar x_j=\sum_i b_i\bar y_i$。

证明思路：

最终目标 $z=v'+\sum_{j\in N}c'_jx_j$，终止条件给出 $c'_j\le0\ (j\in N)$；补定义 $c'_j=0\ (j\in B)$，则 $z=v'+\sum_{j=1}^{n+m}c'_jx_j$。基本解下非基变量为 0、基变量系数为 0，故原目标值 $\sum_j c_j\bar x_j=v'$。
所有松弛型等价，故对任意 $x$，$\sum_{j=1}^n c_jx_j=v'+\sum_{j=1}^{n+m}c'_jx_j$。把松弛变量 $x_{n+i}=b_i-\sum_j a_{ij}x_j$ 代入并用 $c'_{n+i}=-\bar y_i$，整理得 $\sum_{j=1}^n c_jx_j=\Big(v'-\sum_i b_i\bar y_i\Big)+\sum_{j=1}^n\Big(c'_j+\sum_i a_{ij}\bar y_i\Big)x_j .\qquad(29.99)$
对 (29.99) 应用引理 29.3：$v'-\sum_ib_i\bar y_i=0$（29.100），且 $c'_j+\sum_ia_{ij}\bar y_i=c_j$（29.101）。前者说明对偶目标值等于原目标值 $v'$。
对偶可行性：$c'_j\le0$ 对所有 $j$ 成立，故 $c_j=c'_j+\sum_ia_{ij}\bar y_i\le\sum_ia_{ij}\bar y_i$；且 $\bar y_i=-c'_{n+i}\ge0$。由推论 29.9 得最优。

结论：若 LP 可行、INITIALIZE-SIMPLEX 返回可行松弛型、SIMPLEX 不报告无界，则返回的解最优，并同时得到对偶最优解。

29.4 习题概览：29.4-1 写 29.3-5 的对偶；29.4-2 直接对非标准型 LP 写对偶的规则（等式约束→自由对偶变量，自由变量→等式对偶约束等）；29.4-3 最大流 LP 的对偶并解释为最小割；29.4-4 最小费用流 LP 的对偶及其图论解释；29.4-5 对偶的对偶是原问题；29.4-6 第 26 章哪个结论可视为最大流的弱对偶（任一流值 ≤ 任一割容量，引理 26.4 推论）。

29.5 初始基本可行解（The initial basic feasible solution）（PDF p.907–915）

问题

LP 可行，但初始基本解（所有原变量取 0）不一定可行。例： $\max 2x_1-x_2$，s.t. $2x_1-x_2\le2$，$x_1-5x_2\le-4$，$x_1,x_2\ge0$。取 $x_1=x_2=0$ 违反第二个约束。解决办法：构造辅助线性规划（auxiliary linear program）。

引理 29.11：设 $L$ 为标准型 LP，引入新变量 $x_0$，定义 $L_{aux}$（$n+1$ 个变量）：

\max\ -x_0\quad\text{s.t.}\quad \sum_{j=1}^n a_{ij}x_j-x_0\le b_i\ (i=1..m),\qquad x_j\ge0\ (j=0..n).

则 $L$ 可行 $\iff$ $L_{aux}$ 的最优值为 0。证明：$L$ 可行解配 $x_0=0$ 是 $L_{aux}$ 的可行解，目标 0，而 $-x_0\le0$，故最优；反之最优值 0 意味着 $\bar x_0=0$，其余分量满足 $L$。（直观：$x_0$ 是"允许违反约束的量"，最小化它。）

INITIALIZE-SIMPLEX

def INITIALIZE_SIMPLEX(A, b, c):
    k = argmin_i b[i]
    if b[k] >= 0:                      # 初始基本解已可行
        return ({1..n}, {n+1..n+m}, A, b, c, 0)
    # 构造 L_aux：每个约束左边加 -x0，目标改为 -x0
    N, B, A, b, c, v = slack_form_of(L_aux)   # 非基 {0,1..n}，基 {n+1..n+m}
    l = n + k                          # 基本解中最负的那个基变量出基
    N, B, A, b, c, v = PIVOT(N, B, A, b, c, v, l, 0)   # x0 入基，之后基本解可行
    反复执行 SIMPLEX 的 while 循环（第 3–12 行），直到求得 L_aux 的最优解
    if 最优解中 x0 == 0:
        if x0 是基变量:
            做一次（退化）转轴，选任一 e∈N 且 a[0][e] != 0 入基，使 x0 出基
        从约束中删去 x0；恢复 L 的原目标函数，
        并把其中出现的每个基变量替换为其约束的右端表达式
        return 修改后的最终松弛型
    else:
        return "infeasible"

第 1–3 行隐式检查初始松弛型（标准型与松弛型的 $A,b,c$ 相同，不需显式转换）。
关键技巧：仅一次转轴（$x_0$ 入基、最负的 $x_{n+k}$ 出基）就能让 $L_{aux}$ 的基本解可行。
退化情形：$x_0$ 仍是基变量但取值 0，做一次退化转轴把它移出基，不改变任何变量的值。
复杂度：外加一次完整的单纯形求解（同样最坏指数次迭代），每次转轴 $O(mn)$。

例题（PDF p.909–911）

辅助问题 $\max -x_0$，s.t. $2x_1-x_2-x_0\le2$，$x_1-5x_2-x_0\le-4$。松弛型：$z=-x_0$，$x_3=2-2x_1+x_2+x_0$，$x_4=-4-x_1+5x_2+x_0$。基本解 $x_4=-4$ 不可行。

$x_0$ 入基、$x_4$（最负）出基：$z=-4-x_1+5x_2-x_4$，$x_0=4+x_1-5x_2+x_4$，$x_3=6-x_1-4x_2+x_4$；基本解 $(x_0,\dots,x_4)=(4,0,0,6,0)$ 可行。
$x_2$ 入基、$x_0$ 出基：$z=-x_0$，$x_2=\tfrac45-\tfrac{x_0}{5}+\tfrac{x_1}{5}+\tfrac{x_4}{5}$，$x_3=\tfrac{14}{5}+\tfrac{4x_0}{5}-\tfrac{9x_1}{5}+\tfrac{x_4}{5}$。最优值 0，原问题可行。
删去 $x_0$，恢复原目标 $2x_1-x_2$ 并代入 $x_2$ 的表达式：$z=-\tfrac45+\tfrac{9x_1}{5}-\tfrac{x_4}{5}$，约束 $x_2=\tfrac45+\tfrac{x_1}{5}+\tfrac{x_4}{5}$，$x_3=\tfrac{14}{5}-\tfrac{9x_1}{5}+\tfrac{x_4}{5}$。这个松弛型基本解可行，交给 SIMPLEX。

引理 29.12：若 $L$ 不可行，INITIALIZE-SIMPLEX 返回 "infeasible"；否则返回一个基本解可行的合法松弛型。证明要点：

不可行时：$L_{aux}$ 最优值非零，又 $x_0\ge0$ 故为负；且有限（$x_i=0$、$x_0=|\min_i b_i|$ 可行，目标 $-|\min b_i|$），第 11 行检验失败返回 infeasible。
可行且 $b\ge0$：直接返回（习题 29.3-4）。
可行但 $b_k<0$：在 $L_{aux}$ 的标准型中 $x_0$ 系数 $a_{i0}=-1$（对所有 $i$），故 $a_{le}=-1$。转轴后 $\bar x_0=b_l/a_{le}=-b_k>0$；其余 $\bar x_i=b_i-a_{ie}(b_l/a_{le})=b_i-b_l\ge0$（因 $b_l=b_k$ 是最小值）。所以基本解可行；随后解 $L_{aux}$ 得最优值 0，删去 $x_0$ 即得 $L$ 的可行松弛型。

线性规划基本定理

定理 29.13（Fundamental theorem of linear programming）：标准型 LP $L$ 必为以下三者之一：(1) 有有限最优值的最优解；(2) 不可行；(3) 无界。SIMPLEX 在三种情形下分别返回最优解、"infeasible"、"unbounded"。证明把引理 29.12、29.7、定理 29.10、引理 29.2 串起来即可。

29.5 习题概览：29.5-1 写出 INITIALIZE-SIMPLEX 第 5、14 行的详细伪代码；29.5-2 证明在 INITIALIZE-SIMPLEX 内运行主循环不会返回 unbounded（$-x_0\le0$ 有上界）；29.5-3 若 $L$ 与其对偶的初始基本解都可行，则 $L$ 最优值为 0；29.5-4 允许严格不等式时基本定理不成立（例如 $\max x$ s.t. $x<1$ 无最优解）；29.5-5～29.5-8 手算若干 LP（需要两阶段）；29.5-9 单变量 LP $P:\max tx$ s.t. $rx\le s$，$x\ge0$ 与其对偶 $D$，讨论 $r,s,t$ 何时出现"二者都有限最优""P 可行 D 不可行""D 可行 P 不可行""都不可行"四种情况。

第 29 章思考题与章末注记（PDF p.915–918）

思考题：

29-1 线性不等式可行性（linear-inequality feasibility）：(a) 用 LP 算法解可行性问题；(b) 反过来用可行性算法解 LP（提示：把原约束、对偶约束和"原目标 ≥ 对偶目标"放在一起，强对偶保证可行解即最优解）。规模都要是 $n,m$ 的多项式。
29-2 互补松弛性（complementary slackness）：$\bar x,\bar y$ 分别原/对偶可行，则二者同为最优的充要条件是 $\sum_i a_{ij}\bar y_i=c_j\ \text{或}\ \bar x_j=0\ (\forall j);\qquad \sum_j a_{ij}\bar x_j=b_i\ \text{或}\ \bar y_i=0\ (\forall i).$ (a) 在 29.53 例上验证；(b) 一般证明；(c) 原可行解 $\bar x$ 最优 ⟺ 存在对偶可行 $\bar y$，使 $\bar x_j>0$ 处对偶约束取等，原约束有松弛（$<b_i$）处 $\bar y_i=0$。
29-3 整数线性规划（integer linear programming）：判定可行性就是 NP 难的（习题 34.5-3）。(a) 弱对偶仍成立；(b) 强对偶不一定成立；(c) 证明 $IP\le P=D\le ID$（整数间隙，integrality gap）。
29-4 Farkas 引理：$A$ 为 $m\times n$ 矩阵，$c$ 为 $n$ 维向量，则以下两个系统恰有一个有解：$Ax\le0,\ c^Tx>0$；与 $A^Ty=c,\ y\ge0$。
29-5 最小费用循环流（minimum-cost circulation）：没有源汇和需求，只要容量约束和处处守恒，求费用最小的可行流。(a) 写 LP；(b) 所有边费用为正时最优解为零流；(c) 把最大流化为最小费用循环流（加一条 $t\to s$ 的大容量负费用边）；(d) 把单源最短路化为最小费用循环流。

章末注记：单纯形法由 G. Dantzig 于 1947 年提出，至今其变体仍是最常用的 LP 方法。第一个多项式时间算法是 Khachian（1979）的椭球法（ellipsoid algorithm），但实践中不敌单纯形法；Karmarkar 提出第一个内点法（interior-point algorithm）。Klee 与 Minty 构造了让单纯形法走 $2^n-1$ 次迭代的例子；Borgwardt 等证明了在某些概率假设下期望多项式时间；Spielman 与 Teng 提出平滑分析（smoothed analysis）解释其实际高效。网络单纯形法（network simplex）在最短路、最大流、最小费用流等问题上有多项式时间变体。参考书：Chvátal、Gass、Karloff、Schrijver、Vanderbei 等；本章讲法取自 Chvátal。

第 29 章本章要点（就本块覆盖部分）

LP 的标准型 $\max c^Tx,\ Ax\le b,\ x\ge0$ 与松弛型（引入松弛变量，基变量用非基变量表示）是单纯形法的工作语言；基本解 = 非基变量取 0。
单纯形法每轮做一次转轴：选目标系数为正的入基变量，用最小比值检验选出基变量；每次转轴 $O(mn)$，迭代次数最坏指数级。
退化可能导致循环；Bland 规则（最小下标）可防止循环。基变量集合唯一决定松弛型，故不循环则至多 $\binom{n+m}{m}$ 步终止。
对偶：弱对偶 $c^T\bar x\le b^T\bar y$；强对偶由单纯形法构造性证明，对偶最优解就是最终目标中松弛变量系数的相反数（影子价格）。
两阶段思想：用辅助 LP（最小化人工变量 $x_0$）判定可行性并获得初始基本可行解。
基本定理：LP 只有"有限最优""不可行""无界"三种结局。互补松弛、Farkas 引理是对偶理论的重要推论。

第 29 章与量化交易的关联

组合优化：带线性约束和线性目标的组合构建问题（如最大化预期 alpha，约束行业/风格暴露、仓位上下限、换手）就是 LP；风险项若用 L1 或 CVaR（Rockafellar–Uryasev 形式）度量，也可写成 LP。实务中用 HiGHS、Gurobi、CPLEX 等求解器，但理解基、退化、对偶对解读求解器输出很关键。
影子价格与约束归因：对偶变量告诉你每条约束"放松一单位能多挣多少目标值"，可用来评估行业中性、换手上限、单票上限等约束对组合 alpha 的成本，这是组合约束归因（constraint attribution）的基础。
无套利与 Farkas 引理：资产定价基本定理的离散版本——"不存在套利 ⟺ 存在正的状态价格向量"——正是 Farkas 引理/LP 对偶的直接应用；也可用 LP 检验一组期权报价是否存在静态套利、求无模型的价格上下界（超复制问题的对偶）。
执行与交易调度：多商品流、最小费用流可用于多资产、多场所的订单拆分与路由建模；最小费用流也用于头寸转移/资金调拨。
整数规划的困难：手数、最小交易单位、持仓数目上限（基数约束）使问题变为 MILP（NP 难），29-3 中的整数间隙解释了为何"先松弛再取整"可能明显次优。

第 29 章推荐习题

29.3-5、29.5-5：手算单纯形与两阶段法，熟悉转轴和比值检验。
29.4-2、29.4-5：掌握任意形式 LP 的对偶写法——组合优化里解读对偶最常用。
29.4-3：最大流对偶即最小割，理解对偶的组合含义。
思考题 29-2（互补松弛）和 29-4（Farkas 引理）：对偶理论核心，直接连接到无套利定价。
思考题 29-3：整数规划的对偶间隙。

第 30 章多项式与快速傅里叶变换（Polynomials and the FFT）（PDF p.919–940）

第 30 章引言（PDF p.919–921）

两个 $n$ 次多项式相加用直接方法需 $\Theta(n)$，相乘需 $\Theta(n^2)$；快速傅里叶变换（fast Fourier transform, FFT）把乘法降为 $\Theta(n\lg n)$。
FFT 最常见的用途是信号处理：信号在时域（time domain）给出（时间→幅度）；傅里叶分析把它表示为不同频率、带相移的正弦波加权和，权重和相位刻画了信号的频域（frequency domain）。应用包括 MP3 等音视频压缩。
多项式：域 $F$（通常是复数 $\mathbb C$）上的形式和 $A(x)=\sum_{j=0}^{n-1}a_jx^j$，$a_j$ 为系数（coefficients）。最高非零系数为 $a_k$ 时次数（degree）为 $k$；任何严格大于次数的整数都是次数界（degree-bound），次数界为 $n$ 的多项式次数可为 $0..n-1$。
加法：$c_j=a_j+b_j$。例：$A=6x^3+7x^2-10x+9$，$B=-2x^3+4x-5$，$C=4x^3+7x^2-6x+4$。
乘法：$C(x)=A(x)B(x)$ 的次数界为 $2n-1$： $C(x)=\sum_{j=0}^{2n-2}c_jx^j,\qquad c_j=\sum_{k=0}^{j}a_kb_{j-k}.\qquad(30.1,30.2)$ 例：上面 $A\cdot B=-12x^6-14x^5+44x^4-20x^3-75x^2+86x-45$。
$\deg C=\deg A+\deg B$；次数界 $n_a$ 与 $n_b$ 的乘积次数界为 $n_a+n_b-1$，通常就说 $n_a+n_b$。
本章用 $i$ 专指 $\sqrt{-1}$。章节安排：30.1 两种表示（系数表示、点值表示）；30.2 单位复根、DFT 与 FFT；30.3 高效串行与并行实现。

30.1 多项式的表示（Representing polynomials）（PDF p.921–927）

系数表示（coefficient representation）

向量 $a=(a_0,\dots,a_{n-1})$（列向量）。
求值用 Horner 法则（Horner's rule）$\Theta(n)$：$A(x_0)=a_0+x_0(a_1+x_0(a_2+\cdots+x_0(a_{n-2}+x_0a_{n-1})\cdots))$。
加法 $\Theta(n)$；直接乘法 $\Theta(n^2)$。乘积系数向量 $c$ 称为 $a,b$ 的卷积（convolution），记 $c=a\otimes b$。

def horner(a, x0):          # a[0..n-1]
    y = 0
    for j in range(n-1, -1, -1):
        y = a[j] + x0 * y
    return y                # 时间 Θ(n)，空间 O(1)

点值表示（point-value representation）

次数界 $n$ 的多项式由 $n$ 个点值对 $\{(x_0,y_0),\dots,(x_{n-1},y_{n-1})\}$ 表示，$x_k$ 两两不同，$y_k=A(x_k)$。同一多项式有无数种点值表示。
由系数求点值 = 在 $n$ 点求值（evaluation），用 Horner 法需 $\Theta(n^2)$；巧选点可降到 $\Theta(n\lg n)$。逆运算是插值（interpolation）。

定理 30.1（插值多项式唯一性）：对任意 $n$ 个 $x_k$ 互异的点值对，存在唯一次数界为 $n$ 的多项式 $A$ 满足 $y_k=A(x_k)$。证明：$y=A(x_k)$ 等价于矩阵方程 $V(x_0,\dots,x_{n-1})\,a=y$，其中 $V$ 是范德蒙德矩阵（Vandermonde matrix），第 $k$ 行为 $(1,x_k,x_k^2,\dots,x_k^{n-1})$，行列式 $\prod_{0\le j<k\le n-1}(x_k-x_j)\ne0$，故可逆，$a=V^{-1}y$。

用 LU 分解解此方程组需 $O(n^3)$。更快的是拉格朗日公式（Lagrange's formula）： $A(x)=\sum_{k=0}^{n-1}y_k\frac{\prod_{j\ne k}(x-x_j)}{\prod_{j\ne k}(x_k-x_j)},\qquad(30.5)$ 可在 $\Theta(n^2)$ 时间内求出系数（习题 30.1-5）。
脚注提醒：插值在数值稳定性上是出名的棘手问题，输入的小扰动或舍入误差会导致结果大变。

点值表示下的运算

加法：同一组点上 $y_k+y'_k$，$\Theta(n)$。
乘法：逐点相乘 $y_ky'_k$，$\Theta(n)$。但 $C$ 的次数界为 $2n$，需要 $2n$ 个点才能唯一确定（习题 30.1-4），所以 $A,B$ 必须先用 $2n$ 个点的扩展点值表示（extended point-value representation）。
在新点求值：没有比先转回系数表示更简单的办法。

系数形式下的快速乘法（图 30.1）

关键在于能否快速在两种表示间转换。选单位复根（complex roots of unity）作求值点，求值即对系数向量做离散傅里叶变换（DFT），插值即逆 DFT，FFT 可在 $\Theta(n\lg n)$ 内完成两者。设 $n$ 为 2 的幂（不足补高位零系数）：

倍增次数界：$A,B$ 各补 $n$ 个高位零系数，成为次数界 $2n$ 的多项式，$\Theta(n)$；
求值：各做一次 $2n$ 阶 FFT，得到在 $2n$ 次单位复根处的值，$\Theta(n\lg n)$；
逐点相乘：得到 $C$ 在每个 $2n$ 次单位根处的值，$\Theta(n)$；
插值：对 $2n$ 个点值做逆 DFT（同样用 FFT）得 $C$ 的系数，$\Theta(n\lg n)$。

定理 30.2：两个次数界为 $n$ 的多项式可在 $\Theta(n\lg n)$ 时间内相乘，输入输出均为系数表示。

30.1 习题概览：30.1-1 用 (30.1)(30.2) 计算 $(7x^3-x^2+x-10)(8x^3-6x+3)$；30.1-2 用综合除法 $A(x)=q(x)(x-x_0)+r$ 在 $\Theta(n)$ 内求商与余数（$r=A(x_0)$）；30.1-3 由 $A$ 的点值表示推出逆序多项式 $A^{rev}(x)=\sum a_{n-1-j}x^j$ 的点值表示（$A^{rev}(x)=x^{n-1}A(1/x)$）；30.1-4 证明少于 $n$ 个点不能唯一确定次数界 $n$ 的多项式；30.1-5 用拉格朗日公式在 $\Theta(n^2)$ 内插值（先求 $\prod_j(x-x_j)$，再逐个除以 $(x-x_k)$）；30.1-6 说明点值表示下"直接相除 $y$ 值"做多项式除法的问题（整除与不整除两种情形）；30.1-7 笛卡尔和：$A,B$ 各含 $n$ 个 $[0,10n]$ 内的整数，求 $C=\{x+y\}$ 及每个和出现的次数，$O(n\lg n)$（把集合表示为 $\sum_{a\in A}x^a$ 的多项式再相乘——FFT 计数的经典技巧）。

30.2 DFT 与 FFT（The DFT and FFT）（PDF p.927–936）

单位复根（complex roots of unity）

$n$ 次单位复根：满足 $\omega^n=1$ 的复数 $\omega$，恰有 $n$ 个：$e^{2\pi ik/n}$，$k=0,\dots,n-1$。利用 $e^{iu}=\cos u+i\sin u$，它们在复平面单位圆上等距分布（图 30.2 画出 $n=8$ 的情形）。
主 $n$ 次单位根（principal $n$th root of unity）：$\omega_n=e^{2\pi i/n}$（30.6），其他单位根都是它的幂。脚注：信号处理文献多用 $\omega_n=e^{-2\pi i/n}$，数学实质相同。
$\omega_n^0,\dots,\omega_n^{n-1}$ 在乘法下构成群，与加法群 $(\mathbb Z_n,+)$ 同构：$\omega_n^j\omega_n^k=\omega_n^{(j+k)\bmod n}$，$\omega_n^{-1}=\omega_n^{n-1}$。

引理 30.3（消去引理，cancellation lemma）：对整数 $n\ge0,k\ge0,d>0$，$\omega_{dn}^{dk}=\omega_n^k$。证明：$(e^{2\pi i/dn})^{dk}=(e^{2\pi i/n})^k$。

推论 30.4：$n>0$ 为偶数时，$\omega_n^{n/2}=\omega_2=-1$。

引理 30.5（折半引理，halving lemma）：$n>0$ 为偶数，则 $n$ 个 $n$ 次单位复根的平方恰为 $n/2$ 个 $n/2$ 次单位复根（每个出现两次）。证明：$(\omega_n^k)^2=\omega_{n/2}^k$；且 $(\omega_n^{k+n/2})^2=\omega_n^{2k+n}=\omega_n^{2k}=(\omega_n^k)^2$（也可由 $\omega_n^{k+n/2}=-\omega_n^k$ 看出）。折半引理保证分治递归的子问题规模减半。

引理 30.6（求和引理，summation lemma）：$n\ge1$，非零整数 $k$ 不被 $n$ 整除，则 $\sum_{j=0}^{n-1}(\omega_n^k)^j=0$。证明：等比级数 $\frac{(\omega_n^k)^n-1}{\omega_n^k-1}=\frac{(\omega_n^n)^k-1}{\omega_n^k-1}=\frac{1-1}{\omega_n^k-1}=0$，分母因 $n\nmid k$ 不为 0。

DFT 定义

对系数向量 $a=(a_0,\dots,a_{n-1})$，在 $n$ 个 $n$ 次单位复根处求值：

y_k=A(\omega_n^k)=\sum_{j=0}^{n-1}a_j\omega_n^{kj},\quad k=0,\dots,n-1.\qquad(30.8)

$y=(y_0,\dots,y_{n-1})$ 称为 $a$ 的离散傅里叶变换（discrete Fourier transform, DFT），记 $y=\mathrm{DFT}_n(a)$。（脚注：在多项式乘法中这里的 $n$ 就是 30.1 节的 $2n$。）

FFT：分治

设 $n$ 为 2 的幂（非 2 的幂的方法超出本书范围）。按下标奇偶拆分：

A^{[0]}(x)=a_0+a_2x+\cdots+a_{n-2}x^{n/2-1},\qquad A^{[1]}(x)=a_1+a_3x+\cdots+a_{n-1}x^{n/2-1},

A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2).\qquad(30.9)

于是在 $\omega_n^0,\dots,\omega_n^{n-1}$ 求 $A$ 归结为：在 $(\omega_n^k)^2$ 处求两个次数界 $n/2$ 的多项式，再按 (30.9) 合并。由折半引理，这些平方值只是 $n/2$ 个 $n/2$ 次单位根，所以子问题与原问题同型、规模减半。

def RECURSIVE_FFT(a):             # len(a) = n 为 2 的幂；下标从 0 开始
    n = len(a)
    if n == 1:
        return a                  # 单元素的 DFT 是其本身：y0 = a0 * ω_1^0 = a0
    wn = exp(2j * pi / n)
    w = 1
    y0 = RECURSIVE_FFT(a[0::2])   # 偶下标系数
    y1 = RECURSIVE_FFT(a[1::2])   # 奇下标系数
    y = [0] * n
    for k in range(n // 2):
        t = w * y1[k]             # 旋转因子 ω_n^k 乘以 y1[k]
        y[k] = y0[k] + t
        y[k + n//2] = y0[k] - t
        w = w * wn
    return y

正确性：递归给出 $y^{[0]}_k=A^{[0]}(\omega_{n/2}^k)=A^{[0]}(\omega_n^{2k})$，$y^{[1]}_k=A^{[1]}(\omega_n^{2k})$。

$y_k=y^{[0]}_k+\omega_n^ky^{[1]}_k=A^{[0]}(\omega_n^{2k})+\omega_n^kA^{[1]}(\omega_n^{2k})=A(\omega_n^k)$；
$y_{k+n/2}=y^{[0]}_k-\omega_n^ky^{[1]}_k=A^{[0]}(\omega_n^{2k+n})+\omega_n^{k+n/2}A^{[1]}(\omega_n^{2k+n})=A(\omega_n^{k+n/2})$（用了 $\omega_n^{k+n/2}=-\omega_n^k$、$\omega_n^{2k+n}=\omega_n^{2k}$）。

因子 $\omega_n^k$ 以正负两种形式使用，称为旋转因子（twiddle factors）。维护运行变量 $\omega$ 而不是每次重新计算 $\omega_n^k$，节省时间。

复杂度：除递归外每层 $\Theta(n)$，$T(n)=2T(n/2)+\Theta(n)=\Theta(n\lg n)$；空间 $\Theta(n\lg n)$（若每层新建数组且不释放）或 $\Theta(n)$ 辅助空间（递归栈深 $\lg n$，逐层释放）。

在单位复根处插值：逆 DFT

DFT 写成矩阵乘积 $y=V_na$，$V_n$ 为范德蒙德矩阵，$(k,j)$ 元素为 $\omega_n^{kj}$（指数构成乘法表）。

定理 30.7：$V_n^{-1}$ 的 $(j,k)$ 元素为 $\omega_n^{-kj}/n$。证明：$[V_n^{-1}V_n]_{jj'}=\sum_{k=0}^{n-1}\omega_n^{k(j'-j)}/n$，$j'=j$ 时为 1，否则因 $-(n-1)\le j'-j\le n-1$ 不被 $n$ 整除，由求和引理为 0。

于是

a_j=\frac1n\sum_{k=0}^{n-1}y_k\omega_n^{-kj},\quad j=0,\dots,n-1.\qquad(30.11)

与 (30.8) 对比：把 FFT 中 $a$ 与 $y$ 互换、$\omega_n$ 换成 $\omega_n^{-1}$、结果每项除以 $n$，即得逆 DFT，同样 $\Theta(n\lg n)$（习题 30.2-4）。

定理 30.8（卷积定理，convolution theorem）：对长度为 $n$（2 的幂）的向量 $a,b$，

a\otimes b=\mathrm{DFT}_{2n}^{-1}\big(\mathrm{DFT}_{2n}(a)\cdot\mathrm{DFT}_{2n}(b)\big),

其中 $a,b$ 先补零到长度 $2n$，$\cdot$ 为逐分量乘积。

30.2 习题概览：30.2-1 证明推论 30.4；30.2-2 计算 $(0,1,2,3)$ 的 DFT（答案 $(6,-2-2i,-2,-2+2i)$）；30.2-3 用 FFT 方案重做 30.1-1；30.2-4 写逆 DFT 伪代码；30.2-5 $n$ 为 3 的幂时的 FFT 推广，递推 $T(n)=3T(n/3)+\Theta(n)$；30.2-6★ 在模 $m=2^{tn/2}+1$ 的整数环上以 $\omega=2^t$ 作主 $n$ 次单位根，证明 DFT 与逆 DFT 良定义（数论变换）；30.2-7 给定根 $z_0,\dots,z_{n-1}$ 构造以它们为根的多项式，$O(n\lg^2n)$（分治乘积树）；30.2-8★ chirp 变换 $y_k=\sum_ja_jz^{kj}$（任意复数 $z$），利用 $kj=\tfrac{k^2+j^2-(k-j)^2}{2}$ 改写为卷积 $y_k=z^{k^2/2}\sum_j(a_jz^{j^2/2})z^{-(k-j)^2/2}$，$O(n\lg n)$（Bluestein 算法，可处理任意长度 DFT）。

30.3 高效 FFT 实现（Efficient FFT implementations）（PDF p.936–941）

迭代版 FFT

动机：信号处理要求极致速度。迭代版同为 $\Theta(n\lg n)$，但常数可能更小（视实现而定，递归版有时对缓存更友好）。
公共子表达式（common subexpression）：$\omega_n^ky^{[1]}_k$ 在循环中算了两次，存入临时变量 $t$ 只算一次。"乘旋转因子、存入 $t$、与 $y^{[0]}_k$ 相加和相减"称为蝴蝶操作（butterfly operation，图 30.3）。
递归调用树（图 30.4，$n=8$）：根 $(a_0..a_7)$，下一层 $(a_0,a_2,a_4,a_6)$ 与 $(a_1,a_3,a_5,a_7)$，再下一层 $(a_0,a_4),(a_2,a_6),(a_1,a_5),(a_3,a_7)$，叶子顺序 $a_0,a_4,a_2,a_6,a_1,a_5,a_3,a_7$。若一开始就把输入排成叶子顺序，就可以自底向上执行：先两两做 1 次蝴蝶得 $n/2$ 个 2 元 DFT，再两两合并成 $n/4$ 个 4 元 DFT……直到一个 $n$ 元 DFT。
叶子顺序是位逆序置换（bit-reversal permutation）：$a_k$ 放到 $A[\mathrm{rev}(k)]$，$\mathrm{rev}(k)$ 为 $k$ 的 $\lg n$ 位二进制逆序。例：0,4,2,6,1,5,3,7 的二进制 000,100,010,110,001,101,011,111 逆序后恰为 0..7。原因：顶层按最低位分左右子树，每层剥掉最低位继续分。

def BIT_REVERSE_COPY(a, A):
    n = len(a)
    for k in range(n):
        A[rev(k)] = a[k]          # rev(k)：k 的 lg n 位二进制反转

def ITERATIVE_FFT(a):
    n = len(a); A = [0]*n
    BIT_REVERSE_COPY(a, A)
    for s in range(1, lg(n) + 1):         # 第 s 层：合并成 2^s 元 DFT
        m = 2 ** s
        wm = exp(2j * pi / m)
        for k in range(0, n, m):          # 每组
            w = 1
            for j in range(m // 2):       # 组内 m/2 个蝴蝶
                t = w * A[k + j + m//2]
                u = A[k + j]
                A[k + j] = u + t
                A[k + j + m//2] = u - t
                w = w * wm
    return A

复杂度：BIT-REVERSE-COPY 为 $O(n\lg n)$（每次反转 $O(\lg n)$ 位）；实践中 $n$ 预知时可查表做到 $\Theta(n)$，或用思考题 17-1 的摊还逆序二进制计数器。最内层循环执行次数 $L(n)=\sum_{s=1}^{\lg n}\frac{n}{2^s}\cdot2^{s-1}=\sum_{s=1}^{\lg n}\frac n2=\Theta(n\lg n)$。原地计算，额外空间 $O(1)$（除输出数组）。

并行 FFT 电路（图 30.5）

电路先做位逆序置换，再有 $\lg n$ 级（stage），每级 $n/2$ 个蝴蝶并行执行；第 $s$ 级含 $n/2^s$ 组、每组 $2^{s-1}$ 个蝴蝶，使用旋转因子 $\omega_m^0,\dots,\omega_m^{m/2-1}$（$m=2^s$）。
电路深度（任一输出到能到达它的输入之间的最大计算元件数）为 $\Theta(\lg n)$，总共 $\Theta(n\lg n)$ 个蝴蝶。只有蝴蝶上下两根线与之交互，穿过中间的线不受影响。

30.3 习题概览：30.3-1 跟踪 ITERATIVE-FFT 计算 $(0,2,3,-1,4,5,7,9)$ 的 DFT；30.3-2 把位逆序置换放到计算末尾（提示：考虑逆 DFT——即 decimation-in-frequency 版本）；30.3-3 每级计算旋转因子的次数，改写为第 $s$ 级只算 $2^{s-1}$ 次；30.3-4★ 电路中恰有一个加法器故障恒输出 0，如何通过输入输出定位它。

第 30 章思考题与章末注记（PDF p.941–946）

30-1 分治乘法：(a) $(ax+b)(cx+d)$ 只用 3 次乘法（$ac$、$bd$、$(a+b)(c+d)$，中间项 $=(a+b)(c+d)-ac-bd$）；(b) 两种 $\Theta(n^{\lg3})$ 多项式乘法（按高低半拆分、按奇偶拆分）——即 Karatsuba 思想；(c) 两个 $n$ 位整数相乘 $O(n^{\lg 3})$ 位操作。
30-2 Toeplitz 矩阵（$a_{ij}=a_{i-1,j-1}$，沿对角线为常数）：(a) 和是 Toeplitz，积不一定；(b) 用 $2n-1$ 个数表示，$O(n)$ 相加；(c) Toeplitz 矩阵乘向量 $O(n\lg n)$（化为卷积，用 FFT）；(d) 两个 Toeplitz 矩阵相乘的高效算法。
30-3 多维 FFT：$d$ 维 DFT $y_{k_1..k_d}=\sum a_{j_1..j_d}\omega_{n_1}^{j_1k_1}\cdots\omega_{n_d}^{j_dk_d}$。(a) 可依次沿每一维做一维 DFT；(b) 维的顺序无关；(c) 总时间 $O(n\lg n)$，与 $d$ 无关（$n=n_1\cdots n_d$）。
30-4 在一点求多项式所有阶导数：(a) 若已知 $A(x)=\sum b_j(x-x_0)^j$，则 $A^{(t)}(x_0)=t!\,b_t$，$O(n)$；(b) 由 $A(x_0+\omega_n^k)$ 用逆 DFT 在 $O(n\lg n)$ 求 $b_j$；(c) 证明 $A(x_0+\omega_n^k)=\sum_r\frac{\omega_n^{kr}}{r!}\sum_jf(j)g(r-j)$，其中 $f(j)=a_jj!$，$g(l)=x_0^{-l}/(-l)!$（$-(n-1)\le l\le0$），否则 0；(d) 因此可在 $O(n\lg n)$ 内求出全部非平凡导数。
30-5 多点求值：假设多项式取模可 $O(n\lg n)$ 完成（如 $(3x^3+x^2-3x+1)\bmod(x^2+x+2)=-7x+5$）。定义 $P_{ij}(x)=\prod_{k=i}^j(x-x_k)$，$Q_{ij}=A\bmod P_{ij}$。(a) $A(x)\bmod(x-z)=A(z)$；(b) $Q_{kk}=A(x_k)$，$Q_{0,n-1}=A$；(c) $Q_{ik}=Q_{ij}\bmod P_{ik}$；(d) 得到 $O(n\lg^2n)$ 的多点求值算法（余数树）。
30-6 模算术 FFT（数论变换，NTT）：复数 FFT 有舍入误差；整数系数多项式乘法可用模运算得精确结果。(a) 找最小 $k$ 使 $p=kn+1$ 为素数，启发式论证 $k\approx\ln n$，$p$ 的长度约为 $n$ 长度的常数倍（$O(\lg n)$ 位）；(b) 令 $g$ 为 $\mathbb Z_p^*$ 生成元，$w=g^k\bmod p$ 作主 $n$ 次单位根，DFT 与逆 DFT 在模 $p$ 下互逆；(c) $O(n\lg n)$ 实现；(d) 计算 $(0,5,3,7,7,2,1,6)$ 模 17 的 DFT（$g=3$）。

章末注记：参考 Van Loan、Press 等《Numerical Recipes》、Oppenheim–Schafer（含非 2 的幂长度）、Oppenheim–Willsky；图像处理中的多维 FFT 见 Gonzalez–Woods、Pratt。FFT 通常归功于 Cooley 与 Tukey（1960 年代），但此前多次被发现，Heideman 等追溯到 1805 年的高斯。Frigo 与 Johnson 的 FFTW（"fastest Fourier transform in the West"）先运行"规划器"（planner）试跑选择最优分解，适配缓存，小规模子问题用优化的直线代码，且对任意 $n$（包括大素数）都是 $\Theta(n\lg n)$。非等距数据的近似 FFT 见 Ware 的综述。

第 30 章本章要点

多项式有系数表示与点值表示：前者求值/加法快，乘法 $\Theta(n^2)$；后者加法/乘法 $\Theta(n)$。两者之间的转换（求值/插值）选单位复根时可用 FFT 在 $\Theta(n\lg n)$ 完成。
单位复根三条性质（消去、折半、求和引理）是 FFT 和逆 FFT 成立的根本。
FFT 按奇偶下标分治：$A(x)=A^{[0]}(x^2)+xA^{[1]}(x^2)$，$T(n)=2T(n/2)+\Theta(n)$。
逆 DFT 矩阵 $V_n^{-1}=\frac1n\overline{V_n}$，只需换 $\omega_n\to\omega_n^{-1}$ 并除以 $n$。
卷积定理：卷积 = 逐点乘积的逆变换；注意必须补零到 $2n$ 以避免循环卷积的混叠（wrap-around）。
工程实现：迭代版 + 位逆序 + 原地蝴蝶；并行电路深度 $\Theta(\lg n)$。整数精确计算可用模素数的数论变换。

第 30 章与量化交易的关联

快速卷积/相关：滚动加权（如指数或任意核的移动平均）、长序列的自相关函数和互相关函数（例如多资产领先滞后关系、成交量季节性）都可用 FFT 在 $O(n\lg n)$ 内计算；注意需要补零避免循环卷积的边界混叠，这是实操中最常见的错误。
期权定价：Carr–Madan（1999）方法用 FFT 一次性对一整条执行价网格计算欧式期权价格，前提是已知对数价格的特征函数（Heston、Variance Gamma 等模型）；30-2/习题 30.2-8 中的 chirp/Bluestein 技巧对应分数 FFT（fractional FFT），用于独立选择执行价与积分网格；COS 方法也属同一类思想。
谱分析与周期检测：功率谱、周期图（periodogram）用于检测日内成交量的周期模式和季节性；频域滤波可做去噪，但金融序列非平稳、信噪比低，频域滤波在回测中易引入前视偏差（使用了未来数据的全样本变换），要用因果滤波。
概率分布卷积：独立收益的和的分布 = 分布的卷积，可用 FFT 计算组合损失分布（信用组合违约损失的离散分布、聚合风险模型），30.1-7 的"集合笛卡尔和计数"是同一技巧。
系统实现：实际使用 numpy/scipy.fft（底层为 pocketfft）或 FFTW；理解位逆序和蝴蝶有助于在 GPU/FPGA 上做低延迟信号处理。

第 30 章推荐习题

30.2-2：手算小规模 DFT，熟悉单位根。
30.1-7：用多项式乘法做计数卷积，是 FFT 最常用的建模套路。
30.2-8（chirp 变换/Bluestein）：任意长度 DFT 与分数 FFT 的基础，直接用于 FFT 期权定价。
思考题 30-1：Karatsuba 分治思想。
思考题 30-2(c)：Toeplitz 矩阵-向量乘法化为卷积——平稳时间序列协方差矩阵就是 Toeplitz 矩阵。
思考题 30-6：数论变换，理解精确卷积。

第 31 章数论算法（Number-Theoretic Algorithms）（PDF p.947–1005）

第 31 章引言（PDF p.947–948）

数论曾被视为优美但无用的纯数学；如今因基于大素数的密码体制而广泛应用：这些体制可行是因为大素数容易找到，安全是因为不知道如何高效分解大素数之积（或求离散对数）。
章节安排：31.1 整除、同余、唯一分解；31.2 欧几里得算法求 gcd；31.3 模运算（群论）；31.4 用欧几里得算法解 $ax\equiv b\pmod n$；31.5 中国剩余定理；31.6 模幂与反复平方法（素性测试和密码学的核心）；31.7 RSA 公钥密码；31.8 随机化素性测试（用于生成 RSA 密钥）；31.9 小整数分解的启发式（Pollard rho）。分解恰是人们希望它困难的问题。

输入规模与算术代价：

这里"大输入"指"大整数"而非"多个整数"。输入规模按二进制位数度量；输入为 $a_1,\dots,a_k$ 的算法若运行时间是 $\lg a_1,\dots,\lg a_k$ 的多项式，称为多项式时间算法。
大整数的算术不再是单位时间，因此改用位操作（bit operations）计数：普通方法乘两个 $\beta$ 位整数需 $\Theta(\beta^2)$ 位操作；$\beta$ 位数除以较短数或取余也是 $\Theta(\beta^2)$（习题 31.1-12）。更快的方法：分治乘法 $\Theta(\beta^{\lg3})$，已知最快 $\Theta(\beta\lg\beta\lg\lg\beta)$（Schönhage–Strassen）。实践中常用 $\Theta(\beta^2)$ 作分析基准。本章同时按算术运算次数和位操作数分析。

31.1 初等数论概念（Elementary number-theoretic notions）（PDF p.948–954）

整数集 $\mathbb Z$，自然数集 $\mathbb N=\{0,1,2,\dots\}$。

整除与约数：

$d\mid a$（$d$ 整除 $a$）指存在整数 $k$ 使 $a=kd$。每个整数都整除 0。$a>0$ 且 $d\mid a$ 则 $|d|\le|a|$。$a$ 是 $d$ 的倍数（multiple）；不整除记 $d\nmid a$。
$d\mid a$ 且 $d\ge0$ 时称 $d$ 为 $a$ 的约数（divisor）。非零整数 $a$ 的约数在 $1$ 到 $|a|$ 之间。例：24 的约数为 1,2,3,4,6,8,12,24。
1 和 $a$ 是平凡约数（trivial divisors），其余称为因子（factors），如 20 的因子为 2,4,5,10。

素数与合数：$a>1$ 且只有平凡约数，称为素数（prime）。前 20 个素数：2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71。素数有无穷多（习题 31.1-2）。$a>1$ 且非素数为合数（composite），如 39（$3\mid39$）。1 称为单位（unit），既非素数也非合数；0 和负数也都不是。

定理 31.1（带余除法，division theorem）：对任意整数 $a$ 和正整数 $n$，存在唯一整数 $q,r$ 使 $0\le r<n$ 且 $a=qn+r$。$q=\lfloor a/n\rfloor$ 为商（quotient），$r=a\bmod n$ 为余数（remainder/residue）。$n\mid a\iff a\bmod n=0$。

同余类：模 $n$ 的等价类 $[a]_n=\{a+kn:k\in\mathbb Z\}$，例：$[3]_7=\{\dots,-11,-4,3,10,17,\dots\}$，也可记 $[-4]_7$、$[10]_7$。$a\in[b]_n$ 等价于 $a\equiv b\pmod n$。$\mathbb Z_n=\{[a]_n:0\le a\le n-1\}$（31.1），常简写为 $\{0,1,\dots,n-1\}$（31.2），每类用最小非负元代表；说 $-1\in\mathbb Z_n$ 实指 $[n-1]_n$。

公约数与最大公约数：

同时整除 $a,b$ 的 $d$ 为公约数（common divisor）。例：30 的约数 1,2,3,5,6,10,15,30，24 与 30 的公约数为 1,2,3,6。
性质：$d\mid a,\ d\mid b\Rightarrow d\mid(a\pm b)$（31.3），更一般 $d\mid(ax+by)$ 对任意整数 $x,y$（31.4）。若 $a\mid b$ 则 $|a|\le|b|$ 或 $b=0$，故 $a\mid b$ 且 $b\mid a\Rightarrow a=\pm b$（31.5）。
最大公约数 $\gcd(a,b)$（$a,b$ 不全为 0）：如 $\gcd(24,30)=6$，$\gcd(5,7)=1$，$\gcd(0,9)=9$。定义 $\gcd(0,0)=0$ 以使性质普遍成立。
基本性质：$\gcd(a,b)=\gcd(b,a)=\gcd(-a,b)=\gcd(|a|,|b|)$；$\gcd(a,0)=|a|$；$\gcd(a,ka)=|a|$（31.6–31.10）。

定理 31.2：$a,b$ 不全为 0，则 $\gcd(a,b)$ 是线性组合集合 $\{ax+by:x,y\in\mathbb Z\}$ 中最小的正元素。证明：设 $s=ax+by$ 为最小正线性组合，$q=\lfloor a/s\rfloor$，$a\bmod s=a-qs=a(1-qx)+b(-qy)$ 也是线性组合，又 $0\le a\bmod s<s$，故 $a\bmod s=0$，即 $s\mid a$；同理 $s\mid b$，于是 $\gcd(a,b)\ge s$。另一方面 $\gcd(a,b)$ 整除 $a,b$ 故整除 $s$，$s>0$ 推出 $\gcd(a,b)\le s$。

推论 31.3：$d\mid a,\ d\mid b\Rightarrow d\mid\gcd(a,b)$。 推论 31.4：$\gcd(an,bn)=n\gcd(a,b)$（$n\ge0$）。 推论 31.5：$n\mid ab$ 且 $\gcd(a,n)=1\Rightarrow n\mid b$（习题 31.1-5）。

互素：$\gcd(a,b)=1$ 称 $a,b$ 互素（relatively prime），如 8 与 15。 定理 31.6：$\gcd(a,p)=\gcd(b,p)=1\Rightarrow\gcd(ab,p)=1$。证明：$ax+py=1$，$bx'+py'=1$，相乘得 $ab(xx')+p(ybx'+y'ax+pyy')=1$，再用定理 31.2。 $n_1,\dots,n_k$ 两两互素（pairwise relatively prime）：任意 $i\ne j$ 有 $\gcd(n_i,n_j)=1$。

唯一分解： 定理 31.7：素数 $p\mid ab\Rightarrow p\mid a$ 或 $p\mid b$。反证：若都不整除，则 $\gcd(a,p)=\gcd(b,p)=1$，由定理 31.6 得 $\gcd(ab,p)=1$，与 $p\mid ab$ 矛盾。 定理 31.8（唯一分解定理）：合数 $a$ 唯一地写成 $a=p_1^{e_1}p_2^{e_2}\cdots p_r^{e_r}$，$p_1<\cdots<p_r$ 为素数，$e_i$ 为正整数。例：$6000=2^4\cdot3\cdot5^3$。

31.1 习题概览：31.1-1 $a>b>0$、$c=a+b$ 则 $c\bmod a=b$；31.1-2 素数无穷多（$p_1\cdots p_k+1$ 不被任何 $p_i$ 整除）；31.1-3 整除传递性；31.1-4 素数 $p$ 与 $0<k<p$ 互素；31.1-5 证推论 31.5；31.1-6 $p\mid\binom pk$（$0<k<p$），推出 $(a+b)^p\equiv a^p+b^p\pmod p$（"新生之梦"）；31.1-7 $a\mid b$ 时 $(x\bmod b)\bmod a=x\bmod a$；31.1-8 判断 $\beta$ 位整数是否为非平凡幂（$n=a^k$，$k>1$），多项式时间（对每个 $k\le\beta$ 二分查找 $a$）；31.1-9 证 31.6–31.10；31.1-10 gcd 满足结合律；31.1-11★ 证唯一分解；31.1-12 $\beta$ 位数除以短整数与取余的 $\Theta(\beta^2)$ 算法；31.1-13 二进制转十进制，若乘除法耗时 $M(\beta)$，可分治在 $\Theta(M(\beta)\lg\beta)$ 完成。

31.2 最大公约数（Greatest common divisor）（PDF p.954–960）

只考虑非负整数（因 $\gcd(a,b)=\gcd(|a|,|b|)$）。
由素因子分解 $a=\prod p_i^{e_i}$、$b=\prod p_i^{f_i}$（不出现的素数指数记 0）可得 $\gcd(a,b)=\prod p_i^{\min(e_i,f_i)}$（31.13），但分解没有多项式时间算法，此路不通。

定理 31.9（GCD 递归定理）：对非负整数 $a$ 和正整数 $b$，$\gcd(a,b)=\gcd(b,a\bmod b)$。证明：两者互相整除。令 $d=\gcd(a,b)$，$a\bmod b=a-\lfloor a/b\rfloor b$ 是 $a,b$ 的线性组合，故 $d\mid(a\bmod b)$，结合 $d\mid b$ 得 $d\mid\gcd(b,a\bmod b)$；反向同理，因 $a=qb+(a\bmod b)$。再由 (31.5) 得相等。

欧几里得算法（Euclid's algorithm）

《几何原本》（约公元前 300 年）记载，可能更早。

def EUCLID(a, b):          # a, b 为非负整数
    if b == 0:
        return a
    return EUCLID(b, a % b)

例：$\mathrm{EUCLID}(30,21)=\mathrm{EUCLID}(21,9)=\mathrm{EUCLID}(9,3)=\mathrm{EUCLID}(3,0)=3$，递归 3 次。正确性：定理 31.9 + $\gcd(a,0)=a$；第二个参数严格递减且非负，必终止。

运行时间分析：与斐波那契数的联系

不妨设 $a>b\ge0$（若 $b>a$，第一次递归只是交换参数；若 $a=b>0$ 一次即结束）。运行时间正比于递归次数。

引理 31.10：若 $a>b\ge1$ 且 EUCLID$(a,b)$ 执行 $k\ge1$ 次递归调用，则 $a\ge F_{k+2}$、$b\ge F_{k+1}$。证明（对 $k$ 归纳）：$k=1$ 时 $b\ge1=F_2$，$a>b$ 得 $a\ge2=F_3$。每次递归第一个参数严格大于第二个。归纳步：EUCLID$(b,a\bmod b)$ 做 $k-1$ 次递归，由归纳假设 $b\ge F_{k+1}$，$a\bmod b\ge F_k$。又因 $\lfloor a/b\rfloor\ge1$，$b+(a\bmod b)=b+(a-b\lfloor a/b\rfloor)\le a$，故 $a\ge F_{k+1}+F_k=F_{k+2}$。

定理 31.11（Lamé 定理）：对整数 $k\ge1$，若 $a>b\ge1$ 且 $b<F_{k+1}$，则 EUCLID$(a,b)$ 的递归调用少于 $k$ 次。

上界是紧的：EUCLID$(F_{k+1},F_k)$ 恰好递归 $k-1$ 次（$k\ge2$）。归纳：$F_{k+1}\bmod F_k=F_{k-1}$（习题 31.1-1），所以 $\gcd(F_{k+1},F_k)=\gcd(F_k,F_{k-1})$，多一次递归。相邻斐波那契数是欧几里得算法的最坏输入。
因 $F_k\approx\phi^k/\sqrt5$（$\phi=(1+\sqrt5)/2$ 黄金分割比），递归次数为 $O(\lg b)$（更紧的界见习题 31.2-5）。
对两个 $\beta$ 位数：$O(\beta)$ 次算术运算、$O(\beta^3)$ 位操作（每次乘除 $O(\beta^2)$）；思考题 31-2 可证明实际是 $O(\beta^2)$ 位操作。

扩展欧几里得算法（extended form of Euclid's algorithm）

同时求出系数 $x,y$（可为 0 或负数）使

d=\gcd(a,b)=ax+by.\qquad(31.16)

这些系数后面用于求模逆元。

def EXTENDED_EUCLID(a, b):     # 返回 (d, x, y)，d = gcd(a,b) = a*x + b*y
    if b == 0:
        return (a, 1, 0)
    d1, x1, y1 = EXTENDED_EUCLID(b, a % b)
    return (d1, y1, x1 - (a // b) * y1)

推导：递归返回 $d'=bx'+(a\bmod b)y'$，代入 $a\bmod b=a-b\lfloor a/b\rfloor$ 得 $d=ay'+b(x'-\lfloor a/b\rfloor y')$，所以 $x=y'$，$y=x'-\lfloor a/b\rfloor y'$。递归次数与 EUCLID 相同，$O(\lg b)$；递归栈空间 $O(\lg b)$（可改写为迭代版 $O(1)$ 额外空间）。

图 31.1：EXTENDED-EUCLID(99,78) 的递归过程

$a$	$b$	$\lfloor a/b\rfloor$	$d$	$x$	$y$
99	78	1	3	−11	14
78	21	3	3	3	−11
21	15	1	3	−2	3
15	6	2	3	1	−2
6	3	2	3	0	1
3	0	—	3	1	0

结果：$\gcd(99,78)=3=99\cdot(-11)+78\cdot14$。每层返回的 $(d,x,y)$ 成为上一层的 $(d',x',y')$。

31.2 习题概览：31.2-1 由分解式推出 (31.13)；31.2-2 计算 EXTENDED-EUCLID(899,493)（结果 $d=29$）；31.2-3 $\gcd(a,n)=\gcd(a+kn,n)$；31.2-4 改写为只用常数存储的迭代版；31.2-5 递归次数至多 $1+\log_\phi b$，并改进为 $1+\log_\phi(b/\gcd(a,b))$；31.2-6 EXTENDED-EUCLID$(F_{k+1},F_k)$ 返回什么（系数为 $\pm$ 斐波那契数）；31.2-7 多参数 gcd 与顺序无关、求 $\gcd(a_0..a_n)=\sum a_ix_i$ 的系数，除法次数 $O(n+\lg\max a_i)$；31.2-8 用 gcd 计算多个数的最小公倍数 lcm（$\mathrm{lcm}(a,b)=ab/\gcd(a,b)$ 逐个累积）；31.2-9 $n_1..n_4$ 两两互素 ⟺ $\gcd(n_1n_2,n_3n_4)=\gcd(n_1n_3,n_2n_4)=1$，推广为 $\lceil\lg k\rceil$ 对乘积互素即可判定。

31.3 模运算（Modular arithmetic）（PDF p.960–967）

非正式地，模 $n$ 运算就是普通整数运算，但每个结果 $x$ 换成 $x\bmod n$。对加、减、乘足够；正式模型用群论描述。

有限群

群（group）$(S,\oplus)$：集合 $S$ 及其上的二元运算，满足

封闭性（closure）：$a\oplus b\in S$；
单位元（identity）：存在 $e$ 使 $e\oplus a=a\oplus e=a$；
结合律（associativity）：$(a\oplus b)\oplus c=a\oplus(b\oplus c)$；
逆元（inverses）：每个 $a$ 有唯一 $b$ 使 $a\oplus b=b\oplus a=e$。

例：$(\mathbb Z,+)$，单位元 0，$a$ 的逆为 $-a$。满足交换律的称为阿贝尔群（abelian group）；$|S|<\infty$ 称为有限群（finite group）。

模加法群与模乘法群

若 $a\equiv a'$、$b\equiv b'\pmod n$，则 $a+b\equiv a'+b'$，$ab\equiv a'b'\pmod n$。于是定义 $[a]_n+_n[b]_n=[a+b]_n$，$[a]_n\cdot_n[b]_n=[ab]_n$（31.18）；减法类似，除法较复杂。这为"用最小非负代表元计算、结果取 mod"的惯例提供依据。
模 $n$ 加法群 $(\mathbb Z_n,+_n)$，大小 $n$。图 31.2(a) 给出 $(\mathbb Z_6,+_6)$ 运算表。

定理 31.12：$(\mathbb Z_n,+_n)$ 是有限阿贝尔群。封闭性由 (31.18)；结合律、交换律继承自 $+$；单位元 0；$a$ 的逆元为 $-a$（即 $[n-a]_n$）。

模 $n$ 乘法群 $(\mathbb Z_n^*,\cdot_n)$，$\mathbb Z_n^*=\{[a]_n\in\mathbb Z_n:\gcd(a,n)=1\}$。良定义性：由习题 31.2-3，$\gcd(a,n)=1\Rightarrow\gcd(a+kn,n)=1$。例：$\mathbb Z_{15}^*=\{1,2,4,7,8,11,13,14\}$，图 31.2(b) 给出乘法表，如 $8\cdot11\equiv13\pmod{15}$，单位元为 1。

定理 31.13：$(\mathbb Z_n^*,\cdot_n)$ 是有限阿贝尔群。封闭性由定理 31.6；单位元 $[1]_n$；逆元存在性：对 $a\in\mathbb Z_n^*$，EXTENDED-EUCLID$(a,n)$ 返回 $d=1$ 且 $ax+ny=1$（31.19），即 $ax\equiv1\pmod n$，$[x]_n$ 即逆元；又由定理 31.2，$\gcd(x,n)=1$，故 $[x]_n\in\mathbb Z_n^*$。唯一性留到推论 31.26。例：$a=5,n=11$，EXTENDED-EUCLID 返回 $(1,-2,1)$，$1=5\cdot(-2)+11\cdot1$，故 $5^{-1}\equiv-2\equiv9\pmod{11}$。

记号约定：用代表元表示等价类，$+_n,\cdot_n$ 记为普通 $+$、$\cdot$；$ax\equiv b\pmod n$ 等价于 $[a]_n\cdot_n[x]_n=[b]_n$；群简称 $\mathbb Z_n$、$\mathbb Z_n^*$。逆元记 $a^{-1}\bmod n$；$\mathbb Z_n^*$ 中除法定义为 $a/b\equiv ab^{-1}\pmod n$。例：$7^{-1}\equiv13\pmod{15}$（$7\cdot13=91\equiv1$），所以 $4/7\equiv4\cdot13\equiv7\pmod{15}$。

欧拉 φ 函数（Euler's phi function）

$|\mathbb Z_n^*|$ 记为 $\phi(n)$：

\phi(n)=n\prod_{p\text{ 素数},\ p\mid n}\Big(1-\frac1p\Big).\qquad(31.20)

直观：从 $0..n-1$ 中对每个整除 $n$ 的素数 $p$ 划掉其倍数。例：$\phi(45)=45\cdot\tfrac23\cdot\tfrac45=24$。$p$ 为素数时 $\mathbb Z_p^*=\{1,\dots,p-1\}$，$\phi(p)=p-1$（31.21）。$n$ 为合数时 $\phi(n)<n-1$，但有下界：$n\ge3$ 时 $\phi(n)>\dfrac{n}{e^{\gamma}\ln\ln n+\frac{3}{\ln\ln n}}$（31.22），$\gamma=0.5772156649\ldots$ 为欧拉常数；更简单的 $n>5$ 时 $\phi(n)>\dfrac{n}{6\ln\ln n}$（31.23）；且 $\liminf_{n\to\infty}\dfrac{\phi(n)}{n/\ln\ln n}=e^{-\gamma}$（31.24），故 (31.22) 本质上最优。

子群（Subgroups）

$S'\subseteq S$ 且 $(S',\oplus)$ 也是群，称为子群。例：偶数是 $(\mathbb Z,+)$ 的子群。
定理 31.14：有限群的非空、对运算封闭的子集是子群（习题 31.3-3）。例：$\{0,2,4,6\}$ 是 $\mathbb Z_8$ 的子群。
定理 31.15（拉格朗日定理，Lagrange's theorem）：有限群 $S$ 的子群 $S'$ 满足 $|S'|$ 整除 $|S|$。（证明略）
$S'\ne S$ 时称为真子群（proper subgroup）。推论 31.16：真子群大小 $\le|S|/2$。（用于 31.8 节 Miller–Rabin 分析。）

由一个元素生成的子群

定义 $a^{(k)}=\underbrace{a\oplus a\oplus\cdots\oplus a}_{k}$。例：$\mathbb Z_6$ 中 $a=2$，序列为 $2,4,0,2,4,0,\dots$。在 $\mathbb Z_n$ 中 $a^{(k)}=ka\bmod n$，在 $\mathbb Z_n^*$ 中 $a^{(k)}=a^k\bmod n$。
由 $a$ 生成的子群 $\langle a\rangle=\{a^{(k)}:k\ge1\}$，$a$ 称为 $\langle a\rangle$ 的生成元（generator）。因 $a^{(i)}\oplus a^{(j)}=a^{(i+j)}$，$\langle a\rangle$ 封闭，故为子群。例：$\mathbb Z_6$ 中 $\langle0\rangle=\{0\}$，$\langle1\rangle=\{0,..,5\}$，$\langle2\rangle=\{0,2,4\}$；$\mathbb Z_7^*$ 中 $\langle1\rangle=\{1\}$，$\langle2\rangle=\{1,2,4\}$，$\langle3\rangle=\{1,..,6\}$。
$a$ 的阶（order）$\mathrm{ord}(a)$：使 $a^{(t)}=e$ 的最小正整数 $t$。

定理 31.17：$\mathrm{ord}(a)=|\langle a\rangle|$。证明：设 $t=\mathrm{ord}(a)$，$a^{(t+k)}=a^{(k)}$，故 $a^{(t)}$ 之后不出现新元素，$|\langle a\rangle|\le t$；若 $a^{(i)}=a^{(j)}$（$1\le i<j\le t$），则 $a^{(i+(t-j))}=a^{(j+(t-j))}=e$，而 $i+t-j<t$，矛盾，故前 $t$ 个互不相同。

推论 31.18：序列 $a^{(1)},a^{(2)},\dots$ 以 $t=\mathrm{ord}(a)$ 为周期：$a^{(i)}=a^{(j)}\iff i\equiv j\pmod t$。约定 $a^{(0)}=e$，$a^{(i)}=a^{(i\bmod t)}$。

推论 31.19：有限群 $S$ 中任意 $a$ 满足 $a^{(|S|)}=e$。证明：拉格朗日定理给出 $\mathrm{ord}(a)\mid|S|$。

31.3 习题概览：31.3-1 画 $(\mathbb Z_4,+_4)$ 与 $(\mathbb Z_5^*,\cdot_5)$ 运算表并给出同构映射；31.3-2 列出 $\mathbb Z_9$ 与 $\mathbb Z_{13}^*$ 的全部子群；31.3-3 证定理 31.14；31.3-4 $\phi(p^e)=p^{e-1}(p-1)$；31.3-5 对 $a\in\mathbb Z_n^*$，$f_a(x)=ax\bmod n$ 是 $\mathbb Z_n^*$ 上的置换。

31.4 求解模线性方程（Solving modular linear equations）（PDF p.967–971）

问题：给定 $a>0,n>0,b$，求所有满足

ax\equiv b\pmod n\qquad(31.25)

的 $x$（模 $n$），可能 0 个、1 个或多个解。用于 31.7 节 RSA 密钥生成。

设 $\langle a\rangle$ 为 $\mathbb Z_n$ 中 $a$ 生成的子群 $=\{ax\bmod n:x>0\}$，方程有解 ⟺ $[b]\in\langle a\rangle$。

定理 31.20：$d=\gcd(a,n)$，则在 $\mathbb Z_n$ 中 $\langle a\rangle=\langle d\rangle=\{0,d,2d,\dots,(n/d-1)d\}$（31.26），故 $|\langle a\rangle|=n/d$。证明：EXTENDED-EUCLID 给出 $ax'+ny'=d$，即 $ax'\equiv d$，所以 $d\in\langle a\rangle$，进而 $\langle d\rangle\subseteq\langle a\rangle$；反之 $m\in\langle a\rangle$ 时 $m=ax+ny$，$d\mid a$、$d\mid n$ 推出 $d\mid m$，故 $\langle a\rangle\subseteq\langle d\rangle$。$0..n-1$ 中 $d$ 的倍数恰有 $n/d$ 个。

推论 31.21：$ax\equiv b\pmod n$ 有解 ⟺ $d\mid b$。 推论 31.22：方程要么无解，要么模 $n$ 恰有 $d$ 个不同解。证明：序列 $ai\bmod n$ 以 $n/d$ 为周期，$i=0..n-1$ 中长度为 $n/d$ 的块重复 $d$ 次，$b$ 恰出现 $d$ 次。

定理 31.23：设 $d=ax'+ny'$，若 $d\mid b$，则 $x_0=x'(b/d)\bmod n$ 是一个解。证明：$ax_0\equiv ax'(b/d)\equiv d(b/d)=b$。

定理 31.24：若有解且 $x_0$ 为任一解，则恰有 $d$ 个解 $x_i=x_0+i(n/d)$，$i=0,\dots,d-1$。证明：它们模 $n$ 互不相同；$a\cdot i(n/d)=i(a/d)n$ 是 $n$ 的倍数，故都是解；再由推论 31.22 知没有别的解。

def MODULAR_LINEAR_EQUATION_SOLVER(a, b, n):
    d, x1, y1 = EXTENDED_EUCLID(a, n)
    if b % d == 0:
        x0 = (x1 * (b // d)) % n
        for i in range(d):
            print((x0 + i * (n // d)) % n)
    else:
        print("no solutions")

例：$14x\equiv30\pmod{100}$。EXTENDED-EUCLID(14,100) 返回 $(2,-7,1)$；$2\mid30$，$x_0=(-7)(15)\bmod100=95$；两个解 95 和 45。
复杂度：$O(\lg n+\gcd(a,n))$ 次算术运算（扩展欧几里得 $O(\lg n)$，循环 $d$ 次）。

推论 31.25：$n>1$，$\gcd(a,n)=1$ 时 $ax\equiv b\pmod n$ 模 $n$ 有唯一解。 推论 31.26：$n>1$，$\gcd(a,n)=1$ 时 $ax\equiv1\pmod n$ 有唯一解，否则无解。于是记号 $a^{-1}\bmod n$ 良定义，并可直接取 EXTENDED-EUCLID 返回的 $x$（因 $1=ax+ny$）。求模逆元 = 扩展欧几里得，$O(\lg n)$ 次算术运算。

31.4 习题概览：31.4-1 解 $35x\equiv10\pmod{50}$（$d=5$，化为 $7x\equiv2\pmod{10}$ 得 $x\equiv6$，5 个解为 $6,16,26,36,46$）；31.4-2 $\gcd(a,n)=1$ 时 $ax\equiv ay\Rightarrow x\equiv y$，并给出 $\gcd>1$ 的反例；31.4-3 若第 3 行改为 $x_0=x'(b/d)\bmod(n/d)$ 是否可行；31.4-4★ 模素数 $p$ 的 $t$ 次多项式至多有 $t$ 个不同根（先证 $a$ 为根则 $f(x)\equiv(x-a)g(x)$）。

31.5 中国剩余定理（The Chinese remainder theorem）（PDF p.971–975）

历史：约公元 100 年，中国数学家孙子（Sun-Tsŭ）求被 3、5、7 除分别余 2、3、2 的整数，一个解为 23，全部解为 $23+105k$。
定理建立"两两互素模数下的方程组"与"模它们乘积的一个方程"之间的对应。两大用途：(1) 结构定理：$\mathbb Z_n$ 与笛卡尔积 $\mathbb Z_{n_1}\times\cdots\times\mathbb Z_{n_k}$（各分量独立做模 $n_i$ 加乘）结构相同；(2) 算法设计：在各个小模数下计算（位操作更少）比直接模 $n$ 更高效。

定理 31.27（中国剩余定理）：$n=n_1n_2\cdots n_k$，$n_i$ 两两互素。对应

a\leftrightarrow(a_1,\dots,a_k),\quad a\in\mathbb Z_n,\ a_i=a\bmod n_i\qquad(31.27)

是 $\mathbb Z_n$ 与 $\mathbb Z_{n_1}\times\cdots\times\mathbb Z_{n_k}$ 之间的双射；且加、减、乘可在各分量上独立进行： $(a\pm b)\bmod n\leftrightarrow((a_1\pm b_1)\bmod n_1,\dots)$，$(ab)\bmod n\leftrightarrow(a_1b_1\bmod n_1,\dots)$（31.28–31.30）。

证明（构造逆映射）：

正向只需 $k$ 次取模。
反向：令 $m_i=n/n_i$（除 $n_i$ 外所有 $n_j$ 之积），$c_i=m_i\,(m_i^{-1}\bmod n_i)$（31.31）。由定理 31.6，$m_i$ 与 $n_i$ 互素，逆元存在。则 $a\equiv a_1c_1+a_2c_2+\cdots+a_kc_k\pmod n.\qquad(31.32)$
验证：$j\ne i$ 时 $m_j\equiv0\pmod{n_i}$，故 $c_j\equiv0$；且 $c_i\equiv1\pmod{n_i}$。于是 $c_i\leftrightarrow(0,\dots,0,1,0,\dots,0)$，$c_i$ 构成表示的"基"，$a\equiv a_ic_i\equiv a_i\pmod{n_i}$。双向可转换故为双射；运算对应由习题 31.1-7（$x\bmod n_i=(x\bmod n)\bmod n_i$）得出。

def CRT(a_list, n_list):           # n_i 两两互素
    n = prod(n_list); a = 0
    for ai, ni in zip(a_list, n_list):
        mi = n // ni
        _, inv, _ = EXTENDED_EUCLID(mi % ni, ni)   # mi^{-1} mod ni
        ci = mi * (inv % ni)
        a = (a + ai * ci) % n
    return a
# 算术运算 O(k lg n)（每个逆元一次扩展欧几里得）

推论 31.28：$n_i$ 两两互素，则方程组 $x\equiv a_i\pmod{n_i}$（$i=1..k$）模 $n$ 有唯一解。 推论 31.29：$n_i$ 两两互素，则对所有 $i$ 有 $x\equiv a\pmod{n_i}$ ⟺ $x\equiv a\pmod n$。

例：$a\equiv2\pmod5$，$a\equiv3\pmod{13}$，$n=65$。$13^{-1}\equiv2\pmod5$，$5^{-1}\equiv8\pmod{13}$，故 $c_1=13\cdot2=26$，$c_2=5\cdot8=40$，$a\equiv2\cdot26+3\cdot40=52+120\equiv42\pmod{65}$。图 31.3 用 $5\times13$ 的表展示对应：第 $i$ 行第 $j$ 列是满足 $a\bmod5=i$、$a\bmod13=j$ 的 $a$（模 65）；下移一行 $a$ 加 $c_1=26$，右移一列加 $c_2=40$，$a$ 加 1 对应沿对角线右下移动（越界回绕）。如第 4 行第 12 列为 64（即 $-1$）。

31.5 习题概览：31.5-1 解 $x\equiv4\pmod5$、$x\equiv5\pmod{11}$（答案 $x\equiv49\pmod{55}$）；31.5-2 求被 9、8、7 除余 1、2、3 的整数（$x\equiv10\pmod{504}$）；31.5-3 逆元在 CRT 下逐分量对应；31.5-4 多项式 $f(x)\equiv0\pmod n$ 的根数等于各 $\bmod n_i$ 根数之积。

31.6 元素的幂（Powers of an element）（PDF p.975–979）

考察 $a\in\mathbb Z_n^*$ 的幂序列 $a^0,a^1,a^2,\dots\pmod n$（31.33），$a^0\bmod n=1$。例：$3^i\bmod7$：$1,3,2,6,4,5,1,3,\dots$（周期 6）；$2^i\bmod7$：$1,2,4,1,2,4,\dots$（周期 3）。
本节 $\langle a\rangle$ 指 $\mathbb Z_n^*$ 中由乘法生成的子群，$\mathrm{ord}_n(a)$ 为 $a$ 模 $n$ 的阶。如 $\mathbb Z_7^*$ 中 $\langle2\rangle=\{1,2,4\}$，$\mathrm{ord}_7(2)=3$。

定理 31.30（欧拉定理）：$n>1$ 时，对所有 $a\in\mathbb Z_n^*$，$a^{\phi(n)}\equiv1\pmod n$。（推论 31.19 在 $\mathbb Z_n^*$ 上的翻译。） 定理 31.31（费马小定理）：$p$ 为素数时，对所有 $a\in\mathbb Z_p^*$，$a^{p-1}\equiv1\pmod p$。（因 $\phi(p)=p-1$。）对所有 $a\in\mathbb Z_p$（含 0）有 $a^p\equiv a\pmod p$。

原根/生成元（primitive root / generator）：若 $\mathrm{ord}_n(g)=|\mathbb Z_n^*|$，则 $\mathbb Z_n^*$ 中每个元素都是 $g$ 的幂。例：3 是模 7 的原根，2 不是。存在原根时称 $\mathbb Z_n^*$ 为循环群（cyclic）。
定理 31.32：$\mathbb Z_n^*$ 是循环群当且仅当 $n=2,4,p^e,2p^e$（$p>2$ 为素数，$e\ge1$）。（证明略。）
若 $g$ 为原根，对任意 $a\in\mathbb Z_n^*$ 存在 $z$ 使 $g^z\equiv a\pmod n$，称 $z$ 为 $a$ 模 $n$ 以 $g$ 为底的离散对数（discrete logarithm）或指标（index），记 $\mathrm{ind}_{n,g}(a)$。

定理 31.33（离散对数定理）：$g$ 为 $\mathbb Z_n^*$ 的原根，则 $g^x\equiv g^y\pmod n\iff x\equiv y\pmod{\phi(n)}$。证明：充分性由欧拉定理 $g^{y+k\phi(n)}\equiv g^y(g^{\phi(n)})^k\equiv g^y$；必要性由推论 31.18（幂序列以 $\phi(n)$ 为周期）。

定理 31.34：$p$ 为奇素数、$e\ge1$，则 $x^2\equiv1\pmod{p^e}$（31.34）只有 $x\equiv\pm1$ 两个解。证明：等价于 $p^e\mid(x-1)(x+1)$。$p>2$ 不能同时整除 $x-1$ 与 $x+1$（否则整除其差 2）。若 $p\nmid(x-1)$，则 $\gcd(p^e,x-1)=1$，由推论 31.5 得 $p^e\mid(x+1)$，即 $x\equiv-1$；对称地另一种情况 $x\equiv1$。

1 的非平凡平方根（nontrivial square root of 1 modulo $n$）：$x^2\equiv1\pmod n$ 但 $x\not\equiv\pm1$。例：$6^2=36\equiv1\pmod{35}$。

推论 31.35：若存在模 $n$ 的 1 的非平凡平方根，则 $n$ 是合数。证明：由定理 31.34 的逆否命题，$n$ 不是奇素数或奇素数幂；模 2 时 1 的平方根都是平凡的；且需 $n>1$。（这是 Miller–Rabin 测试的核心依据。）

反复平方法求模幂（repeated squaring）

求 $a^b\bmod n$（模取幂，modular exponentiation）是素性测试和 RSA 的基本操作。设 $b$ 的二进制为 $\langle b_k,b_{k-1},\dots,b_0\rangle$（$k+1$ 位，$b_k$ 为最高位）。

def MODULAR_EXPONENTIATION(a, b, n):
    c = 0                     # 仅用于证明的辅助变量
    d = 1
    for i in range(k, -1, -1):        # 从最高位到最低位
        c = 2 * c
        d = (d * d) % n               # 平方
        if bit(b, i) == 1:
            c = c + 1
            d = (d * a) % n           # 乘 a
    return d

循环不变式：每次迭代前，(1) $c$ 等于 $b$ 的二进制前缀 $\langle b_k,\dots,b_{i+1}\rangle$；(2) $d=a^c\bmod n$。初始化 $c=0,d=1=a^0$；保持：$b_i=0$ 时 $d'=d^2=a^{2c}$，$b_i=1$ 时 $d'=d^2a=a^{2c+1}$；终止时 $i=-1$，$c=b$，$d=a^b\bmod n$。
例（图 31.4）：$a=7$，$b=560=\langle1000110000\rangle$，$n=561$。各步 $c$：1,2,4,8,17,35,70,140,280,560；$d$：7,49,157,526,160,241,298,166,67,1。最终 $7^{560}\equiv1\pmod{561}$（561 是 Carmichael 数，后文用到）。
复杂度：输入为 $\beta$ 位数时，$O(\beta)$ 次算术运算、$O(\beta^3)$ 位操作；空间 $O(1)$ 个大整数。

31.6 习题概览：31.6-1 列出 $\mathbb Z_{11}^*$ 每个元素的阶，取最小原根 $g$（$g=2$）并列出 $\mathrm{ind}_{11,g}(x)$ 表；31.6-2 从右到左扫描 $b$ 的二进制位的模幂算法；31.6-3 已知 $\phi(n)$ 时用模幂求逆元：$a^{-1}\equiv a^{\phi(n)-1}\pmod n$。

31.7 RSA 公钥密码系统（The RSA public-key cryptosystem）（PDF p.979–986）

公钥密码系统（public-key cryptosystem）

作用：(1) 加密通信，窃听者无法解读；(2) 数字签名（digital signature）：电子版手写签名，任何人可验证、无人可伪造、消息改动任意一位即失效——同时认证签名者身份和消息内容，适用于电子合同、电子支票、电子订单等。
RSA 依赖"找大素数容易"与"分解两个大素数之积困难"之间的巨大差距。
每个参与者有公钥（public key）和私钥（secret key）；惯例用 Alice 与 Bob：$P_A,S_A$ 与 $P_B,S_B$。私钥保密，公钥可公开（如放在公共目录）。
设 $\mathcal D$ 为允许的消息集合（如所有有限长比特串）。最初、最简单的表述要求公钥、私钥确定 $\mathcal D$ 上的一一映射（置换）$P_A(\cdot)$、$S_A(\cdot)$，给定密钥可高效计算，且互为逆： $M=S_A(P_A(M)),\qquad M=P_A(S_A(M))\qquad(31.35,31.36)$
关键要求：除 Alice 外无人能在实际可行时间内计算 $S_A(\cdot)$，即便所有人都知道 $P_A$、能高效计算其逆函数 $P_A(\cdot)$。设计难点就是公开一个变换而不泄露其逆变换的算法。

加密流程（图 31.5）：Bob 取得 Alice 公钥 $P_A$；计算密文 $C=P_A(M)$ 发送；Alice 用 $S_A(C)=S_A(P_A(M))=M$ 解密。只有 Alice 能算 $S_A$，故只有她能读懂。

签名流程（图 31.6）：Alice 对消息 $M'$ 计算签名 $\sigma=S_A(M')$，发送 $(M',\sigma)$；Bob 用 Alice 公钥验证 $M'=P_A(\sigma)$（$M'$ 中应含 Alice 的名字以便知道用谁的公钥）。成立则接受，否则说明传输出错或是伪造。签名可被任何持有公钥者验证并转交第三方验证（如 Bob 把 Alice 签名的电子支票交给银行）。

签名消息不一定加密（可以明文）。签名并加密：先附签名，再用接收者公钥加密整个"消息+签名"；接收者用私钥解密后再用签名者公钥验签——相当于签了字的文件装进只有收件人能拆的信封。

RSA 系统

密钥生成：

随机选两个不同的大素数 $p,q$（例如各 1024 位）；
计算 $n=pq$；
选一个与 $\phi(n)=(p-1)(q-1)$ 互素的小奇数 $e$；
计算 $d=e^{-1}\bmod\phi(n)$（推论 31.26 保证存在且唯一，用 31.4 节扩展欧几里得求）；
公开 $P=(e,n)$ 作为公钥；
保密 $S=(d,n)$ 作为私钥。

消息域 $\mathcal D=\mathbb Z_n$：

P(M)=M^e\bmod n,\qquad(31.37)\qquad S(C)=C^d\bmod n.\qquad(31.38)

加密与签名都用这两个式子（签名时对消息施加私钥，验证时对签名施加公钥）。

复杂度：用 MODULAR-EXPONENTIATION。若 $\lg e=O(1)$，$\lg d\le\beta$，$\lg n\le\beta$，则施加公钥需 $O(1)$ 次模乘、$O(\beta^2)$ 位操作；施加私钥需 $O(\beta)$ 次模乘、$O(\beta^3)$ 位操作。

定理 31.36（RSA 正确性）：(31.37)(31.38) 定义 $\mathbb Z_n$ 上互逆的变换。证明：$P(S(M))=S(P(M))=M^{ed}\pmod n$。由 $ed=1+k(p-1)(q-1)$：若 $M\not\equiv0\pmod p$，$M^{ed}\equiv M(M^{p-1})^{k(q-1)}\equiv M\cdot1\equiv M\pmod p$（费马小定理）；若 $M\equiv0\pmod p$ 也显然成立。同理模 $q$ 成立。由中国剩余定理推论 31.29，$M^{ed}\equiv M\pmod n$ 对所有 $M$ 成立。

安全性与实用

安全性主要依赖大整数分解的困难：若能分解 $n$，就能像密钥创建者一样由 $p,q$ 推出私钥。反方向（"分解难⟹破解 RSA 难"）未被证明，但二十多年研究未找到比分解更容易的攻击。两个随机 1024 位素数相乘即可得到当前技术无法破解的公钥；需按推荐标准小心实现。写作时（2009）RSA 模数常用 768–2048 位。因此必须能高效找大素数（31.8 节）。
混合模式（hybrid / key-management mode）：RSA 很慢，故 Alice 随机选一个快速对称密码（加解密密钥相同）的短密钥 $K$，用 $K$ 加密长消息得 $C$，再用 Bob 的 RSA 公钥加密 $K$，发送 $(C,P_B(K))$。
哈希签名：结合抗碰撞哈希函数（collision-resistant hash function）$h$（易算，但找 $h(M)=h(M')$ 的两条消息在计算上不可行），$h(M)$ 是短"指纹"（如 256 位）。Alice 发送 $(M,S_A(h(M)))$；Bob 计算 $h(M)$ 并检查 $P_A(S_A(h(M)))=h(M)$。
证书（certificates）：受信任机构 $T$（公钥人人皆知）签发"Alice 的公钥是 $P_A$"的签名消息，Alice 随签名消息附上证书，收件人即可确认公钥归属。

31.7 习题概览：31.7-1 $p=11,q=29,n=319,e=3$，求 $d$（$\phi=280$，$d=187$）及 $M=100$ 的密文；31.7-2 若 $e=3$ 且敌手得到 $d$，可在多项式时间内分解 $n$（去掉 $e=3$ 的条件也成立，见 Miller）；31.7-3★ RSA 的乘法同态性 $P_A(M_1)P_A(M_2)\equiv P_A(M_1M_2)\pmod n$，并据此证明：能解密 1% 密文的过程可被随机化放大为高概率解密所有密文。

★31.8 素性测试（Primality testing）（PDF p.986–996）

素数的密度

素数分布函数（prime distribution function）$\pi(n)$：$\le n$ 的素数个数，如 $\pi(10)=4$。
定理 31.37（素数定理，prime number theorem）：$\lim_{n\to\infty}\dfrac{\pi(n)}{n/\ln n}=1$。即使 $n$ 较小也较准：$n=10^9$ 时 $\pi(n)=50{,}847{,}534$，$n/\ln n\approx48{,}254{,}942$，误差不到 6%。
随机选 $n$ 判素可看作伯努利试验，成功概率约 $1/\ln n$；由几何分布，期望约 $\ln n$ 次试验找到同长度的素数。例：找 1024 位素数约需测试 $\ln2^{1024}\approx710$ 个随机数（只选奇数可减半）。
记 $n=p_1^{e_1}\cdots p_r^{e_r}$（31.39），$n$ 为素数 ⟺ $r=1$ 且 $e_1=1$。
试除法（trial division）：用 $2,3,\dots,\lfloor\sqrt n\rfloor$ 去除，最坏 $\Theta(\sqrt n)=\Theta(2^{\beta/2})$，是输入长度的指数；只适合小 $n$ 或有小因子的 $n$，优点是能顺便给出一个因子。
本节只判断是否为素数，不求分解。判素远比分解容易，这有些出人意料。

伪素数测试（Pseudoprimality testing）

$\mathbb Z_n^+=\{1,\dots,n-1\}$；$n$ 为素数时 $\mathbb Z_n^+=\mathbb Z_n^*$。
$n$ 为合数且 $a^{n-1}\equiv1\pmod n$（31.40），称 $n$ 为以 $a$ 为基的伪素数（base-$a$ pseudoprime）。由费马小定理，若某个 $a$ 不满足 (31.40)，$n$ 必为合数。

def PSEUDOPRIME(n):          # n 为大于 2 的奇数
    if MODULAR_EXPONENTIATION(2, n - 1, n) != 1:
        return "COMPOSITE"    # 一定正确
    return "PRIME"            # 希望正确

只会犯一类错误：说合数一定对；说素数时仅当 $n$ 是基 2 伪素数时出错。$10{,}000$ 以下只错 22 个，前四个是 341、561、645、1105。随机 $\beta$ 位数上的出错概率随 $\beta\to\infty$ 趋于 0；按 Pomerance 估计，随机 512 位数被判素而实为基 2 伪素数的概率 $<10^{-20}$，1024 位 $<10^{-41}$。所以随机找大素数时几乎不会错；但输入非随机时需要更好的方法。
不能靠多试几个基彻底消除错误：存在 Carmichael 数——合数 $n$ 对所有 $a\in\mathbb Z_n^*$ 都满足 (31.40)。（$\gcd(a,n)>1$ 时 (31.40) 不成立，但若 $n$ 只有大素因子，很难碰到这样的 $a$。）前三个 Carmichael 数：561、1105、1729；$10^8$ 以下只有 255 个。

Miller–Rabin 随机化素性测试

两点改进：(1) 随机试多个基 $a$；(2) 在模幂的最后一串平方中寻找模 $n$ 的 1 的非平凡平方根，找到即判合数（推论 31.35）。

令 $n-1=2^tu$，$t\ge1$，$u$ 为奇数（$n-1$ 的二进制 = $u$ 的二进制后接 $t$ 个 0），于是 $a^{n-1}\equiv(a^u)^{2^t}$：先算 $a^u$，再连续平方 $t$ 次。

def WITNESS(a, n):            # True 表示 a 是 n 为合数的"证据"
    t, u = 分解 n-1 = 2^t * u（u 奇数, t ≥ 1）
    x_prev = MODULAR_EXPONENTIATION(a, u, n)       # x0 = a^u mod n
    for i in range(1, t + 1):
        x = x_prev * x_prev % n                   # x_i = x_{i-1}^2
        if x == 1 and x_prev != 1 and x_prev != n - 1:
            return True                           # 找到 1 的非平凡平方根
        x_prev = x
    if x_prev != 1:
        return True                               # a^{n-1} ≠ 1，费马测试失败
    return False

def MILLER_RABIN(n, s):       # n > 2 为奇数，s 为试验次数
    for j in range(s):
        a = RANDOM(1, n - 1)
        if WITNESS(a, n):
            return "COMPOSITE"   # 一定正确
    return "PRIME"               # 几乎肯定正确

序列满足 $x_i\equiv a^{2^iu}\pmod n$，$x_t\equiv a^{n-1}$。WITNESS 返回 TRUE 时，$a$ 连同返回原因（第 6 行：$x_{i-1}$ 是非平凡平方根；第 8 行：费马测试失败）构成 $n$ 为合数的证明。
用序列 $X=\langle x_0,\dots,x_t\rangle$ 描述（若中途某 $x_i=1$，后面都视为 1），四种情况：
1. $X=\langle\dots,d\rangle$，$d\ne1$：返回 TRUE（费马）；
2. $X=\langle1,1,\dots,1\rangle$：返回 FALSE，$a$ 不是证据；
3. $X=\langle\dots,-1,1,\dots,1\rangle$（末尾为 1，最后一个非 1 为 $-1$）：返回 FALSE；
4. $X=\langle\dots,d,1,\dots,1\rangle$，$d\ne\pm1$：返回 TRUE（$d$ 是非平凡平方根）。
例：$n=561$（Carmichael 数），$n-1=560=2^4\cdot35$，$t=4,u=35$。取 $a=7$，由图 31.4，$x_0=7^{35}\equiv241$，$X=\langle241,298,166,67,1\rangle$。最后一次平方发现 $7^{280}\equiv67$ 而 $7^{560}\equiv1$，67 是非平凡平方根，故返回 COMPOSITE。
复杂度：$\beta$ 位的 $n$ 需 $O(s\beta)$ 次算术运算、$O(s\beta^3)$ 位操作（不超过 $s$ 次模幂）。空间 $O(1)$ 个大整数。

Miller–Rabin 的错误率

与 PSEUDOPRIME 不同，出错概率不依赖于 $n$（没有坏输入），只取决于 $s$ 和抽样运气。

定理 31.38：若 $n$ 为奇合数，则 $n$ 为合数的证据数至少为 $(n-1)/2$。

证明：证明非证据（nonwitness）至多 $(n-1)/2$ 个。

非证据都在 $\mathbb Z_n^*$ 中：非证据 $a$ 满足 $a^{n-1}\equiv1$，即 $a\cdot a^{n-2}\equiv1$，方程 $ax\equiv1$ 有解，由推论 31.21 $\gcd(a,n)=1$。
再证非证据都落在 $\mathbb Z_n^*$ 的某个真子群 $B$ 中，由推论 31.16，$|B|\le|\mathbb Z_n^*|/2\le(n-1)/2$。
情形 1（$n$ 不是 Carmichael 数，实践中的主要情形）：存在 $x\in\mathbb Z_n^*$ 使 $x^{n-1}\not\equiv1$。取 $B=\{b\in\mathbb Z_n^*:b^{n-1}\equiv1\}$，非空（含 1）且乘法封闭，由定理 31.14 是子群；非证据都在 $B$ 中；$x\notin B$，故为真子群。
情形 2（$n$ 是 Carmichael 数，对所有 $x\in\mathbb Z_n^*$ 有 $x^{n-1}\equiv1$，31.41）：
- $n$ 不是素数幂：若 $n=p^e$（$e>1$，$p$ 奇素数），由定理 31.32 $\mathbb Z_n^*$ 循环，有生成元 $g$，$\mathrm{ord}_n(g)=\phi(n)=(p-1)p^{e-1}$；由 (31.41) 和离散对数定理得 $(p-1)p^{e-1}\mid p^e-1$，但左边被 $p$ 整除而右边不被 $p$ 整除，矛盾。
- 于是可把 $n$ 写成两个互素的、大于 1 的奇数之积 $n=n_1n_2$（如 $n_1=p_1^{e_1}$，$n_2$ 为其余部分）。
- 称 $(v,j)$ 为可接受对（acceptable pair）：$v\in\mathbb Z_n^*$，$j\in\{0..t\}$，$v^{2^ju}\equiv-1\pmod n$。$(n-1,0)$ 就是一个（$u$ 为奇数）。取使可接受对存在的最大 $j$，固定对应的 $v$。令 $B=\{x\in\mathbb Z_n^*:x^{2^ju}\equiv\pm1\pmod n\}$，乘法封闭，是子群。每个非证据都在 $B$ 中：非证据的序列 $X$ 要么全为 1，要么在不晚于第 $j$ 个位置出现 $-1$（由 $j$ 的最大性）。
- 构造 $w\in\mathbb Z_n^*-B$：由 $v^{2^ju}\equiv-1\pmod n$ 得 $\pmod{n_1}$ 也成立。由 CRT 取 $w\equiv v\pmod{n_1}$、$w\equiv1\pmod{n_2}$，则 $w^{2^ju}\equiv-1\pmod{n_1}$、$\equiv1\pmod{n_2}$，所以 $w^{2^ju}\not\equiv1$ 且 $\not\equiv-1\pmod n$，$w\notin B$。又 $\gcd(w,n_1)=\gcd(v,n_1)=1$、$\gcd(w,n_2)=1$，由定理 31.6 得 $w\in\mathbb Z_n^*$。故 $B$ 是真子群。

定理 31.39：对任意奇数 $n>2$ 和正整数 $s$，MILLER-RABIN$(n,s)$ 出错概率至多 $2^{-s}$。证明：$n$ 为合数时每轮至少以 $1/2$ 概率找到证据，$s$ 轮全部错过的概率 $\le2^{-s}$。$n$ 为素数时总是报告 PRIME。

结合先验概率（贝叶斯分析）：随机取 $\beta$ 位 $n$，事件 $A$ = "$n$ 为素数"，$\Pr\{A\}\approx1/\ln n\approx1.443/\beta$；事件 $B$ = "MILLER-RABIN 返回 PRIME"，$\Pr\{B\mid A\}=1$，$\Pr\{B\mid\bar A\}\le2^{-s}$。由贝叶斯公式

\Pr\{A\mid B\}=\frac{\Pr\{A\}\Pr\{B\mid A\}}{\Pr\{A\}\Pr\{B\mid A\}+\Pr\{\bar A\}\Pr\{B\mid\bar A\}}\approx\frac{1}{1+2^{-s}(\ln n-1)}.

$s$ 超过 $\lg(\ln n-1)$ 之前该概率不超过 1/2——需要这么多次试验才能抵消"$n$ 多半是合数"的先验偏向。1024 位数约需 $\lg(\beta/1.443)\approx9$ 次。实际中 $s=50$ 足以应付几乎任何应用。

实际情况更好：对随机奇合数，非证据的期望数目远小于 $(n-1)/2$，取 $s=3$ 就很少出错（未证明）。若 $n$ 不是随机选取的，改进版定理 31.38 能证明非证据至多 $(n-1)/4$，且这一界可达到。

31.8 习题概览：31.8-1 奇数 $n>1$ 既非素数也非素数幂时存在 1 的非平凡平方根；31.8-2★ 欧拉定理加强为 $a^{\lambda(n)}\equiv1$，$\lambda(n)=\mathrm{lcm}(\phi(p_1^{e_1}),\dots,\phi(p_r^{e_r}))$（Carmichael 函数，31.42），证明 $\lambda(n)\mid\phi(n)$；合数 $n$ 为 Carmichael 数 ⟺ $\lambda(n)\mid n-1$；最小者 $561=3\cdot11\cdot17$，$\lambda=\mathrm{lcm}(2,10,16)=80\mid560$；证明 Carmichael 数必无平方因子且至少是三个素数之积（故罕见）；31.8-3 若 $x$ 是模 $n$ 的 1 的非平凡平方根，则 $\gcd(x-1,n)$ 与 $\gcd(x+1,n)$ 都是 $n$ 的非平凡约数（由此 Miller–Rabin 的情形 4 可直接给出因子）。

★31.9 整数分解（Integer factorization）（PDF p.996–1001）

素性测试只告诉我们 $n$ 是合数，不给出因子。分解似乎难得多：现有超级计算机和最好算法也无法分解任意 1024 位数。

Pollard rho 启发式（Pollard's rho heuristic）

试除到 $R$ 可完全分解 $R^2$ 以内的数；同样工作量，POLLARD-RHO 可（若不太倒霉）分解 $R^4$ 以内的数。它只是启发式：运行时间和成功都不保证，但实践中非常有效；只用常数个存储单元（可用可编程计算器实现）。

def POLLARD_RHO(n):
    i = 1
    x = RANDOM(0, n - 1)        # x_1
    y = x
    k = 2
    while True:
        i += 1
        x = (x * x - 1) % n     # x_i = (x_{i-1}^2 - 1) mod n    (31.43)
        d = gcd(y - x, n)
        if d != 1 and d != n:
            print(d)            # 输出的一定是 n 的非平凡约数
        if i == k:
            y = x               # 保存下标为 2 的幂的 x：x_1, x_2, x_4, x_8, ...
            k = 2 * k

序列 $x_1,x_2,\dots$（31.44）只需保存最新值，故常数空间。$y$ 依次保存 $x_1,x_2,x_4,x_8,x_{16},\dots$，$k$ 总是下一个要保存的下标。
输出的数一定正确，但可能什么也不输出。期望在 $\Theta(\sqrt p)$ 次迭代后输出因子 $p$；$n$ 为合数时，除最大素因子外其余素因子都小于 $\sqrt n$，因此约 $n^{1/4}$ 次更新可完全分解 $n$。

分析：

$\mathbb Z_n$ 有限且每个值只依赖前一个，序列终将重复；一旦 $x_i=x_j$（$j<i$）就进入循环。把 $x_1..x_{j-1}$ 画成"尾巴"、$x_j..x_i$ 画成"圈"，形如希腊字母 ρ，故名。
假设 $f_n(x)=(x^2-1)\bmod n$ 表现得像随机函数（并非真随机，但与观测一致），由 5.4.1 节生日悖论，期望 $\Theta(\sqrt n)$ 步出现重复。
关键修改：取 $n$ 的非平凡因子 $p$ 且 $\gcd(p,n/p)=1$（例如 $p=p_1^{e_1}$；$e_1=1$ 时就是最小素因子）。令 $x'_i=x_i\bmod p$，则 $x'_{i+1}=((x_i^2-1)\bmod n)\bmod p=(x_i^2-1)\bmod p=((x'_i)^2-1)\bmod p=f_p(x'_i),$ 即模 $p$ 的序列服从同一递推（用到习题 31.1-7）。它期望 $\Theta(\sqrt p)$ 步重复，若 $p\ll n$，比模 $n$ 的序列早得多重复：只要两个 $x_i$ 模 $p$ 同余即可。
设 $t$ 为 $\langle x'_i\rangle$ 第一个重复值的下标，$u>0$ 为圈长，即 $x'_{t+i}=x'_{t+u+i}$（$i\ge0$）的最小 $t,u$，二者期望都是 $\Theta(\sqrt p)$。此时 $p\mid(x_{t+u+i}-x_{t+i})$，故 $\gcd(x_{t+u+i}-x_{t+i},n)>1$。一旦保存的 $y=x_k$ 满足 $k\ge t$，$y\bmod p$ 就在圈上；当 $k>u$ 后，算法绕圈一整圈而不改 $y$，于是 $x_i\equiv y\pmod p$ 时发现因子。期望步数 $\Theta(\sqrt p)$。发现的通常是 $p$，偶尔是 $p$ 的倍数。
图 31.7 例：$n=1387=19\cdot73$，$x_1=2$，序列 $x_1..x_{10}$ 为 2, 3, 8, 63, 1194, 1186, 177, 814, 996, 310, …（之后进入圈，最终回到 1186）；在 $x_7=177$ 时计算 $\gcd(63-177,1387)=19$ 得到因子 19（$y=x_4=63$）。第一个会重复的值是 1186，但在重复前已发现 19。(b) 模 19 的序列 2,3,8,6,16,8,…：$x_4=63$ 与 $x_7=177$ 模 19 都是 6；(c) 模 73 的序列。由 CRT，(a) 中每个节点对应 (b)(c) 中各一个节点。
两个可能的问题：(1) 启发式分析不严格，模 $p$ 的圈可能远大于 $\sqrt p$，此时结果正确但慢（实践中似乎不是问题）；(2) 可能只得到平凡因子 $n$：如 $n=pq$ 时 $p$ 与 $q$ 的 $t,u$ 恰好相同，两个因子在同一次 gcd 中一起出现。必要时换递推 $x_{i+1}=(x_i^2-c)\bmod n$ 重启（避免 $c=0$ 和 $c=2$）。
结论：它是找大数小素因子的首选方法。完全分解 $\beta$ 位合数只需找出所有 $<\lfloor n^{1/2}\rfloor$ 的素因子，期望至多 $n^{1/4}=2^{\beta/4}$ 次算术运算、$n^{1/4}\beta^2=2^{\beta/4}\beta^2$ 次位操作。最吸引人的是以 $\Theta(\sqrt p)$ 次运算找到小因子 $p$。

31.9 习题概览：31.9-1 在图 31.7(a) 的执行中，POLLARD-RHO 何时打印因子 73；31.9-2 给定 $f$ 与 $x_0$，精确求出 ρ 的尾长 $t$ 和圈长 $u$ 的高效算法（如 Floyd 判圈/Brent 算法）；31.9-3 找形如 $p^e$ 的因子期望需要多少步；31.9-4★ 批量 gcd：把连续若干个 $(y-x_i)$ 的乘积模 $n$ 累积后再做一次 gcd，说明实现、正确性及 $\beta$ 位 $n$ 的最佳批大小。

第 31 章思考题与章末注记（PDF p.1002–1005）

31-1 二进制 gcd 算法（binary gcd）：减法、奇偶测试、折半比取余快。(a) $a,b$ 都偶：$\gcd(a,b)=2\gcd(a/2,b/2)$；(b) $a$ 奇 $b$ 偶：$\gcd(a,b)=\gcd(a,b/2)$；(c) 都奇：$\gcd(a,b)=\gcd((a-b)/2,b)$；(d) 设计 $O(\lg a)$ 的算法（每种基本操作单位时间）。
31-2 欧几里得算法的位操作分析：(a) 长除法 $a\div b$ 需 $O((1+\lg q)\lg b)$ 位操作；(b) 定义 $\mu(a,b)=(1+\lg a)(1+\lg b)$，把 $\gcd(a,b)$ 归约为 $\gcd(b,a\bmod b)$ 的位操作至多 $c(\mu(a,b)-\mu(b,a\bmod b))$；(c) 因此 EUCLID 总共 $O(\mu(a,b))$，两个 $\beta$ 位输入时 $O(\beta^2)$（望远镜求和）。
31-3 斐波那契数的三种算法：(a) 朴素递归是指数时间；(b) 备忘录 $O(n)$；(c) 用矩阵 $\begin{pmatrix}0&1\\1&1\end{pmatrix}$ 的幂（反复平方）$O(\lg n)$；(d) 若加法 $\Theta(\beta)$、乘法 $\Theta(\beta^2)$（$F_n$ 有 $\Theta(n)$ 位），重新分析三者的运行时间。
31-4 二次剩余（quadratic residues）：$p$ 为奇素数，$a\in\mathbb Z_p^*$ 若 $x^2\equiv a\pmod p$ 有解则为二次剩余。(a) 恰有 $(p-1)/2$ 个；(b) 勒让德符号（Legendre symbol）$\left(\frac ap\right)\equiv a^{(p-1)/2}\pmod p$（欧拉判别法），据此高效判定；(c) $p=4k+3$ 时 $a^{k+1}\bmod p$ 是 $a$ 的平方根；(d) 随机找二次非剩余的算法（期望 2 次尝试）。

章末注记：入门书 Niven–Zuckerman；Knuth 第 2 卷讨论 gcd 等基本数论算法；Bach、Riesel、Dixon、Pomerance、Bach–Shallit 等综述计算数论。欧几里得算法见《几何原本》第 7 卷命题 1、2（约公元前 300 年），可能源自约公元前 375 年的 Eudoxus，或许是最古老的非平凡算法（只有古埃及乘法算法可与之相比）。中国剩余定理的特例归于孙子（约公元前 200 年至公元 200 年间），希腊的 Nichomachus 约公元 100 年给出同样特例，秦九韶（Chhin Chiu-Shao）1247 年推广，欧拉 1734 年给出完整陈述和证明。Miller–Rabin 测试来自 Miller 和 Rabin，是已知最快的随机化素性测试（常数因子内）；定理 31.39 的证明改编自 Bach，Monier 证明了更强结果。2002 年 Agrawal、Kayal、Saxena 给出确定性多项式时间素性测试（AKS），此前最快确定性算法（Cohen–Lenstra）为 $(\lg n)^{O(\lg\lg\lg n)}$；但实践中随机化测试仍更高效。找大"随机"素数见 Beauchemin 等。公钥密码概念来自 Diffie 与 Hellman；RSA 由 Rivest、Shamir、Adleman 于 1977 年提出。Goldwasser–Micali 证明随机化对安全公钥加密有效；Goldwasser–Micali–Rivest 给出伪造与分解同样困难的签名方案；Menezes 等综述应用密码学。rho 方法由 Pollard 发明，本书版本是 Brent 的变体。大数分解最好的算法运行时间约随数长度的立方根指数增长：一般数域筛法（general number-field sieve）估计为 $L(1/3,n)^{1.902+o(1)}$，其中 $L(\alpha,n)=e^{(\ln n)^\alpha(\ln\ln n)^{1-\alpha}}$；Lenstra 的椭圆曲线法（elliptic-curve method）像 rho 一样能快速找小因子 $p$，时间约 $L(1/2,p)^{\sqrt2+o(1)}$。

第 31 章本章要点

数论算法的输入规模按位数计，多项式时间指 $\mathrm{poly}(\lg a)$；大整数运算要按位操作计费（乘法 $\Theta(\beta^2)$）。
gcd 是最小正线性组合；欧几里得算法递归次数 $O(\lg b)$，相邻斐波那契数是最坏输入（Lamé 定理）；扩展欧几里得同时给出 Bézout 系数，用于求模逆。
$\mathbb Z_n$（加法群）、$\mathbb Z_n^*$（乘法群，大小 $\phi(n)$）；拉格朗日定理 ⟹ 元素阶整除群大小 ⟹ 欧拉定理/费马小定理。
$ax\equiv b\pmod n$ 有解 ⟺ $\gcd(a,n)\mid b$，有解时恰 $d$ 个，间隔 $n/d$。
中国剩余定理：$\mathbb Z_n\cong\mathbb Z_{n_1}\times\cdots\times\mathbb Z_{n_k}$，可构造性地重建（$c_i=m_i(m_i^{-1}\bmod n_i)$）。
反复平方模幂 $O(\beta)$ 次模乘；RSA：$ed\equiv1\pmod{\phi(n)}$，正确性靠费马小定理 + CRT，安全性依赖分解困难。
费马测试会被 Carmichael 数欺骗；Miller–Rabin 再检测 1 的非平凡平方根，每轮误判概率 ≤ 1/2，$s$ 轮 ≤ $2^{-s}$；结合素数定理的先验需做贝叶斯修正。
Pollard rho 利用生日悖论，期望 $\Theta(\sqrt p)$ 步找到因子 $p$，常数空间。

第 31 章与量化交易的关联

直接关联有限：数论算法本身不进入因子研究、风险模型或组合优化，如实说明。
系统实现与安全：交易系统与券商/交易所 API 的认证（RSA/ECDSA 签名、TLS 握手、HMAC 请求签名）、私钥管理与证书链都建立在本章内容之上；理解"签名 = 私钥作用于哈希"有助于正确实现 API 签名和排查验签失败。
随机数与模拟：线性同余生成器（LCG）$x_{i+1}=(ax_i+c)\bmod m$ 的周期分析依赖本章的模运算、阶和原根理论；蒙特卡洛定价和回测随机化需要知道生成器周期与质量（实务应使用 PCG、Mersenne Twister 或 Philox 等）。Pollard rho 中的 Floyd/Brent 判圈法也可用于检测伪随机序列或状态机的循环。
哈希与数据工程：模素数的哈希（如 Rabin–Karp 的滚动哈希，见第 32 章）、一致性哈希、布隆过滤器常用于行情去重、订单 ID 分片；CRT 思想可用于多模数哈希降低碰撞。
贝叶斯修正的思想：Miller–Rabin 的"先验概率 × 检验功效"分析与多重检验下的因子发现完全同构——在大量候选因子中真正有效的先验比例很低时，即便单次检验很严格，"通过检验"的因子仍可能多数是假阳性，需要更严格的阈值（类比需要 $s>\lg(\ln n-1)$ 次检验才抵消先验）。

第 31 章推荐习题

31.2-2、31.4-1、31.5-2：手算扩展欧几里得、模线性方程、中国剩余定理。
31.6-2、31.6-3：右到左模幂、用欧拉定理求逆元。
31.7-1：完整走一遍 RSA 密钥生成与加密。
31.8-3：非平凡平方根如何给出因子，理解 Miller–Rabin 与分解的联系。
思考题 31-1（二进制 gcd）与 31-3（矩阵快速幂求斐波那契）：快速幂思想在线性递推（如 AR 模型多步预测、马尔可夫链 $n$ 步转移）中通用。
31.9-2：判圈算法。

第 32 章字符串匹配（String Matching）（PDF p.1006–1029）

第 32 章引言（PDF p.1006–1009）

应用：文本编辑器查找单词、DNA 序列中查找模式、搜索引擎查找相关网页。
形式化：文本 $T[1..n]$，模式 $P[1..m]$（$m\le n$），元素取自有限字母表（alphabet）$\Sigma$（如 $\{0,1\}$ 或 $\{a..z\}$），字符数组称为字符串（strings）。
若 $0\le s\le n-m$ 且 $T[s+1..s+m]=P[1..m]$，称 $P$ 在 $T$ 中以偏移（shift）$s$ 出现（即从位置 $s+1$ 开始出现），$s$ 为有效偏移（valid shift），否则为无效偏移。字符串匹配问题：找出所有有效偏移。图 32.1：$P=abaa$，$T=abcabaabcabac$，只在 $s=3$ 出现。

图 32.2 各算法预处理与匹配时间：

算法	预处理时间	匹配时间
朴素算法（Naive）	0	$O((n-m+1)m)$
Rabin–Karp	$\Theta(m)$	$O((n-m+1)m)$（期望好得多）
有限自动机（Finite automaton）	$O(m\lvert\Sigma\rvert)$	$\Theta(n)$
Knuth–Morris–Pratt	$\Theta(m)$	$\Theta(n)$

总时间 = 预处理 + 匹配。Rabin–Karp 最坏不比朴素好，但平均和实际表现好得多，且易推广到其他模式匹配问题。

记号与术语：

$\Sigma^*$：$\Sigma$ 上所有有限长字符串的集合，含空串 $\varepsilon$。$|x|$ 为长度，$xy$ 为连接，长度 $|x|+|y|$。
前缀（prefix）$w\sqsubset x$：$x=wy$；后缀（suffix）$w\sqsupset x$：$x=yw$。二者都推出 $|w|\le|x|$。例：$ab\sqsubset abcca$，$cca\sqsupset abcca$。$\varepsilon$ 是任何串的前缀和后缀。$x\sqsupset y\iff xa\sqsupset ya$。两关系均传递。

引理 32.1（重叠后缀引理，overlapping-suffix lemma）：$x\sqsupset z$ 且 $y\sqsupset z$，若 $|x|\le|y|$ 则 $x\sqsupset y$；若 $|x|\ge|y|$ 则 $y\sqsupset x$；若相等则 $x=y$。（图 32.3 图示证明：二者都是 $z$ 的尾部，短的必是长的尾部。）

记 $P_k=P[1..k]$（$P_0=\varepsilon$，$P_m=P$），$T_k$ 同理。问题等价于找所有 $0\le s\le n-m$ 使 $P\sqsupset T_{s+m}$。
约定：等长字符串比较 "x == y" 为原语，从左到右比较、遇到不匹配即停，耗时 $\Theta(t+1)$，$t$ 为最长公共前缀长度。

32.1 朴素字符串匹配算法（The naive string-matching algorithm）（PDF p.1009–1011）

def NAIVE_STRING_MATCHER(T, P):
    n, m = len(T), len(P)
    for s in range(0, n - m + 1):
        if P[1..m] == T[s+1..s+m]:        # 隐式逐字符比较
            print("Pattern occurs with shift", s)

相当于把模式当"模板"在文本上滑动（图 32.4：$P=aab$，$T=acaabc$，试 $s=0,1,2,3$，在 $s=2$ 找到）。
时间 $O((n-m+1)m)$ 且最坏情况紧：$T=a^n$，$P=a^m$ 时每个偏移都比较 $m$ 次，$\Theta((n-m+1)m)$，$m=\lfloor n/2\rfloor$ 时为 $\Theta(n^2)$。无预处理，空间 $O(1)$。
低效原因：完全丢弃了一个偏移上获得的文本信息。例：$P=aaab$ 且 $s=0$ 有效，则 $s=1,2,3$ 都不可能有效（因 $T[4]=b$）。后续各节利用这类信息。

32.1 习题概览：32.1-1 写出 $P=0001$ 在 $T=000010001010001$ 上的比较过程；32.1-2 若 $P$ 中字符互不相同，可加速到 $O(n)$（失配时直接跳过已匹配部分）；32.1-3 $d$ 元字母表上随机 $P,T$，朴素算法期望比较次数为 $(n-m+1)\frac{1-d^{-m}}{1-d^{-1}}\le2(n-m+1)$，即随机串上朴素算法相当高效；32.1-4 模式含可匹配任意串（含空串）的间隙字符 $\diamond$（如 $ab\diamond ba\diamond c$ 在 $cabccbacbacab$ 中有两种出现方式），给出多项式时间判定算法（贪心地逐段匹配）。

32.2 Rabin–Karp 算法（PDF p.1011–1016）

实践表现好，并可推广到二维模式匹配等问题。预处理 $\Theta(m)$，最坏 $\Theta((n-m+1)m)$，在一定假设下平均更好。用到模同余等初等数论概念。
设 $\Sigma=\{0,..,9\}$（一般情形把字符看作 $d=|\Sigma|$ 进制数字），长度 $k$ 的串对应 $k$ 位十进制数，如 31415 对应 31,415。
$p$ = $P[1..m]$ 的十进制值；$t_s$ = $T[s+1..s+m]$ 的值（$s=0..n-m$）。$t_s=p\iff s$ 有效。若能 $\Theta(m)$ 算出 $p$、$\Theta(n-m+1)$ 算出所有 $t_s$，就能在 $\Theta(n)$ 内找出全部有效偏移。（脚注：写 $\Theta(n-m+1)$ 是因为 $m=n$ 时仍要 $\Theta(1)$。）
Horner 法则：$p=P[m]+10(P[m-1]+10(P[m-2]+\cdots+10(P[2]+10P[1])\cdots))$，$\Theta(m)$；$t_0$ 同理。
滚动更新： $t_{s+1}=10\,(t_s-10^{m-1}T[s+1])+T[s+m+1]\qquad(32.1)$ 减去最高位、左移一位、加入新的最低位。例：$m=5$，$t_s=31415$，新低位 2：$t_{s+1}=10(31415-10000\cdot3)+2=14152$。预计算 $10^{m-1}$（反复平方 $O(\lg m)$，直接算 $O(m)$ 也够），每次更新常数次运算。
数值过大问题：$m$ 位数不能假设单位时间运算。解决：对合适的模 $q$ 取模。选素数 $q$ 使 $10q$（一般为 $dq$）恰能放进一个机器字，全部用单精度运算： $t_{s+1}=\big(d(t_s-T[s+1]h)+T[s+m+1]\big)\bmod q,\qquad h\equiv d^{m-1}\pmod q.\qquad(32.2)$ $h$ 是 $m$ 位窗口最高位上数字 "1" 的值。
取模不完美：$t_s\equiv p\pmod q$ 不能推出 $t_s=p$；但 $t_s\not\equiv p$ 一定说明 $s$ 无效。所以把同余检验当作快速启发式过滤，同余时再显式比较 $P[1..m]=T[s+1..s+m]$ 排除伪命中（spurious hit）。$q$ 足够大时伪命中很少。

图 32.5 例：模 13。(a) 文本 2359023141526739921，阴影长 5 窗口 31415 的值模 13 为 7；(b) 对每个长 5 窗口算出模 13 的值，模式 $P=31415\equiv7\pmod{13}$，找到两个值为 7 的窗口：从文本位置 7 开始的是真匹配，从位置 13 开始的（67399）是伪命中；(c) 滚动计算：$31415\to14152$，模 13 下 $(7-3\cdot3)\cdot10+2\equiv8\pmod{13}$（$10^4\equiv3\pmod{13}$）。

def RABIN_KARP_MATCHER(T, P, d, q):
    n, m = len(T), len(P)
    h = pow(d, m - 1, q)
    p = 0; t = 0
    for i in range(1, m + 1):               # 预处理
        p = (d * p + P[i]) % q
        t = (d * t + T[i]) % q
    for s in range(0, n - m + 1):           # 匹配
        if p == t:                          # 命中（可能是伪命中）
            if P[1..m] == T[s+1..s+m]:
                print("Pattern occurs with shift", s)
        if s < n - m:
            t = (d * (t - T[s+1] * h) + T[s+m+1]) % q

不变式：每次执行第 10 行时 $t_s=T[s+1..s+m]\bmod q$。
复杂度：预处理 $\Theta(m)$；匹配最坏 $\Theta((n-m+1)m)$（如 $P=a^m,T=a^n$，每个偏移都有效都要验证）。空间 $O(1)$。
期望分析：许多应用中有效偏移只有常数 $c$ 个，期望匹配时间 $O((n-m+1)+cm)=O(n+m)$ 加上处理伪命中的时间。启发式假设"模 $q$ 像从 $\Sigma^*$ 到 $\mathbb Z_q$ 的随机映射"（参见 11.3.1 节除法散列；形式化可假设 $q$ 随机选取），则任一 $t_s\equiv p$ 的概率约 $1/q$，伪命中期望 $O(n/q)$ 个。期望匹配时间 $O(n)+O\big(m(v+n/q)\big),$ $v$ 为有效偏移数。若 $v=O(1)$ 且 $q\ge m$，则为 $O(n+m)=O(n)$。

32.2 习题概览：32.2-1 模 $q=11$ 在 $T=3141592653589793$ 中找 $P=26$ 会遇到多少伪命中；32.2-2 同时查找 $k$ 个模式（先设等长：把 $k$ 个哈希值放进哈希表；再推广到不等长）；32.2-3 在 $n\times n$ 字符阵列中找 $m\times m$ 模式（可平移不可旋转）——二维滚动哈希；32.2-4 Alice 与 Bob 比较两个 $n$ 位文件是否相同：选素数 $q>1000n$ 和随机 $x$，比较 $A(x)=\sum a_ix^i\bmod q$ 与 $B(x)$，不同文件碰撞概率 $\le1/1000$（多项式 $A-B$ 模 $q$ 至多 $n-1$ 个根，见 31.4-4）——多项式指纹（fingerprinting）。

32.3 利用有限自动机进行字符串匹配（String matching with finite automata）（PDF p.1016–1023）

字符串匹配自动机每个文本字符只检查一次、常数时间，匹配时间 $\Theta(n)$；但 $\Sigma$ 大时构造自动机代价大（32.4 节绕过此问题）。

有限自动机

有限自动机（finite automaton）$M=(Q,q_0,A,\Sigma,\delta)$：$Q$ 有限状态集；$q_0\in Q$ 初始状态；$A\subseteq Q$ 接受状态集；$\Sigma$ 有限输入字母表；$\delta:Q\times\Sigma\to Q$ 为转移函数（transition function）。

从 $q_0$ 开始逐个读字符，在状态 $q$ 读到 $a$ 就转移到 $\delta(q,a)$；当前状态属于 $A$ 时称已接受（accepted）已读串，否则拒绝（rejected）。
图 32.6：两状态自动机 $Q=\{0,1\}$，$q_0=0$，$\Sigma=\{a,b\}$，$\delta(0,a)=1,\delta(0,b)=0,\delta(1,a)=0,\delta(1,b)=0$，唯一接受态 1。它接受以奇数个 $a$ 结尾的串（$x=yz$，$y=\varepsilon$ 或以 $b$ 结尾，$z=a^k$，$k$ 奇）。输入 $abaaa$ 状态序列 $\langle0,1,0,1,0,1\rangle$，接受；$abbaa$ 为 $\langle0,1,0,0,1,0\rangle$，拒绝。
终态函数（final-state function）$\phi:\Sigma^*\to Q$：$\phi(\varepsilon)=q_0$，$\phi(wa)=\delta(\phi(w),a)$；$M$ 接受 $w\iff\phi(w)\in A$。

字符串匹配自动机

后缀函数（suffix function）$\sigma:\Sigma^*\to\{0..m\}$： $\sigma(x)=\max\{k:P_k\sqsupset x\},\qquad(32.3)$ 即 $P$ 中同时是 $x$ 后缀的最长前缀的长度。良定义因 $P_0=\varepsilon$ 是任何串的后缀。例：$P=ab$，$\sigma(\varepsilon)=0$，$\sigma(ccaca)=1$，$\sigma(ccab)=2$。$\sigma(x)=m\iff P\sqsupset x$；$x\sqsupset y\Rightarrow\sigma(x)\le\sigma(y)$。
模式 $P[1..m]$ 的自动机：状态集 $\{0..m\}$，初始 0，唯一接受态 $m$；转移函数 $\delta(q,a)=\sigma(P_qa).\qquad(32.4)$
直观：希望状态 $q$ 记录"$P$ 的、同时是已读文本 $T_i$ 后缀的最长前缀长度"，即维持不变式 $\phi(T_i)=\sigma(T_i)$（32.5）。读入 $a=T[i+1]$ 后应转到 $\sigma(T_ia)$；由于 $P_q$ 已概括了 $T_i$ 中有用的部分，$\sigma(T_ia)=\sigma(P_qa)$（引理 32.3）。两种情形：$a=P[q+1]$ 时沿"脊"（spine）前进到 $q+1$；否则要找更短的、同时是 $T_i$ 后缀的前缀——由预处理时模式与自身匹配得到。
例（图 32.7，$P=ababaca$）：$\delta(5,c)=6$（继续匹配）；$\delta(5,b)=4$，因 $P_5b=ababab$ 的最长"是 $P$ 前缀的后缀"为 $P_4=abab$。

图 32.7(b) 转移函数表（$P=ababaca$）：

状态	a	b	c	$P[q+1]$
0	1	0	0	a
1	1	2	0	b
2	3	0	0	a
3	1	4	0	b
4	5	0	0	a
5	1	4	6	c
6	7	0	0	a
7	1	2	0	—

图 32.7(c)：$T=abababacaba$，逐字符状态为 1,2,3,4,5,4,5,6,7,2,3，在位置 9 结束处找到一次出现（偏移 2）。

def FINITE_AUTOMATON_MATCHER(T, delta, m):
    n = len(T); q = 0
    for i in range(1, n + 1):
        q = delta[q][T[i]]
        if q == m:
            print("Pattern occurs with shift", i - m)
# 匹配时间 Θ(n)，额外空间 O(1)（不含转移表 O(m|Σ|)）

正确性

引理 32.2（后缀函数不等式）：对任意串 $x$ 和字符 $a$，$\sigma(xa)\le\sigma(x)+1$。证明：令 $r=\sigma(xa)$，$r=0$ 平凡；$r>0$ 时 $P_r\sqsupset xa$，去掉末尾 $a$ 得 $P_{r-1}\sqsupset x$，故 $r-1\le\sigma(x)$（图 32.8）。

引理 32.3（后缀函数递归引理）：若 $q=\sigma(x)$，则 $\sigma(xa)=\sigma(P_qa)$。证明：$P_q\sqsupset x$ ⟹ $P_qa\sqsupset xa$。令 $r=\sigma(xa)$，由引理 32.2，$r\le q+1=|P_qa|$。$P_r\sqsupset xa$、$P_qa\sqsupset xa$ 且 $|P_r|\le|P_qa|$，由引理 32.1 得 $P_r\sqsupset P_qa$，故 $r\le\sigma(P_qa)$；又 $P_qa\sqsupset xa$ 推出 $\sigma(P_qa)\le\sigma(xa)$。故相等（图 32.9）。

定理 32.4：对 $i=0..n$，$\phi(T_i)=\sigma(T_i)$。证明：对 $i$ 归纳。$i=0$ 时都为 0。设 $q=\phi(T_i)=\sigma(T_i)$，$a=T[i+1]$：$\phi(T_{i+1})=\phi(T_ia)=\delta(\phi(T_i),a)=\delta(q,a)=\sigma(P_qa)=\sigma(T_ia)=\sigma(T_{i+1})$（倒数第二步用引理 32.3 与归纳假设）。于是第 5 行 $q=m$ 当且仅当刚扫描完一次 $P$ 的出现，算法正确。

计算转移函数

def COMPUTE_TRANSITION_FUNCTION(P, Sigma):
    m = len(P)
    for q in range(0, m + 1):
        for a in Sigma:
            k = min(m + 1, q + 2)
            while True:                     # repeat ... until
                k = k - 1
                if P_k 是 P_q + a 的后缀:  # P_k ⊐ P_q a
                    break
            delta[q][a] = k
    return delta

从可能的最大值 $\min(m,q+1)$ 开始递减，直到 $P_k\sqsupset P_qa$（$k=0$ 时必成立）。
时间 $O(m^3|\Sigma|)$：外层 $m|\Sigma|$，repeat 至多 $m+1$ 次，每次后缀检验至多比较 $m$ 个字符。利用模式的前缀函数（习题 32.4-8）可降到 $O(m|\Sigma|)$。于是总体为 $O(m|\Sigma|)$ 预处理 + $\Theta(n)$ 匹配。空间 $O(m|\Sigma|)$。

32.3 习题概览：32.3-1 构造 $P=aabab$ 的自动机并在 $T=aaababaabaababaab$ 上演示；32.3-2 画 $P=ababbabbababbababbabb$ 在 $\{a,b\}$ 上的状态转移图；32.3-3 不可重叠模式（nonoverlappable，$P_k\sqsupset P_q$ 仅当 $k=0$ 或 $k=q$）的自动机形态；32.3-4★ 同时查找两个模式 $P,P'$ 的自动机，尽量减少状态数；32.3-5 含间隙字符的模式的自动机，$O(n)$ 匹配。

★32.4 Knuth–Morris–Pratt 算法（PDF p.1023–1033）

线性时间匹配，完全避免计算转移函数 $\delta$。只需辅助数组 $\pi[1..m]$（$\Theta(m)$ 时间预计算），匹配时可"即时"按摊还意义高效模拟 $\delta$。$\pi[q]$ 包含计算 $\delta(q,a)$ 所需的、与 $a$ 无关的信息；$\pi$ 只有 $m$ 项而 $\delta$ 有 $\Theta(m|\Sigma|)$ 项，预处理节省 $|\Sigma|$ 倍。

模式的前缀函数（prefix function）

动机（图 32.10）：$P=ababaca$ 在偏移 $s$ 处前 $q=5$ 个字符匹配、第 6 个失配。已匹配的 5 个字符决定了对应文本字符，因此 $s+1$ 必然无效（$P[1]=a$ 将对齐一个已知为 $b$ 的字符），而 $s'=s+2$ 让前三个模式字符与已知必然匹配的文本对齐。
一般问题：已知 $P[1..q]=T[s+1..s+q]$，求最小 $s'>s$ 使对某 $k<q$ 有 $P[1..k]=T[s'+1..s'+k]$，且 $s'+k=s+q$（32.6）。即求 $P_q$ 的、同时是 $T_{s+q}$ 后缀的最长真前缀 $P_k$。新偏移 $s'=s+(q-k)$；最好情况 $k=0$，一次跳过 $q-1$ 个偏移。新偏移处前 $k$ 个字符无需再比较。
由于 $T[s'+1..s'+k]$ 是已知文本、是 $P_q$ 的后缀，问题化为求最大的 $k<q$ 使 $P_k\sqsupset P_q$——只需模式与自身比较（图 32.10(c)：$P_5$ 的最长真后缀前缀为 $P_3$，$\pi[5]=3$）。
定义：$\pi:\{1..m\}\to\{0..m-1\}$， $\pi[q]=\max\{k:k<q\text{ 且 }P_k\sqsupset P_q\},$ 即 $P$ 中是 $P_q$ 真后缀的最长前缀的长度。

图 32.11(a)：$P=ababaca$ 的 $\pi$：

$i$	1	2	3	4	5	6	7
$P[i]$	a	b	a	b	a	c	a
$\pi[i]$	0	0	1	2	3	0	1

KMP 算法与前缀函数计算

def KMP_MATCHER(T, P):
    n, m = len(T), len(P)
    pi = COMPUTE_PREFIX_FUNCTION(P)
    q = 0                                   # 已匹配字符数
    for i in range(1, n + 1):               # 从左到右扫描文本
        while q > 0 and P[q+1] != T[i]:
            q = pi[q]                       # 下一个字符不匹配：回退
        if P[q+1] == T[i]:
            q = q + 1                       # 下一个字符匹配
        if q == m:                          # P 全部匹配
            print("Pattern occurs with shift", i - m)
            q = pi[q]                       # 继续寻找下一个匹配

def COMPUTE_PREFIX_FUNCTION(P):
    m = len(P)
    pi = [0] * (m + 1)
    pi[1] = 0
    k = 0
    for q in range(2, m + 1):
        while k > 0 and P[k+1] != P[q]:
            k = pi[k]
        if P[k+1] == P[q]:
            k = k + 1
        pi[q] = k
    return pi

两个过程结构相同：KMP-MATCHER 把 $T$ 与 $P$ 匹配，COMPUTE-PREFIX-FUNCTION 把 $P$ 与自身匹配。

运行时间（聚合分析）

COMPUTE-PREFIX-FUNCTION 为 $\Theta(m)$：难点是 while 循环（第 6–7 行）总共执行 $O(m)$ 次。$k$ 从 0 开始，只在第 9 行增加，每次 for 迭代至多加 1，总增量 $\le m-1$；进入循环时 $k<q$ 且 $q$ 每次递增，故始终 $k<q$，于是 $\pi[q]<q$，每次 while 迭代都使 $k$ 严格减小；$k$ 永不为负。因此总减量 ≤ 总增量 ≤ $m-1$，while 至多执行 $m-1$ 次。
同理 KMP-MATCHER 匹配时间 $\Theta(n)$（习题 32.4-4；也可用势函数，习题 32.4-5）。
与自动机相比：预处理从 $O(m|\Sigma|)$ 降到 $\Theta(m)$，匹配仍为 $\Theta(n)$。空间 $\Theta(m)$。

前缀函数计算的正确性

定义迭代 $\pi^{(0)}[q]=q$，$\pi^{(i)}[q]=\pi[\pi^{(i-1)}[q]]$，$\pi^*[q]=\{\pi[q],\pi^{(2)}[q],\dots,\pi^{(t)}[q]\}$，到 $\pi^{(t)}[q]=0$ 为止。

引理 32.5（前缀函数迭代引理）：对 $q=1..m$，$\pi^*[q]=\{k:k<q\text{ 且 }P_k\sqsupset P_q\}$。即反复迭代 $\pi$ 可以枚举出所有是 $P_q$ 真后缀的前缀。图 32.11(b)：$q=5$ 时 $\pi^*[5]=\{3,1,0\}$（$\pi[5]=3,\pi[3]=1,\pi[1]=0$），恰是把模板右滑时 $P_k$ 与 $P_5$ 某后缀匹配的 $k$。证明：(⊆) 对迭代次数 $u$ 归纳，利用 $\pi[i]<i$、$P_{\pi[i]}\sqsupset P_i$ 及 $\sqsupset$ 的传递性。(⊇) 反证：设差集非空，取其中最大元 $j$。因 $\pi[q]$ 是右边集合的最大元且在 $\pi^*[q]$ 中，$j<\pi[q]$；令 $j'$ 为 $\pi^*[q]$ 中大于 $j$ 的最小元。$P_j\sqsupset P_q$、$P_{j'}\sqsupset P_q$，由引理 32.1 得 $P_j\sqsupset P_{j'}$，且 $j$ 是小于 $j'$ 具有此性质的最大值，故 $\pi[j']=j$，从而 $j\in\pi^*[q]$，矛盾。

引理 32.6：若 $\pi[q]>0$，则 $\pi[q]-1\in\pi^*[q-1]$。证明：$r=\pi[q]$，$P_r\sqsupset P_q$，两端去掉末字符得 $P_{r-1}\sqsupset P_{q-1}$，由引理 32.5 得证。

对 $q=2..m$ 定义 $E_{q-1}=\{k\in\pi^*[q-1]:P[k+1]=P[q]\}=\{k:k<q-1\text{ 且 }P_{k+1}\sqsupset P_q\}$：那些能把 $P_k$ 扩展一个字符成为 $P_q$ 真后缀的 $k$。

推论 32.7：

\pi[q]=\begin{cases}0,&E_{q-1}=\emptyset\\1+\max\{k\in E_{q-1}\},&E_{q-1}\ne\emptyset\end{cases}

证明：非空时由定义 $\pi[q]\ge1+\max E_{q-1}$（32.8）；又令 $r=\pi[q]-1$，有 $P_{r+1}\sqsupset P_q$ 故 $P[r+1]=P[q]$，且由引理 32.6 $r\in\pi^*[q-1]$，于是 $r\in E_{q-1}$，得 $\pi[q]\le1+\max E_{q-1}$（32.9）。

算法正确性：每次 for 迭代开始时 $k=\pi[q-1]$。while 循环按降序遍历 $\pi^*[q-1]$ 直到找到 $P[k+1]=P[q]$，此时 $k=\max E_{q-1}$，置 $\pi[q]=k+1$；若找不到，$k$ 降为 0，第 8 行再检查 $P[1]=P[q]$ 决定 $\pi[q]$ 为 1 或 0。

KMP 正确性：模拟有限自动机

证明 KMP-MATCHER 第 $i$ 次迭代中、与 $m$ 比较时的 $q$ 与 FINITE-AUTOMATON-MATCHER 相同。
直观：若 $a=P[q+1]$，$\delta(q,a)=q+1$，KMP 不用 $\pi$ 直接加 1。否则 $0\le\delta(q,a)\le q$，KMP 的 while 循环按降序遍历 $\pi^*[q]$ 中的状态，停在某 $q'$ 使 $a=P[q'+1]$（然后转到 $q'+1$）或降到 0。
例（$P=ababaca$，$q=5$，$\pi^*[5]=\{3,1,0\}$）：读 $c$：直接到 6；读 $b$：while 执行一次到 $q'=\pi[5]=3$，$P[4]=b$ 匹配，到 4 $=\delta(5,b)$；读 $a$：$P[6]=c\ne a\to3$，$P[4]=b\ne a\to1$，$P[2]=b\ne a\to0$，再检查 $P[1]=a$，到 1 $=\delta(5,a)$。
形式化：对迭代次数归纳，设迭代开始时 $q'=\sigma(T_{i-1})$。三种情况：$\sigma(T_i)=0$（while 遍历完 $\pi^*[q']$ 都找不到，$q=0$）；$\sigma(T_i)=q'+1$（第一次 while 检验即失败，第 9 行加 1）；$0<\sigma(T_i)\le q'$（while 停在 $\pi^*[q']$ 中第一个满足 $P[q+1]=T[i]$ 的 $q$，此时 $q+1=\sigma(P_{q'}T[i])=\sigma(T_{i-1}T[i])=\sigma(T_i)$，用引理 32.3）。
第 12 行 $q=\pi[q]$ 必不可少：否则找到匹配后第 6 行会访问 $P[m+1]$。其正确性依赖 $\delta(m,a)=\delta(\pi[m],a)$（习题 32.4-8 的提示）。

32.4 习题概览：32.4-1 计算 $ababbabbabbababbabb$ 的前缀函数；32.4-2 $|\pi^*[q]|$ 的上界（$q$，如 $P=a^m$ 时紧）；32.4-3 通过串 $PT$ 的 $\pi$ 函数找出 $P$ 在 $T$ 中的出现（$\pi$ 值等于 $m$ 的位置，需注意 $\pi$ 不能超过 $m$，常在中间插入分隔符）；32.4-4 聚合分析证 KMP 匹配 $\Theta(n)$；32.4-5 势函数证同一结论；32.4-6 用 $\pi'$ 替换第 7 行的 $\pi$：$\pi'[q]=0$（若 $\pi[q]=0$），$=\pi'[\pi[q]]$（若 $P[\pi[q]+1]=P[q+1]$），$=\pi[q]$（否则），避免回退到注定失配的状态；32.4-7 线性时间判断 $T$ 是否为 $T'$ 的循环移位（如 arc 与 car）：在 $T'T'$ 中查找 $T$；32.4-8★ 利用 $\delta(q,a)=\delta(\pi[q],a)$（当 $q=m$ 或 $P[q+1]\ne a$）在 $O(m|\Sigma|)$ 内计算转移函数。

第 32 章思考题与章末注记（PDF p.1033–1034）

32-1 基于重复因子的匹配：$y^i$ 表示 $y$ 自身连接 $i$ 次，如 $(ab)^3=ababab$。若 $x=y^r$（$r>0$），称 $x$ 有重复因子 $r$；$\rho(x)$ 为最大的 $r$。(a) 计算所有 $\rho(P_i)$ 的高效算法；(b) $\rho^*(P)=\max_i\rho(P_i)$，随机二进制模式的期望 $\rho^*(P)=O(1)$；(c) 证明 Galil–Seiferas 的 REPETITION-MATCHER 在 $O(\rho^*(P)n+m)$ 内正确找到所有出现：失配或全匹配时偏移增加 $\max(1,\lceil q/k\rceil)$（$k=1+\rho^*(P)$），并重置 $q=0$。Galil–Seiferas 进一步得到只用 $O(1)$ 额外空间的线性时间算法。
章末注记：字符串匹配与有限自动机理论的关系见 Aho–Hopcroft–Ullman；KMP 由 Knuth 与 Pratt、以及 Morris 独立发明并联合发表；Reingold 等给出另一种讲法；Rabin–Karp 由 Karp 与 Rabin 提出；Galil–Seiferas 的常数额外空间线性算法。

第 32 章本章要点

字符串匹配：求所有有效偏移。朴素算法 $O((n-m+1)m)$，浪费了已匹配信息。
Rabin–Karp：把窗口视为 $d$ 进制数，模素数 $q$ 滚动哈希 $O(1)$ 更新；哈希相等再显式验证以排除伪命中；期望 $O(n+m)$，最坏 $\Theta((n-m+1)m)$。
有限自动机：状态 = 已匹配的最长前缀长度，$\delta(q,a)=\sigma(P_qa)$；匹配 $\Theta(n)$，预处理 $O(m|\Sigma|)$。核心引理：$\sigma(xa)=\sigma(P_{\sigma(x)}a)$。
KMP：前缀函数 $\pi[q]$ = $P_q$ 最长"真后缀即前缀"的长度；用 $\pi$ 摊还地模拟 $\delta$，预处理 $\Theta(m)$、匹配 $\Theta(n)$；聚合分析的关键是 $k$ 的总减量不超过总增量。

第 32 章与量化交易的关联

数据工程：日志、公告、新闻、监管文件的大规模关键字/模式检索（多模式匹配可用 Aho–Corasick，即 KMP 的多模式推广，或 Rabin–Karp 多模式哈希）是文本类另类数据管道的基础组件；Python 的 str.find、正则引擎内部也用类似思想。
滚动哈希与去重：Rabin–Karp 的滚动哈希用于行情/新闻的近重复检测（shingling + MinHash 的基础）、数据文件一致性校验（习题 32.2-4 的多项式指纹），以及增量同步（rsync 的滚动校验和）。
K 线形态的离散化匹配：把价格序列离散成符号串（如涨跌平符号化、SAX 表示）后，形态识别就是字符串匹配；自动机/KMP 能在流式行情上 $O(1)$ 每 tick 地识别预定形态。但需警惕：形态的统计显著性要单独检验，算法高效不等于信号有效。
流式处理：有限自动机只需常数状态、逐字符处理，非常适合低延迟的流式事件检测（如订单状态机、协议解析）。

第 32 章推荐习题

32.2-1：手算 Rabin–Karp 的伪命中，理解模 $q$ 选择的影响。
32.2-4：多项式指纹与随机化相等性检验。
32.3-1：手工构造字符串匹配自动机。
32.4-1、32.4-3、32.4-7：前缀函数的计算与灵活运用（循环移位检测）。
32.4-6、32.4-8：KMP 的改进以及由 $\pi$ 构造 $\delta$。

第 33 章计算几何（Computational Geometry）（PDF p.1035–1068）

第 33 章引言（PDF p.1035–1036）

计算几何研究解决几何问题的算法，应用于计算机图形学、机器人、VLSI 设计、CAD、分子建模、冶金、制造、纺织排样、林业、统计等。输入通常是几何对象集合（点集、线段集、按逆时针给出顶点的多边形），输出是查询回答（如是否有线段相交）或新几何对象（如点集的凸包，最小包围凸多边形）。
本章只讨论二维平面。对象用点集 $\{p_1,p_2,\dots\}$ 表示，$p_i=(x_i,y_i)$，$x_i,y_i\in\mathbb R$。$n$ 顶点多边形 $P$ 用边界上依次出现的顶点序列 $\langle p_0,\dots,p_{n-1}\rangle$ 表示。
安排：33.1 线段基本问题（顺/逆时针、转向、相交）；33.2 "扫描"技术，$O(n\lg n)$ 判断 $n$ 条线段中是否有相交；33.3 两个"旋转扫描"凸包算法：Graham 扫描 $O(n\lg n)$、Jarvis 步进 $O(nh)$（$h$ 为凸包顶点数）；33.4 $O(n\lg n)$ 分治求最近点对。

33.1 线段的性质（Line-segment properties）（PDF p.1036–1042）

凸组合（convex combination）：对不同点 $p_1,p_2$，$p_3=\alpha p_1+(1-\alpha)p_2$，$0\le\alpha\le1$，即直线 $p_1p_2$ 上位于两点之间（含端点）的点。线段 $\overline{p_1p_2}$ 是 $p_1,p_2$ 所有凸组合的集合，$p_1,p_2$ 为端点（endpoints）。考虑顺序时称有向线段 $\overrightarrow{p_1p_2}$；若 $p_1$ 为原点，可把它当作向量 $p_2$。
三个问题：(1) 两个共端点 $p_0$ 的有向线段 $\overrightarrow{p_0p_1}$ 是否在 $\overrightarrow{p_0p_2}$ 的顺时针方向？(2) 依次走 $\overline{p_0p_1}$、$\overline{p_1p_2}$，在 $p_1$ 处是否左转？(3) $\overline{p_1p_2}$ 与 $\overline{p_3p_4}$ 是否相交？
都能 $O(1)$ 回答，且只用加、减、乘和比较，不用除法和三角函数（代价高、易有舍入误差）。例如"求两直线方程 $y=mx+b$、求交点、再判断是否在线段上"的直接方法要用除法，线段接近平行时对除法精度非常敏感；叉积方法精确得多。

叉积（cross product）

$p_1\times p_2$ 可理解为由 $(0,0),p_1,p_2,p_1+p_2$ 构成的平行四边形的有向面积（图 33.1(a)）。更有用的定义： $p_1\times p_2=\det\begin{pmatrix}x_1&x_2\\y_1&y_2\end{pmatrix}=x_1y_2-x_2y_1=-p_2\times p_1.$
若 $p_1\times p_2>0$，则相对原点 $p_1$ 在 $p_2$ 的顺时针方向；$<0$ 则在逆时针方向（习题 33.1-1）；$=0$ 则共线（同向或反向）。脚注：叉积本是三维概念（按右手法则垂直于两向量、模为 $|x_1y_2-x_2y_1|$），本章只取该标量值。
共端点 $p_0$ 时平移：$(p_1-p_0)\times(p_2-p_0)=(x_1-x_0)(y_2-y_0)-(x_2-x_0)(y_1-y_0)$，正则 $\overrightarrow{p_0p_1}$ 在 $\overrightarrow{p_0p_2}$ 的顺时针方向，负则逆时针。

判断连续线段的转向

判断 $\angle p_0p_1p_2$ 的转向不必计算角度（图 33.2）：计算 $(p_2-p_0)\times(p_1-p_0)$。负 ⟹ $\overrightarrow{p_0p_2}$ 在 $\overrightarrow{p_0p_1}$ 逆时针方向，在 $p_1$ 左转；正 ⟹ 顺时针、右转；0 ⟹ 三点共线。

判断两线段是否相交

线段 $\overline{p_1p_2}$ 跨越（straddles）一条直线：$p_1$、$p_2$ 分别在直线两侧（端点恰在直线上是边界情况）。
两线段相交 ⟺ 以下至少一条成立：(1) 每条线段都跨越另一条所在直线；(2) 一条线段的某端点落在另一条线段上（边界情况）。

def DIRECTION(pi, pj, pk):
    return cross(pk - pi, pj - pi)        # (pk - pi) × (pj - pi)

def ON_SEGMENT(pi, pj, pk):               # 已知 pk 与 pi pj 共线
    return (min(xi, xj) <= xk <= max(xi, xj) and
            min(yi, yj) <= yk <= max(yi, yj))

def SEGMENTS_INTERSECT(p1, p2, p3, p4):
    d1 = DIRECTION(p3, p4, p1)
    d2 = DIRECTION(p3, p4, p2)
    d3 = DIRECTION(p1, p2, p3)
    d4 = DIRECTION(p1, p2, p4)
    if ((d1 > 0 and d2 < 0) or (d1 < 0 and d2 > 0)) and \
       ((d3 > 0 and d4 < 0) or (d3 < 0 and d4 > 0)):
        return True                       # 互相跨越
    elif d1 == 0 and ON_SEGMENT(p3, p4, p1): return True
    elif d2 == 0 and ON_SEGMENT(p3, p4, p2): return True
    elif d3 == 0 and ON_SEGMENT(p1, p2, p3): return True
    elif d4 == 0 and ON_SEGMENT(p1, p2, p4): return True
    else: return False
# 时间 O(1)，空间 O(1)

$d_1,d_2$ 异号 ⟺ $\overline{p_1p_2}$ 跨越 $\overline{p_3p_4}$ 所在直线；$d_3,d_4$ 异号同理。图 33.3：(a) 互相跨越，相交；(b) $\overline{p_3p_4}$ 跨越 $\overline{p_1p_2}$ 的直线但反之不然，不相交；(c) $p_3$ 与 $\overline{p_1p_2}$ 共线且在其间，相交（第 12 行返回）；(d) 共线但不在其间，不相交。若所有 $d_k\ne0$ 则不存在边界情况。$d_k=0$ 表示 $p_k$ 与另一线段共线，此时它在线段上 ⟺ 它位于两端点之间（ON-SEGMENT 用坐标包围盒判断）。

叉积的其他应用

33.3 节按相对原点的极角排序可用叉积比较（习题 33.1-3）；33.2 节红黑树维护线段的竖直次序时，不存显式键值，而用叉积判断与某竖直线相交的两条线段谁在上方。

33.1 习题概览：33.1-1 证明叉积符号与顺/逆时针的对应；33.1-2 ON-SEGMENT 只检 $x$ 坐标是错的（竖直线段反例）；33.1-3 用叉积比较、$O(n\lg n)$ 按极角排序（例：$(3,5)$ 相对 $(2,4)$ 极角 45°，$(3,3)$ 相对 $(2,4)$ 为 315°）；33.1-4 $O(n^2\lg n)$ 判断 $n$ 点中是否有三点共线（对每点按极角排序找相同角度）；33.1-5 多边形定义（分段线性闭曲线、边、顶点、简单多边形、内部/边界/外部、凸多边形：内部任两点连线仍在多边形内，凸多边形顶点不能写成两个不同点的凸组合）；判断点序列是否构成凸多边形——"所有连续转角不同时含左转和右转"线性但不总对（星形多边形可绕多圈），修正为同时检查总转角；33.1-6 $O(1)$ 判断从 $p_0$ 出发向右的水平射线是否与线段相交（化为线段相交）；33.1-7 射线法 $\Theta(n)$ 判断点是否在简单多边形内部（奇数次交点），注意射线过顶点和与边重合；33.1-8 $\Theta(n)$ 计算简单多边形面积（鞋带公式 $\frac12|\sum x_iy_{i+1}-x_{i+1}y_i|$）。

33.2 确定任意一对线段是否相交（Determining whether any pair of segments intersects）（PDF p.1042–1050）

只判断是否存在相交，不输出全部交点（最坏可能有 $\Theta(n^2)$ 个交点，习题 33.2-1），$O(n\lg n)$ 时间。
扫描（sweeping）技术：一条假想的竖直扫描线从左向右穿过几何对象，把 $x$ 维当作时间。扫描为几何对象提供一种排序（通常放进动态数据结构）并利用其关系。本算法按左到右顺序考察所有端点，每遇到一个端点就检查相交。
两个简化假设：没有竖直线段；没有三条线段交于一点。（习题 33.2-8、33.2-9 说明稍作修改即可去掉——处理边界条件往往是实现计算几何算法最难的部分。）

线段排序

线段 $s_1,s_2$ 在 $x$ 处可比较（comparable）：横坐标 $x$ 的扫描线与二者都相交。$s_1$ 在 $x$ 处位于 $s_2$ 之上，记 $s_1\succeq_x s_2$：二者可比较且 $s_1$ 与扫描线的交点更高，或二者在该扫描线上相交。图 33.4(a)：$a\succeq_rc$，$a\succeq_tb$，$b\succeq_tc$，$a\succeq_tc$，$b\succeq_uc$；$d$ 与其他线段都不可比较。
对给定 $x$，$\succeq_x$ 是与扫描线相交线段上的全预序（total preorder）：传递，且任两条相交线段至少一个方向成立（在扫描线上相交时两个方向都成立）；自反，但既不对称也不反对称。不同 $x$ 的全预序可能不同：左端点被扫到时线段进入，右端点时离开。
扫描线经过两线段交点时，二者在预序中交换位置（图 33.4(b)：交点左侧 $e\succeq_vf$，右侧 $f\succeq_we$）。因为没有三线共点，必存在某条扫描线 $z$ 使相交的 $e,f$ 在 $\succeq_z$ 中相邻（阴影区内的扫描线都如此）。

移动扫描线

扫描算法通常维护两组数据：(1) 扫描线状态（sweep-line status）：与扫描线相交对象之间的关系；(2) 事件点调度（event-point schedule）：按 $x$ 坐标从左到右排列的事件点，扫描线到达事件点时暂停处理，状态只在事件点改变。有些算法（如习题 33.2-7）事件点动态产生；本算法事先确定：每个端点都是事件点。
端点按 $x$ 升序；同 $x$（covertical）时左端点先于右端点，同类中 $y$ 小者优先。遇左端点插入线段、遇右端点删除线段；两线段首次相邻时检查是否相交。
扫描线状态是全预序 $T$，需支持：INSERT$(T,s)$、DELETE$(T,s)$、ABOVE$(T,s)$（紧邻其上的线段）、BELOW$(T,s)$。两线段在当前扫描线上相交时可能互为"上方"，在 $T$ 中次序任意。用红黑树实现，每个操作 $O(\lg n)$，键比较改为叉积比较（习题 33.2-2）。

def ANY_SEGMENTS_INTERSECT(S):
    T = 空的红黑树（全预序）
    将 S 中所有线段端点从左到右排序：按 (x, e, y) 字典序，e=0 为左端点、e=1 为右端点
    for p in 排序后的端点:
        if p 是线段 s 的左端点:
            INSERT(T, s)
            if (ABOVE(T, s) 存在且与 s 相交) or (BELOW(T, s) 存在且与 s 相交):
                return True
        if p 是线段 s 的右端点:
            if ABOVE(T, s) 与 BELOW(T, s) 都存在且二者相交:
                return True
            DELETE(T, s)
    return False

遇右端点时先检查 $s$ 上下两邻是否相交：若删除 $s$ 后它们会变成相邻，相交就必须在此时发现。边界：若 $p$ 落在另一线段 $s'$ 上，只要求 $s$ 与 $s'$ 被相邻地放入 $T$。
图 33.5：6 条线段 $a..f$，每条虚线为某事件点处的扫描线，下面列出处理后 $T$ 的次序；最右的扫描线在处理 $c$ 的右端点时，包围 $c$ 的 $d$ 与 $b$ 相交，返回 TRUE。

正确性：

返回 TRUE 只在真找到相交时（第 7、10 行），所以 TRUE 必正确。
反之，设存在相交，令 $p$ 为最左交点（并列取 $y$ 最小），$a,b$ 为交于 $p$ 的线段。$p$ 左侧无交点，故 $T$ 给出的次序在 $p$ 左侧处处正确。因无三线共点（脚注：否则可能有 $c$ 夹在 $a,b$ 之间，习题 33.2-8 处理），$a,b$ 在某扫描线 $z$（在 $p$ 左侧或过 $p$）上相邻；设 $q$ 为使它们变为相邻的事件点（$p$ 在 $z$ 上时 $q=p$）。由于按字典序处理事件点（$p$ 是最左交点中最低者；即使 $p$ 是 $a$ 的左端点又是 $b$ 的右端点，左端点事件先处理，$b$ 仍在 $T$ 中），处理 $q$ 之前 $T$ 的次序正确。若 $q$ 被处理，只有两种情况：$a$ 或 $b$ 被插入且另一个是其上/下邻（第 4–7 行发现）；或二者已在 $T$ 中、夹在中间的线段被删除使其相邻（第 8–11 行发现）。若 $q$ 未被处理，说明算法已提前找到相交返回 TRUE。

运行时间：第 2 行排序 $O(n\lg n)$（归并或堆排序）；for 循环至多 $2n$ 次，每次红黑树操作 $O(\lg n)$、相交测试 $O(1)$。总计 $O(n\lg n)$，空间 $O(n)$。

33.2 习题概览：33.2-1 $n$ 条线段可有 $\Theta(n^2)$ 个交点；33.2-2 $O(1)$ 判断可比较的两线段在 $x$ 处谁在上（不相交时用叉积；相交时也只用加减乘）；33.2-3 Mason 教授把"返回"改成"打印并继续"，是否能按从左到右打印出所有交点（不能：交换次序后的相邻关系没有被维护，且会漏交点）；33.2-4 $O(n\lg n)$ 判断 $n$ 顶点多边形是否简单；33.2-5 $O(n\lg n)$ 判断两个简单多边形是否相交；33.2-6 $O(n\lg n)$ 判断 $n$ 个圆盘中是否有两个相交；33.2-7 $O((n+k)\lg n)$ 输出全部 $k$ 个交点（Bentley–Ottmann，交点作为动态事件）；33.2-8 三条及以上线段共点时算法仍正确；33.2-9 把竖直线段的下端点当左端点、上端点当右端点即可处理竖直线段。

33.3 求凸包（Finding the convex hull）（PDF p.1050–1060）

点集 $Q$ 的凸包（convex hull）$CH(Q)$：使 $Q$ 中每点都在其边界或内部的最小凸多边形。假设点互不相同且至少有三点不共线。直观：把每个点看作木板上的钉子，凸包就是紧紧箍住所有钉子的橡皮筋形状（图 33.6：13 个点 $p_0..p_{12}$ 及其凸包）。
两个算法都按逆时针顺序输出凸包顶点：Graham 扫描（Graham's scan）$O(n\lg n)$；Jarvis 步进（Jarvis's march）$O(nh)$，$h$ 为凸包顶点数。凸包的每个顶点都是 $Q$ 中的点，两算法都利用这点决定保留/舍弃哪些点。
其他 $O(n\lg n)$ 方法：
- 增量法（incremental method）：先按 $x$ 排序，第 $i$ 步由前 $i-1$ 点的凸包加入第 $i$ 点（习题 33.3-6，$O(n\lg n)$）；
- 分治法（divide-and-conquer）：$\Theta(n)$ 分成左 $\lceil n/2\rceil$、右 $\lfloor n/2\rfloor$ 两半，递归求凸包，再用巧妙方法 $O(n)$ 合并，$T(n)=2T(n/2)+O(n)=O(n\lg n)$；
- 剪枝搜索法（prune-and-search）：类似 9.3 节线性时间中位数选择，反复丢弃常数比例的点，直到只剩凸包"上链"，再对下链做同样处理；渐近最快，$O(n\lg h)$。
凸包也是其他问题的第一步。例：平面最远点对（farthest-pair）必是凸包顶点（习题 33.3-3），而凸多边形的最远顶点对可 $O(n)$ 求出（旋转卡壳），所以最远点对总体 $O(n\lg n)$。

Graham 扫描

维护候选点栈 $S$：每个点入栈一次，不是凸包顶点的点最终出栈；结束时 $S$ 自底向顶恰为逆时针顺序的凸包顶点。TOP$(S)$ 返回栈顶、NEXT-TO-TOP$(S)$ 返回次栈顶（都不修改栈）。

def GRAHAM_SCAN(Q):                      # |Q| >= 3
    p0 = Q 中 y 最小的点（并列取最左者）
    p1..pm = 其余点按相对 p0 的极角逆时针排序
             （同极角者只保留离 p0 最远的一个）
    if m < 2:
        return "convex hull is empty"
    S = 空栈
    PUSH(p0, S); PUSH(p1, S); PUSH(p2, S)
    for i in range(3, m + 1):
        while NEXT_TO_TOP(S), TOP(S), p_i 构成的角不是左转:   # 用叉积判断
            POP(S)
        PUSH(p_i, S)
    return S

$p_0$（最低、最左）必是凸包顶点。极角排序用叉积比较（习题 33.1-3）；同极角的较近点是 $p_0$ 与最远点的凸组合，直接删去。极角都在 $[0,\pi)$ 内，所以排序即相对 $p_0$ 的逆时针顺序。$p_1$ 与 $p_m$ 也是凸包顶点（习题 33.3-1）。
沿凸包逆时针走每个顶点都应左转；遇到非左转（nonleft turn，右转或共线）就弹出栈顶。检查"非左转"而不只是右转，排除了凸包上的平角顶点（凸多边形顶点不能是其他顶点的凸组合）。
图 33.7：(a) 按极角编号的点与初始栈 $p_0,p_1,p_2$；(b)–(k) 每次 for 迭代后的栈，虚线为导致弹栈的非左转。例如 (h) 中 $\angle p_7p_8p_9$ 右转弹出 $p_8$，随后 $\angle p_6p_7p_9$ 右转弹出 $p_7$；(l) 为最终凸包，与图 33.6 一致。

定理 33.1（Graham 扫描的正确性）：$|Q|\ge3$ 时，终止时栈 $S$ 自底向顶恰为 $CH(Q)$ 的顶点，按逆时针顺序。证明：令 $Q_i=\{p_0,p_1,\dots,p_i\}$。被删去的同极角点不在 $CH(Q)$ 中，故 $CH(Q_m)=CH(Q)$；$p_0,p_1,p_i$ 都是 $CH(Q_i)$ 的顶点。循环不变式：每次 for 迭代开始时，$S$ 自底向顶恰为 $CH(Q_{i-1})$ 的顶点（逆时针）。

初始化：$S=\{p_0,p_1,p_2\}=Q_2$，三点构成自己的凸包且为逆时针。
保持：设 while 结束后栈顶为 $p_j$，其下为 $p_k$，此时 $S$ 与第 $j$ 次迭代后相同，即 $CH(Q_j)$。$p_i$ 的极角大于 $p_j$，且 $\angle p_kp_jp_i$ 左转（否则 $p_j$ 已被弹出），故压入 $p_i$ 后 $S$ 恰为 $CH(Q_j\cup\{p_i\})$（图 33.8(a)）。再证 $CH(Q_j\cup\{p_i\})=CH(Q_i)$：任何在第 $i$ 次迭代被弹出的 $p_t$（弹出时其下为 $p_r$）满足 $\angle p_rp_tp_i$ 非左转、$p_t$ 极角大于 $p_r$，所以 $p_t$ 位于三角形 $p_0p_rp_i$ 内部或边上（但不是顶点），不可能是 $CH(Q_i)$ 的顶点（图 33.8(b)），故 $CH(Q_i-\{p_t\})=CH(Q_i)$（33.1）。对第 $i$ 次迭代弹出的点集 $P_i$ 反复应用得 $CH(Q_i-P_i)=CH(Q_i)$，而 $Q_i-P_i=Q_j\cup\{p_i\}$。
终止：$i=m+1$，$S$ 为 $CH(Q_m)=CH(Q)$。

运行时间 $O(n\lg n)$（$n=|Q|$）：第 1 行 $\Theta(n)$；第 2 行用归并/堆排序加叉积比较 $O(n\lg n)$（删同极角点总共 $O(n)$）；初始化 $O(1)$；for 循环至多 $n-3$ 次，不计 while 共 $O(n)$。while 用聚合分析（同 17.1 节 MULTIPOP）：每点恰入栈一次，弹出次数不超过入栈次数，且 $p_0,p_1,p_m$ 永不弹出，至多 $m-2$ 次 POP，故 while 总计 $O(n)$。空间 $O(n)$。

Jarvis 步进

技术称为包装法（package wrapping / gift wrapping），$O(nh)$；当 $h=o(\lg n)$ 时渐近快于 Graham 扫描。直观：把纸的一端贴在最低点 $p_0$，向右拉紧再向上拉，直到碰到一个点——它也是凸包顶点；保持拉紧绕一圈回到 $p_0$。
形式化：构造凸包顶点序列 $H=\langle p_0,p_1,\dots,p_{h-1}\rangle$。$p_1$ 是相对 $p_0$ 极角最小的点（并列取最远者），$p_2$ 是相对 $p_1$ 极角最小的点……到达最高顶点 $p_k$（并列取最远）时得到右链（right chain）；然后从 $p_k$ 开始，以负 $x$ 轴为基准取最小极角构造左链（left chain），直到回到 $p_0$（图 33.9）。
也可一次绕完整个凸包（记录上一条边的角度，要求边角度在 $0$ 到 $2\pi$ 内严格递增），但分两条链的好处是不必显式计算角度，用 33.1 节的叉积比较即可。

def JARVIS_MARCH(Q):
    p0 = Q 中 y 最小（并列取最左）的点
    H = [p0]; cur = p0
    # 右链：以正 x 轴为基准
    while cur 不是最高点:
        nxt = argmin_{q in Q, q != cur} 相对 cur 的极角（叉积比较；并列取最远）
        H.append(nxt); cur = nxt
    # 左链：以负 x 轴为基准
    while True:
        nxt = argmin_{q in Q, q != cur} 相对 cur、从负 x 轴量起的极角
        if nxt == p0: break
        H.append(nxt); cur = nxt
    return H
# 每个凸包顶点做一次 O(n) 的求最小值，总时间 O(nh)，空间 O(n)

33.3 习题概览：33.3-1 证明 $p_1$ 与 $p_m$ 是凸包顶点；33.3-2 在支持加、比较、乘且排序下界为 $\Omega(n\lg n)$ 的模型中，按顺序输出凸包顶点也需 $\Omega(n\lg n)$（把排序归约为求抛物线 $y=x^2$ 上点的凸包）；33.3-3 最远点对必为凸包顶点；33.3-4 星形多边形（star-shaped polygon）：存在内部点 $p$ 位于边界上每一点 $q$ 的"影子"中（$\overline{qp}$ 完全在多边形内），所有这样的 $p$ 构成核（kernel）；给定逆时针顶点的 $n$ 顶点星形多边形，$O(n)$ 求凸包（图 33.10：(a) 星形；(b) 非星形，两点影子不相交、核为空）；33.3-5 在线凸包：逐点到达，每次输出当前凸包，总共 $O(n^2)$；33.3-6★ 增量法 $O(n\lg n)$ 实现。

33.4 寻找最近点对（Finding the closest pair of points）（PDF p.1060–1065）

在 $n\ge2$ 个点的集合 $Q$ 中找欧氏距离 $\sqrt{(x_1-x_2)^2+(y_1-y_2)^2}$ 最小的两点（可以重合，距离 0）。应用：空中/海上交通管制系统找最近的两个交通工具以检测潜在碰撞。
暴力法检查 $\binom n2=\Theta(n^2)$ 对。分治法满足 $T(n)=2T(n/2)+O(n)$，即 $O(n\lg n)$。

分治算法

每次递归输入子集 $P\subseteq Q$ 以及数组 $X$、$Y$，二者都包含 $P$ 的全部点，$X$ 按 $x$ 坐标单调递增排序，$Y$ 按 $y$ 坐标单调递增排序。注意不能在每次递归中重新排序，否则 $T(n)=2T(n/2)+O(n\lg n)=O(n\lg^2n)$（用习题 4.6-2 的主方法版本）；改用预排序（presorting）。

若 $|P|\le3$：暴力检查所有 $\binom{|P|}{2}$ 对。
分解：找一条竖直线 $l$ 把 $P$ 平分为 $P_L$（$\lceil|P|/2\rceil$ 个点，都在 $l$ 上或左侧）和 $P_R$（$\lfloor|P|/2\rfloor$ 个点，都在 $l$ 上或右侧）；把 $X$ 分为 $X_L,X_R$（仍按 $x$ 排序），把 $Y$ 分为 $Y_L,Y_R$（仍按 $y$ 排序）。
解决：分别对 $(P_L,X_L,Y_L)$、$(P_R,X_R,Y_R)$ 递归，得最近距离 $\delta_L,\delta_R$，令 $\delta=\min(\delta_L,\delta_R)$。
合并：最近点对要么是递归找到的距离为 $\delta$ 的点对，要么一点在 $P_L$、一点在 $P_R$。若后者距离小于 $\delta$，两点都在距 $l$ 不超过 $\delta$ 的范围内，即在以 $l$ 为中心、宽 $2\delta$ 的竖直带内（图 33.11(a)）。步骤：
1. 构造数组 $Y'$：从 $Y$ 中删去带外的点，仍按 $y$ 排序；
2. 对 $Y'$ 中每个点 $p$，只计算它与 $Y'$ 中后续 7 个点的距离，记录带内最近距离 $\delta'$；
3. 若 $\delta'<\delta$，返回带内的这一对及 $\delta'$；否则返回递归得到的点对及 $\delta$。

def CLOSEST_PAIR(P, X, Y):               # X 按 x 排序，Y 按 y 排序
    if len(P) <= 3:
        return 暴力求最近点对(P)
    mid = len(X) // 2 的上取整位置; l = X[mid].x
    PL = X 的前 ceil(|P|/2) 个点; PR = 其余
    XL, XR = 按 PL/PR 拆分 X（保持 x 有序）
    YL, YR = 线性扫描 Y，属于 PL 的追加到 YL，否则追加到 YR（保持 y 有序）
    dL = CLOSEST_PAIR(PL, XL, YL); dR = CLOSEST_PAIR(PR, XR, YR)
    d = min(dL, dR)
    Y1 = [p for p in Y if abs(p.x - l) < d]       # 2δ 宽的带，仍按 y 有序
    for i in range(len(Y1)):
        for j in range(i + 1, min(i + 8, len(Y1))):  # 只看后面 7 个点
            d = min(d, dist(Y1[i], Y1[j]))
    return d
# 主程序：先对 Q 按 x、按 y 各排序一次（预排序 O(n lg n)），再调用 CLOSEST_PAIR

正确性

递归在 $|P|\le3$ 时触底，保证不会出现只有一个点的子问题。
只需检查后续 7 个点：设某层最近点对为 $p_L\in P_L$、$p_R\in P_R$，距离 $\delta'<\delta$。二者到 $l$ 的距离都小于 $\delta$，竖直方向相距也小于 $\delta$，所以都在以 $l$ 为中心的 $\delta\times2\delta$ 矩形内。该矩形左半个 $\delta\times\delta$ 正方形中，$P_L$ 的点两两距离至少 $\delta$，至多容纳 4 个点（四个角，图 33.11(b)）；右半同理至多 4 个 $P_R$ 的点。故矩形内至多 8 个点（$l$ 上的点可能属于任一侧，最多可有 4 个点在 $l$ 上：矩形上、下边与 $l$ 的交点处各有一对重合点，每对一个属于 $P_L$、一个属于 $P_R$）。设 $p_L$ 在 $Y'$ 中位于 $p_R$ 之前，即使 $p_L$ 尽量靠前、$p_R$ 尽量靠后，$p_R$ 也在 $p_L$ 之后的 7 个位置之内。

实现与运行时间

关键：每次调用都要从一个有序数组中线性时间地得到有序子数组（$X_L,X_R,Y_L,Y_R,Y'$），相当于归并排序 MERGE 的逆操作——把一个有序数组拆成两个有序数组：

YL, YR = [], []
for p in Y:                      # 按 y 顺序扫描
    if p in PL: YL.append(p)     # 追加到末尾，保持有序
    else:       YR.append(p)

（判断 $p\in P_L$ 可在预处理时给每点打标记或比较 $(x,\text{id})$ 与分割点。）

预排序只在第一次递归前做一次，额外 $O(n\lg n)$；之后每层递归除子调用外线性时间。设 $T(n)$ 为递归部分时间、$T'(n)$ 为总时间：$T'(n)=T(n)+O(n\lg n)$，$T(n)=2T(n/2)+O(n)$（$n>3$），$T(n)=O(1)$（$n\le3$）。故 $T(n)=T'(n)=O(n\lg n)$，空间 $O(n)$（每层临时数组，可 $O(n\lg n)$ 若不复用）。

33.4 习题概览：33.4-1 Williams 教授提议把 $l$ 上的点都归入 $P_L$ 从而只看后续 5 个点，其缺陷（$l$ 上可能有很多点，破坏 $|P_L|=\lceil|P|/2\rceil$ 的平分，递归不再平衡）；33.4-2 实际只需检查后续 5 个位置；33.4-3 改用 $L_1$（曼哈顿）距离 $|x_1-x_2|+|y_1-y_2|$（一般 $L_m$ 距离为 $(|x_1-x_2|^m+|y_1-y_2|^m)^{1/m}$）；33.4-4 改用 $L_\infty$ 距离 $\max(|x_1-x_2|,|y_1-y_2|)$；33.4-5 有 $\Omega(n)$ 个点 $x$ 坐标相同时如何划分 $P_L,P_R$ 并判断 $Y$ 中点的归属，保持 $O(n\lg n)$；33.4-6 不预排序 $Y$，而在递归返回时归并 $Y_L,Y_R$ 得到有序 $Y$，仍 $O(n\lg n)$。

第 33 章思考题与章末注记（PDF p.1065–1068）

33-1 凸层（convex layers）：第 1 层是 $CH(Q)$ 的顶点，删去前 $i-1$ 层后剩余点 $Q_i$ 的凸包为第 $i$ 层。(a) $O(n^2)$ 算法；(b) 在排序需 $\Omega(n\lg n)$ 的模型中求凸层需 $\Omega(n\lg n)$。
33-2 极大层（maximal layers）：$(x,y)$支配（dominates）$(x',y')$ 若 $x\ge x'$ 且 $y\ge y'$；不被任何点支配的点称为极大点（maximal）。第 1 极大层 $L_1$ 为极大点集，删去前 $i-1$ 层后的极大点为 $L_i$。设共 $k$ 层，$y_i$ 为 $L_i$ 最左点的纵坐标（先假设坐标互不相同）。(a) $y_1>y_2>\cdots>y_k$；(b) 在所有点左侧加入新点 $(x,y)$：令 $j$ 为满足 $y_j<y$ 的最小下标（若 $y<y_k$ 则 $j=k+1$），则 $j\le k$ 时新点成为 $L_j$ 的新最左点，$j=k+1$ 时新增一层 $L_{k+1}=\{(x,y)\}$；(c) 从右向左扫描、二分查找，$O(n\lg n)$ 求所有极大层；(d) 坐标可相同时的问题与解决办法。（极大点即帕累托前沿。）
33-3 捉鬼敢死队与鬼：$n$ 个捉鬼者与 $n$ 个鬼配对，射线不能交叉（无三点共线）。(a) 存在一条过一个捉鬼者和一个鬼的直线，使其一侧捉鬼者数等于鬼数，$O(n\lg n)$ 找到（极角排序）；(b) $O(n^2\lg n)$ 的不交叉配对算法（递归）。
33-4 拾木棒：$n$ 根三维木棒（端点为 $(x,y,z)$，无竖直者），只能拾取上面没有其他木棒压着的木棒。(a) 判断 $a$ 在 $b$ 之上、之下或无关；(b) 判断能否全部拾起并给出合法顺序（构造"压在上面"关系的有向图做拓扑排序，有环则不行）。
33-5 稀疏凸包分布（sparse-hulled distributions）：$n$ 个点凸包期望大小为 $O(n^{1-\epsilon})$。例：单位圆盘内均匀分布 $\Theta(n^{1/3})$；$k$ 边凸多边形内均匀分布（$k$ 为常数）$\Theta(\lg n)$；二维正态分布 $\Theta(\sqrt{\lg n})$。(a) 两个分别有 $n_1,n_2$ 个顶点的凸多边形（可重叠），$O(n_1+n_2)$ 求全部点的凸包；(b) 对稀疏凸包分布，分治（前一半、后一半各求凸包再合并）得平均 $O(n)$ 时间。

章末注记：计算几何专著有 Preparata–Shamos、Edelsbrunner、O'Rourke。几何学古已有之，但几何算法的发展较新；Preparata 与 Shamos 指出最早的问题复杂度概念由 E. Lemoine 于 1902 年提出：他研究尺规作图，定义了五种基本操作（圆规一脚放在给定点上、放在给定线上、画圆、直尺边过给定点、画线），把完成某作图所需的操作数称为该作图的"简单度"（simplicity）。33.2 节算法来自 Shamos 与 Hoey。Graham 扫描原始版本由 Graham 给出，包装法由 Jarvis 提出。Yao 在决策树模型下证明任何凸包算法最坏需 $\Omega(n\lg n)$；考虑凸包顶点数 $h$ 时，Kirkpatrick–Seidel 的 $O(n\lg h)$ 剪枝搜索算法渐近最优。$O(n\lg n)$ 最近点对分治算法由 Shamos 提出，Preparata–Shamos 证明它在决策树模型下渐近最优。

第 33 章本章要点

叉积 $p_1\times p_2=x_1y_2-x_2y_1$ 是计算几何的核心原语：判断顺/逆时针、左/右转、线段相交，只用加减乘与比较，避免除法与三角函数的精度问题。
线段相交 = 互相跨越或端点落在另一线段上（共线边界情况用包围盒判断）。
扫描线技术：扫描线状态（红黑树维护全预序，比较用叉积）+ 事件点调度；只需在两线段首次相邻时检查，$O(n\lg n)$ 判定是否存在相交。
凸包：Graham 扫描（极角排序 + 栈 + 非左转弹栈，$O(n\lg n)$，聚合分析）；Jarvis 步进（礼品包装，$O(nh)$）；其他方法（增量、分治、剪枝搜索 $O(n\lg h)$）；下界 $\Omega(n\lg n)$。
最近点对分治：$2\delta$ 带内每点只需比较后续 7 个点；预排序 + 线性拆分有序数组，$O(n\lg n)$。
处理退化情况（共线、竖直、多线共点、重合点）是几何算法实现与证明中最难的部分。

第 33 章与量化交易的关联

组合优化与有效前沿：思考题 33-2 的"极大层"就是多目标（如高收益、低风险）下的帕累托前沿分层；在均值-方差平面上筛选非被支配的组合或策略可用 $O(n\lg n)$ 扫描。可行组合集合的凸性、有效前沿的凸包结构也与本章直接相关（例如在若干候选组合的（风险，收益）点集中，上凸包即可达的有效前沿近似）。
凸包与定价/风险：期权价格对执行价必须凸（否则存在蝶式套利），对报价点取下凸包可做无套利平滑；隐含波动率曲面的无套利修正、CVaR/风险预算中的分段线性凸函数也常用凸包。订单簿深度曲线、市场冲击成本函数的凸化同理。
最近点对/近邻：高维下最近点对的分治思想推广为 KD 树、近似最近邻（ANN），用于相似股票/相似行情片段检索、聚类和配对交易候选筛选（但金融中常为高维，平面算法不能直接套用）。
扫描线：区间重叠/事件排序问题（如检测交易时间窗口冲突、订单生命周期重叠统计、回测中按时间戳合并多源事件）本质上是一维扫描线，按"开始先于结束"等规则排序的细节与本章事件点排序规则一致。
数值稳健性：只用加减乘避免除法的思想，对低延迟系统中的定点数计算和跨平台结果一致性有借鉴意义。

第 33 章推荐习题

33.1-3：用叉积做极角排序（凸包与扫描算法的基础）。
33.1-7、33.1-8：点在多边形内判定、鞋带公式求面积。
33.2-3：理解扫描线为何只保证"是否存在"而不能直接列出所有交点。
33.3-2：凸包的 $\Omega(n\lg n)$ 下界（归约思想，与第 34 章呼应）。
33.4-2、33.4-6：最近点对的常数改进与免预排序版本。
思考题 33-2：极大层（帕累托前沿分层），与多目标策略筛选直接相关。

第 34 章 NP 完全性（NP-Completeness）（PDF p.1069 起；续见下一块）

第 34 章引言（PDF p.1069–1074）

前面几乎所有算法都是多项式时间算法：规模 $n$ 的输入最坏运行时间为 $O(n^k)$（$k$ 为常数）。并非所有问题都能多项式时间解决：有的问题（如图灵停机问题，Halting Problem）任何计算机都无法解决；有的可解但无法在任何 $O(n^k)$ 时间内解决。一般把多项式时间可解的问题视为易处理的（tractable）、"容易"，需要超多项式时间的视为难处理的（intractable）、"困难"。
本章研究一类状态未知的问题——NP 完全（NP-complete）问题：尚未找到任何一个的多项式时间算法，也无人能证明它们不存在多项式时间算法。这就是自 1971 年提出以来理论计算机科学最深刻的开放问题之一：P ≠ NP?
几对表面相似、难度却截然不同的问题（每对中一个属于 P、另一个 NP 完全）：
- 最短 vs. 最长简单路径：即使有负权边，单源最短路也能 $O(VE)$ 求出；但仅判断图中是否存在至少含给定边数的简单路径就是 NP 完全的。
- 欧拉回路 vs. 哈密顿回路：连通有向图的欧拉回路（Euler tour）恰好经过每条边一次（顶点可重复），$O(E)$ 可判定并求出（思考题 22-3）；哈密顿回路（hamiltonian cycle）是包含每个顶点的简单回路，判定有向图是否有哈密顿回路是 NP 完全的（本章后面证明无向图版本 NP 完全）。
- 2-CNF 可满足性 vs. 3-CNF 可满足性：布尔公式由 0/1 变量、联结词 ∧（AND）、∨（OR）、¬（NOT）和括号组成；若存在使其值为 1 的赋值，则称可满足（satisfiable）。$k$ 合取范式（$k$-CNF）是若干子句的 AND，每个子句是恰好 $k$ 个变量或其否定的 OR。例：$(x_1\lor\neg x_2)\land(\neg x_1\lor x_3)\land(\neg x_2\lor\neg x_3)$ 是 2-CNF，有满足赋值 $x_1=1,x_2=0,x_3=1$。2-CNF 可满足性多项式时间可判定，3-CNF 可满足性 NP 完全。

P、NP 与 NPC（非正式）

P 类：多项式时间可解的问题，即存在常数 $k$，可在 $O(n^k)$ 内求解（$n$ 为输入规模）。前面各章的大多数问题都在 P 中。
NP 类：多项式时间可验证的问题：若给出一个解的"证书"（certificate），能在输入规模的多项式时间内验证其正确。例：哈密顿回路问题的证书是 $|V|$ 个顶点的序列 $\langle v_1,\dots,v_{|V|}\rangle$，只需检查 $(v_i,v_{i+1})\in E$（$i=1..|V|-1$）且 $(v_{|V|},v_1)\in E$；3-CNF 可满足性的证书是一组变量赋值。
P ⊆ NP（P 中问题不需证书就能多项式时间解出）。开放问题是 P 是否为 NP 的真子集。
NPC 类（NP 完全）：属于 NP 且与 NP 中任何问题"一样难"。不加证明地给出：若任何一个 NP 完全问题能多项式时间求解，则 NP 中所有问题都有多项式时间算法。大多数理论计算机科学家相信 NP 完全问题是难处理的——如此多被深入研究的 NP 完全问题无一找到多项式算法，若它们全部可多项式求解将令人震惊；但证明其难处理的努力同样没有定论，因此不能排除其多项式可解的可能。
实践意义：若能证明一个问题 NP 完全，就为其难处理性提供了有力证据；工程上应转而设计近似算法（第 35 章）或求解易处理的特殊情形，而不是寻找精确快速算法。许多看起来不比排序、图搜索、网络流更难的自然问题其实是 NP 完全的。

证明 NP 完全性的思路概览

证明 NP 完全与本书其余部分的算法设计技术根本不同：它陈述的是问题有多难，不是证明高效算法存在，而是证明高效算法不太可能存在（类似 8.1 节比较排序 $\Omega(n\lg n)$ 下界，但技术不同于决策树）。依赖三个关键概念：

判定问题 vs. 最优化问题：最优化问题（optimization problem）中每个可行解有一个值，要找值最优的可行解。如 SHORTEST-PATH：给定无向图 $G$ 和顶点 $u,v$，求边数最少的 $u$–$v$ 路径（无权无向图的单对最短路）。NP 完全性只直接适用于判定问题（decision problem），答案只有"是/否"（1/0）。通常可给优化目标加一个界把最优化问题转成判定问题，如 PATH：给定图 $G$、顶点 $u,v$ 和整数 $k$，是否存在至多 $k$ 条边的 $u$–$v$ 路径？判定问题在某种意义上"更容易"（至少不更难）：解出 SHORTEST-PATH 后比较边数与 $k$ 即可回答 PATH。因此，若能证明判定问题困难，也就证明了对应最优化问题困难。
归约（reductions）：设要多项式时间解判定问题 $A$，问题的输入称为该问题的实例（instance）（如 PATH 的实例是特定的 $G,u,v,k$）。若已知判定问题 $B$ 可多项式时间求解，且有过程把 $A$ 的任意实例 $\alpha$ 变换为 $B$ 的实例 $\beta$，满足：(i) 变换耗时多项式；(ii) 答案相同（$\alpha$ 为"是"⟺ $\beta$ 为"是"）——称之为多项式时间归约算法（polynomial-time reduction algorithm）。则（图 34.1）：把 $\alpha$ 归约为 $\beta$，用 $B$ 的多项式算法判定 $\beta$，以其答案作为 $\alpha$ 的答案；三步都多项式，所以 $A$ 可多项式时间判定——用 $B$ 的"容易"证明 $A$ 的"容易"。NP 完全性则反过来用：若已知 $A$ 不存在多项式算法，且 $A$ 多项式归约到 $B$，则 $B$ 也不可能有多项式算法（否则 $A$ 就有了），这是反证法。对 NP 完全性，不能假设 $A$ 绝对没有多项式算法，但方法类似：在假定 $A$ 是 NP 完全的前提下证明 $B$ 是 NP 完全的。
第一个 NP 完全问题：归约需要一个已知 NP 完全的问题作起点。本章用电路可满足性（circuit-satisfiability）问题：给定由 AND、OR、NOT 门组成的布尔组合电路，是否存在一组布尔输入使输出为 1？34.3 节证明它是 NP 完全的。

章节安排：34.1 形式化"问题"的概念，定义多项式时间可解判定问题的复杂度类 P，并纳入形式语言框架；34.2 定义解可多项式验证的判定问题类 NP，正式提出 P ≠ NP 问题；34.3 用多项式时间归约联系问题、定义 NP 完全性，并概要证明电路可满足性 NP 完全；34.4 用归约方法更简单地证明其他问题 NP 完全，以两个公式可满足性问题为例；34.5 通过更多归约证明一系列其他问题 NP 完全。

34.1 多项式时间（Polynomial time）（PDF p.1074–1082）

为什么把多项式时间可解视为易处理（哲学而非数学理由）

虽然 $\Theta(n^{100})$ 的问题可以合理地视为难处理，但实际中很少有问题需要这么高次的多项式；经验表明一旦找到第一个多项式算法，往往很快会出现更高效的算法。
对许多合理的计算模型，在一个模型中多项式时间可解的问题在另一个模型中也多项式时间可解：如本书使用的串行随机访问机（RAM）与抽象图灵机（脚注：参见 Hopcroft–Ullman 或 Lewis–Papadimitriou）多项式可解的问题类相同；处理器数随输入规模多项式增长的并行计算机也相同。
多项式时间可解问题类有良好的封闭性：多项式在加法、乘法、复合下封闭。例如一个多项式时间算法的输出作为另一个的输入，复合算法仍为多项式；习题 34.1-5：调用常数次多项式时间子程序并另做多项式时间工作的算法仍是多项式时间的。

抽象问题（Abstract problems）

抽象问题 $Q$ 定义为问题实例集合 $I$ 与问题解集合 $S$ 上的二元关系。例：SHORTEST-PATH 的实例是（图，两个顶点）三元组，解是图中的顶点序列（空序列表示不存在路径）；问题本身是把每个实例关联到连接两顶点的最短路径的关系。最短路径不一定唯一，故一个实例可有多个解。
NP 完全性理论只关注判定问题，此时抽象判定问题可看作从实例集 $I$ 到解集 $\{0,1\}$ 的函数。例：对 PATH 的实例 $i=\langle G,u,v,k\rangle$，若 $u$ 到 $v$ 的最短路径至多 $k$ 条边则 PATH$(i)=1$，否则为 0。许多抽象问题是最优化问题而非判定问题，但通常可重述为不更难的判定问题。

编码（Encodings）

程序要求解抽象问题，必须用程序能理解的方式表示实例。抽象对象集合 $S$ 的编码（encoding）是从 $S$ 到二进制串集合的映射 $e$（脚注：任何至少两个符号的有限字母表上的串都可以）。例：自然数编码为 $\{0,1,10,11,100,\dots\}$，$e(17)=10001$；ASCII 中 A 编码为 1000001。复合对象（多边形、图、函数、有序对、程序）都可通过组合其组成部分的表示编码为二进制串。
实例集合为二进制串集合的问题称为具体问题（concrete problem）。若对长度 $n=|i|$ 的实例 $i$，算法能在 $O(T(n))$ 时间内给出解，则称它在 $O(T(n))$ 时间内解决该具体问题（脚注：假设输出与输入分开；每输出一位至少一步，所以输出规模为 $O(T(n))$）。若存在 $O(n^k)$ 算法（$k$ 为常数），称具体问题多项式时间可解（polynomial-time solvable）。
复杂度类 P：多项式时间可解的具体判定问题的集合。
编码 $e:I\to\{0,1\}^*$ 把抽象判定问题 $Q$ 诱导为具体判定问题 $e(Q)$（脚注：$\{0,1\}^*$ 表示所有 0/1 串的集合）：$e(i)$ 的解就是 $Q(i)$；不表示任何有意义实例的二进制串约定映射到 0。
我们希望多项式时间可解性与具体编码无关，但效率强烈依赖编码：设算法唯一输入为整数 $k$、运行时间 $\Theta(k)$。若 $k$ 用一元（unary，$k$ 个 1）表示，长度为 $n$ 时运行时间 $O(n)$，是多项式；若用自然的二进制表示，输入长度 $n=\lfloor\lg k\rfloor+1$，运行时间 $\Theta(k)=\Theta(2^n)$，是输入规模的指数。（这正是"伪多项式时间"的来源，例如背包问题的 $O(nW)$ 动态规划。）
实际中排除一元这类"昂贵"编码后，具体编码对多项式可解性影响很小：如三进制与二进制可在多项式时间互转。
函数 $f:\{0,1\}^*\to\{0,1\}^*$ 多项式时间可计算（polynomial-time computable）：存在多项式时间算法 $A$，对任意输入 $x$ 输出 $f(x)$。对实例集 $I$，两种编码 $e_1,e_2$ 多项式相关（polynomially related）：存在多项式时间可计算的 $f_{12},f_{21}$，对任意 $i\in I$ 有 $f_{12}(e_1(i))=e_2(i)$、$f_{21}(e_2(i))=e_1(i)$（脚注：还要求把非实例映射为非实例）。

引理 34.1：$Q$ 为实例集 $I$ 上的抽象判定问题，$e_1,e_2$ 为 $I$ 上多项式相关的编码，则 $e_1(Q)\in P\iff e_2(Q)\in P$。证明（只证正向，反向对称）：设 $e_1(Q)$ 可 $O(n^k)$ 求解，且由 $e_2(i)$ 计算 $e_1(i)$ 需 $O(n^c)$（$n=|e_2(i)|$）。求解 $e_2(Q)$：先算 $e_1(i)$，再对它运行 $e_1(Q)$ 的算法。转换耗时 $O(n^c)$，故 $|e_1(i)|=O(n^c)$（串行计算机的输出不会长于其运行时间）；求解耗时 $O(|e_1(i)|^k)=O(n^{ck})$，是多项式。

结论：实例用二进制还是三进制编码不影响问题的"复杂度"（是否多项式可解），但用一元编码可能会改变。此后默认实例用任何合理、简洁的编码：整数编码与其二进制表示多项式相关；有限集编码与"花括号括起、逗号分隔的元素列表"编码多项式相关（ASCII 即一例）。在此"标准"编码基础上可得元组、图、公式等的合理编码，用尖括号表示标准编码，如 $\langle G\rangle$ 表示图 $G$ 的标准编码。只要隐式使用与标准编码多项式相关的编码，就可直接谈论抽象问题，并通常忽略抽象问题与具体问题的区别；但须留意实际中标准编码不明显、编码确实有影响的问题。

形式语言框架（A formal-language framework）

字母表（alphabet）$\Sigma$：有限符号集。$\Sigma$ 上的语言（language）$L$：由 $\Sigma$ 中符号组成的任意串集合。例：$\Sigma=\{0,1\}$ 时 $L=\{10,11,101,111,1011,1101,10001,\dots\}$ 是素数的二进制表示构成的语言。空串 $\varepsilon$，空语言 $\emptyset$，$\Sigma$ 上所有串构成的语言 $\Sigma^*$（如 $\{0,1\}^*=\{\varepsilon,0,1,00,01,10,11,000,\dots\}$）。每个语言都是 $\Sigma^*$ 的子集。
语言运算：并、交（按集合论定义）；补 $\bar L=\Sigma^*-L$；连接 $L_1L_2=\{x_1x_2:x_1\in L_1,x_2\in L_2\}$；闭包或 Kleene 星 $L^*=\{\varepsilon\}\cup L\cup L^2\cup L^3\cup\cdots$（$L^k$ 为 $L$ 自身连接 $k$ 次）。
从语言角度，任何判定问题 $Q$ 的实例集就是 $\Sigma^*$（$\Sigma=\{0,1\}$）。$Q$ 完全由回答为 1 的实例刻画，所以可把 $Q$ 看作语言 $L=\{x\in\Sigma^*:Q(x)=1\}$。例：
$\text{PATH}=\{\langle G,u,v,k\rangle:G=(V,E)\text{ 为无向图},u,v\in V,k\ge0\text{ 为整数},G\text{ 中存在至多 }k\text{ 条边的 }u\text{–}v\text{ 路径}\}.$ （方便时同一名字既指判定问题又指对应语言。）
算法 $A$ 接受（accepts）串 $x$：输入 $x$ 时输出 $A(x)=1$；拒绝（rejects）：$A(x)=0$。$A$ 接受的语言 $L=\{x\in\{0,1\}^*:A(x)=1\}$。即使 $L$ 被 $A$ 接受，$A$ 对 $x\notin L$ 也不一定拒绝（可能永远循环）。若 $L$ 中每个串都被 $A$ 接受、不在 $L$ 中的每个串都被 $A$ 拒绝，称 $L$ 被 $A$ 判定（decided）。若 $L$ 被 $A$ 接受，且存在常数 $k$ 使任意长度为 $n$ 的 $x\in L$ 都在 $O(n^k)$ 时间内被接受，称 $L$ 被 $A$ 在多项式时间内接受；若存在常数 $k$ 使对任意长度 $n$ 的 $x\in\{0,1\}^*$，$A$ 都在 $O(n^k)$ 时间内正确判定 $x$ 是否属于 $L$，称 $L$ 被 $A$ 在多项式时间内判定。接受只要求对 $L$ 中的串给出答案，判定则必须对每个串正确地接受或拒绝。
例：PATH 可在多项式时间内被接受：验证 $G$ 编码一个无向图、$u,v$ 是其顶点，用广度优先搜索求 $u$ 到 $v$ 的最短路径并把边数与 $k$ 比较；满足则输出 1 停机，否则永远运行。它接受但不判定 PATH（对最短路超过 $k$ 条边的实例没有显式输出 0）。判定算法很容易设计：不满足时输出 0 停机（输入编码有误时也输出 0 停机）。而对图灵停机问题这样的问题，存在接受算法但不存在判定算法。
复杂度类（complexity class）非正式地说是一个语言集合，其成员资格由判定"给定串 $x$ 是否属于 $L$"的算法的某种复杂度度量（如运行时间）决定（严格定义更技术化，参见 Hartmanis–Stearns 的奠基论文）。
语言框架下 P 的另一定义：
$\mathrm P=\{L\subseteq\{0,1\}^*:\text{存在在多项式时间内判定 }L\text{ 的算法 }A\}.$

定理 34.2：$\mathrm P=\{L:L\text{ 被某个多项式时间算法接受}\}$。证明：多项式时间判定的语言类显然包含于多项式时间接受的语言类，只需证反向。设 $A$ 在 $O(n^k)$ 内接受 $L$，则存在常数 $c$ 使 $A$ 在至多 $cn^k$ 步内接受 $L$。构造 $A'$（经典的"模拟"论证）：对任意输入 $x$，模拟 $A$ 运行 $cn^k$ 步，若 $A$ 已接受则输出 1，否则输出 0 拒绝。模拟开销至多使运行时间增加多项式倍，所以 $A'$ 是判定 $L$ 的多项式时间算法。

注意此证明是非构造性的：对给定 $L\in\mathrm P$，我们可能并不知道接受 $L$ 的算法 $A$ 的运行时间界，但知道该界存在，从而存在能检查该界的 $A'$，尽管未必容易找到它。

34.1 习题概览：34.1-1 最优化问题 LONGEST-PATH-LENGTH（无向图两顶点间最长简单路径的边数）可多项式时间求解 ⟺ 判定问题 LONGEST-PATH（是否存在至少 $k$ 条边的简单路径）∈ P（二分/逐一询问 $k$）；34.1-2 给出无向图最长简单环问题的形式定义、相关判定问题及其语言；34.1-3 用邻接矩阵和邻接表把有向图编码为二进制串，并论证二者多项式相关；34.1-4 0-1 背包的动态规划算法（习题 16.2-2，$O(nW)$）是否多项式时间——否，$W$ 以二进制编码时它是输入长度的指数（伪多项式）；34.1-5 常数次调用多项式子程序加多项式额外工作仍是多项式，但多项式次调用可能导致指数时间（如每次调用使输出长度翻倍）；34.1-6 P 在并、交、连接、补、Kleene 星下封闭。

34.2 多项式时间验证（Polynomial-time verification）（PDF p.1082–1087）

研究"验证语言成员资格"的算法。对 PATH 实例 $\langle G,u,v,k\rangle$ 若同时给出一条 $u$–$v$ 路径 $p$，很容易检查 $p$ 是否为 $G$ 中路径、长度是否 $\le k$，$p$ 可看作实例属于 PATH 的"证书"。但 PATH ∈ P（甚至线性时间可解），证书帮助不大。下面看一个未知多项式判定算法、但给定证书后验证很容易的问题。

哈密顿回路

无向图 $G=(V,E)$ 的哈密顿回路（hamiltonian cycle）：包含 $V$ 中每个顶点的简单回路。含哈密顿回路的图称为哈密顿图（hamiltonian），否则为非哈密顿图（nonhamiltonian）。名称纪念 W. R. Hamilton，他描述了一个正十二面体上的数学游戏（图 34.2(a)）：一方在任意五个连续顶点插上钉子，另一方须补全路径形成包含所有顶点的回路（脚注引用 Hamilton 1856 年 10 月 17 日致 John T. Graves 的信，提到"Icosion"游戏：再插入 15 枚钉子循环覆盖其余各点，并终止于对手起点的旁边）。十二面体是哈密顿图；但并非所有图都是，如图 34.2(b) 中顶点数为奇数的二分图必为非哈密顿图（习题 34.2-2）。
形式语言：$\text{HAM-CYCLE}=\{\langle G\rangle:G\text{ 是哈密顿图}\}$。
朴素判定算法：列出顶点的所有排列逐一检查是否构成哈密顿回路。若用"合理"的邻接矩阵编码，顶点数 $m=\Omega(\sqrt n)$（$n=|\langle G\rangle|$ 为编码长度），共 $m!$ 种排列，运行时间 $\Omega(m!)=\Omega(\sqrt n!)=\Omega(2^{\sqrt n})$，不是任何 $O(n^k)$。事实上 HAM-CYCLE 是 NP 完全的（34.5 节证明）。

验证算法

稍简单的问题：朋友声称 $G$ 是哈密顿图，并给出沿回路依次的顶点来证明。验证很容易：检查所给序列是否为 $V$ 的一个排列，以及沿回路的每条相邻边是否都在图中，可在 $O(n^2)$ 时间完成。所以"图中存在哈密顿回路"的证明可在多项式时间内验证。
验证算法（verification algorithm）：两个参数的算法 $A$，一个是普通输入串 $x$，另一个是称为证书（certificate）的二进制串 $y$。若存在证书 $y$ 使 $A(x,y)=1$，称 $A$ 验证输入串 $x$。$A$ 验证的语言为 $L=\{x\in\{0,1\}^*:\text{存在 }y\in\{0,1\}^*\text{ 使 }A(x,y)=1\}.$
直观：对任意 $x\in L$，存在证书 $y$ 可让 $A$ 证明 $x\in L$；对任意 $x\notin L$，不存在能"证明" $x\in L$ 的证书。哈密顿回路问题中证书是某条哈密顿回路的顶点列表；非哈密顿图则没有任何顶点列表能骗过仔细检查"回路"的验证算法。

复杂度类 NP

NP：可被多项式时间算法验证的语言类（脚注：NP 意为"非确定性多项式时间"（nondeterministic polynomial time），最初在非确定性的背景下研究；本书用更简单但等价的"验证"概念；Hopcroft–Ullman 用非确定性计算模型讲述 NP 完全性）。精确地说，$L\in\mathrm{NP}$ 当且仅当存在两输入多项式时间算法 $A$ 和常数 $c$，使 $L=\{x\in\{0,1\}^*:\text{存在证书 }y,\ |y|=O(|x|^c),\ A(x,y)=1\}.$ 称 $A$ 在多项式时间内验证 $L$。注意证书长度必须是输入长度的多项式。
HAM-CYCLE ∈ NP（"知道一个重要集合非空总是好的"）。若 $L\in\mathrm P$ 则 $L\in\mathrm{NP}$：多项式判定算法可改成忽略证书的两参数验证算法。故 $\mathrm P\subseteq\mathrm{NP}$。
是否 P = NP 未知，但多数研究者相信不相等。直观上 P 是能被快速求解的问题，NP 是解能被快速验证的问题；经验上从零求解往往比验证清晰给出的解难得多（尤其在时间压力下），理论计算机科学家普遍相信这一类比可推广到 P 与 NP，即 NP 包含不在 P 中的语言。更有力（但非决定性）的证据是 NP 完全语言的存在（34.3 节）。
其他悬而未决的基本问题（图 34.3）：甚至不知道 NP 在补运算下是否封闭，即 $L\in\mathrm{NP}$ 是否推出 $\bar L\in\mathrm{NP}$。定义 co-NP $=\{L:\bar L\in\mathrm{NP}\}$，问题等价于是否 NP = co-NP。由于 P 对补封闭（习题 34.1-6），由习题 34.2-9 得 $\mathrm P\subseteq\mathrm{NP}\cap\text{co-NP}$；但也不知道 $\mathrm P=\mathrm{NP}\cap\text{co-NP}$ 是否成立，或 $\mathrm{NP}\cap\text{co-NP}-\mathrm P$ 是否非空。
图 34.3 四种可能（区域包含表示真子集）：(a) P = NP = co-NP，多数研究者认为最不可能；(b) NP 对补封闭，NP = co-NP，但未必 P = NP；(c) P = NP ∩ co-NP，但 NP 对补不封闭；(d) NP ≠ co-NP 且 P ≠ NP ∩ co-NP，多数研究者认为最可能。
结论：我们对 P 与 NP 精确关系的理解极不完整；但即便无法证明某问题难处理，若能证明它 NP 完全，也获得了关于它的宝贵信息。

34.2 习题概览：34.2-1 GRAPH-ISOMORPHISM $=\{\langle G_1,G_2\rangle:G_1,G_2\text{ 同构}\}\in\mathrm{NP}$（证书是顶点双射）；34.2-2 顶点数为奇数的无向二分图是非哈密顿图（回路在两侧交替，长度必为偶数）；34.2-3 若 HAM-CYCLE ∈ P，则按顺序列出哈密顿回路顶点也可多项式时间完成（逐条删边测试——判定到搜索的自归约）；34.2-4 NP 对并、交、连接、Kleene 星封闭，讨论补；34.2-5 NP 中任何语言都可在 $2^{O(n^k)}$ 时间内判定（枚举所有多项式长度的证书）；34.2-6 哈密顿路径（hamiltonian path，恰访问每个顶点一次的简单路径）语言 HAM-PATH $=\{\langle G,u,v\rangle\}\in\mathrm{NP}$；34.2-7 有向无环图上的哈密顿路径问题多项式可解（拓扑排序后检查相邻顶点间是否都有边）；34.2-8 由变量、¬、∧、∨ 和括号构成的公式若对所有赋值都为 1 则为永真式（tautology），TAUTOLOGY ∈ co-NP；34.2-9 证明 $\mathrm P\subseteq\text{co-NP}$；34.2-10 若 NP ≠ co-NP 则 P ≠ NP；34.2-11 连通无向图 $G$（至少 3 个顶点），把 $G$ 中距离至多 3 的顶点对全部连边得 $G^3$，证明 $G^3$ 是哈密顿图（构造生成树并归纳）。

34.3 NP 完全性与可归约性（NP-completeness and reducibility）（PDF p.1088–1093）（本块读至 PDF p.1093，电路可满足性部分续见下一块）

理论计算机科学家相信 P ≠ NP 的最有力理由来自 NP 完全问题类的存在：若任何一个 NP 完全问题能在多项式时间内解决，则 NP 中每个问题都有多项式时间解，即 P = NP；而多年研究从未为任何 NP 完全问题找到多项式算法。
HAM-CYCLE 就是 NP 完全问题：若能多项式时间判定 HAM-CYCLE，就能多项式时间解 NP 中每个问题；若 NP − P 非空，则可以肯定 HAM-CYCLE ∈ NP − P。NP 完全语言在某种意义上是 NP 中"最难"的语言。本节用精确的"多项式时间可归约性"比较语言的相对难度，正式定义 NP 完全语言，并概要证明 CIRCUIT-SAT 是 NP 完全的；34.4、34.5 节用归约证明更多问题 NP 完全。

可归约性（Reducibility）

直观：若问题 $Q$ 的任何实例都能"容易地改写"为 $Q'$ 的实例，且后者的解给出前者的解，则 $Q$ 可归约到 $Q'$，即 $Q$ "不比 $Q'$ 更难"。例：解一元线性方程 $ax+b=0$ 可归约到解二次方程：改写为 $0x^2+ax+b=0$。
多项式时间可归约（polynomial-time reducible）：语言 $L_1$ 多项式时间可归约到 $L_2$，记 $L_1\le_{\mathrm P}L_2$，若存在多项式时间可计算函数 $f:\{0,1\}^*\to\{0,1\}^*$，使对所有 $x\in\{0,1\}^*$， $x\in L_1\iff f(x)\in L_2.\qquad(34.1)$ $f$ 称为归约函数（reduction function），计算 $f$ 的多项式时间算法 $F$ 称为归约算法（reduction algorithm）。
图 34.4：两个语言都是 $\{0,1\}^*$ 的子集；$f$ 把 $L_1$ 内的点映到 $L_2$ 内，把 $L_1$ 外的点映到 $L_2$ 外。于是"$f(x)\in L_2$?"的答案直接就是"$x\in L_1$?"的答案。

引理 34.3：若 $L_1,L_2\subseteq\{0,1\}^*$ 满足 $L_1\le_{\mathrm P}L_2$，则 $L_2\in\mathrm P\Rightarrow L_1\in\mathrm P$。证明：设 $A_2$ 是多项式时间判定 $L_2$ 的算法，$F$ 是计算 $f$ 的多项式时间归约算法。构造 $A_1$（图 34.5）：对输入 $x$，用 $F$ 计算 $f(x)$，再用 $A_2$ 判定 $f(x)\in L_2$，把其输出作为自己的输出。正确性由 (34.1)；多项式时间由习题 34.1-5（复合）。

def A1(x):                 # 判定 L1
    y = F(x)               # 多项式时间归约：x ∈ L1 ⟺ y ∈ L2
    return A2(y)           # 多项式时间判定 L2
# 若 F 为 O(n^c)、A2 为 O(m^k)，则 |y| = O(n^c)，总时间 O(n^c + n^{ck})

NP 完全性

多项式时间归约给出"一个问题至少与另一个一样难（相差多项式因子）"的形式化手段：$L_1\le_{\mathrm P}L_2$ 表示 $L_1$ 至多比 $L_2$ 难一个多项式因子，"≤"记号因此便于记忆。
NP 完全（NP-complete）：语言 $L\subseteq\{0,1\}^*$ 满足
1. $L\in\mathrm{NP}$；
2. 对每个 $L'\in\mathrm{NP}$，$L'\le_{\mathrm P}L$。
只满足性质 2（不一定满足 1）的称为 NP 难（NP-hard）。NP 完全语言类记为 NPC。

定理 34.4：若任何一个 NP 完全问题多项式时间可解，则 P = NP。等价地，若 NP 中有任何问题不能多项式时间求解，则所有 NP 完全问题都不能多项式时间求解。证明：设 $L\in\mathrm P\cap\mathrm{NPC}$。对任意 $L'\in\mathrm{NP}$，由 NP 完全定义性质 2 有 $L'\le_{\mathrm P}L$，再由引理 34.3 得 $L'\in\mathrm P$。第二个陈述是第一个的逆否命题。

因此 P ≠ NP 问题的研究集中于 NP 完全问题。多数理论计算机科学家相信 P ≠ NP，对应图 34.6 的关系：P 与 NPC 都完全包含在 NP 中，且 $\mathrm P\cap\mathrm{NPC}=\emptyset$。但也许某天有人为某个 NP 完全问题找到多项式算法从而证明 P = NP；在此之前，证明一个问题 NP 完全就是它难处理的极好证据。

电路可满足性（Circuit satisfiability）（本块读到的部分）

至此只定义了 NP 完全，尚未证明任何问题 NP 完全；一旦证明一个，就可用多项式归约证明其他问题。第一个问题是电路可满足性。形式证明所需技术细节超出本书范围，书中给出依赖布尔组合电路基本知识的非正式证明。
布尔组合元件（boolean combinational element）：有常数个布尔输入和输出、执行确定函数的电路元件；布尔值取自 $\{0,1\}$，0 表示 FALSE，1 表示 TRUE。电路可满足性问题中用的元件计算简单布尔函数，称为逻辑门（logic gates）。图 34.7 的三种基本门及其真值表（truth table，给出每种输入组合下的输出）：
- NOT 门（反相器，inverter）：单输入 $x$，输出 $z=\neg x$（$0\to1$，$1\to0$）；
- AND 门：$z=x\land y$，仅 $x=y=1$ 时输出 1；
- OR 门：$z=x\lor y$，仅 $x=y=0$ 时输出 0，如 $0\lor1=1$。 AND、OR 可推广为多输入：AND 全为 1 时输出 1，OR 任一为 1 时输出 1。
布尔组合电路（boolean combinational circuit）：由导线互连的一个或多个布尔组合元件。导线可把一个元件的输出接到另一个元件的输入。一根导线至多连接一个元件输出，但可馈入多个元件输入，馈入的输入数称为该导线的扇出（fan-out）。没有元件输出连接的导线是电路输入（circuit input），接受外部输入值；没有元件输入连接的导线是电路输出（circuit output），把计算结果输出给外界（内部导线也可扇出到电路输出）。定义电路可满足性问题时，限定只有 1 个电路输出（实际硬件可以有多个）。
组合电路无环：为每个元件建一个顶点，对扇出为 $k$ 的导线建 $k$ 条有向边（导线把元件 $u$ 的输出接到元件 $v$ 的输入就有边 $(u,v)$），所得有向图必须无环。
图 34.8：两个只差一个门的电路。(a) 输入 $\langle x_1=1,x_2=1,x_3=0\rangle$ 使输出为 1（图中标出各导线取值），可满足；(b) 无论怎样赋值输出都是 0（习题 34.3-1），不可满足。
真值赋值（truth assignment）：一组布尔输入值。单输出电路可满足（satisfiable）：存在满足赋值（satisfying assignment），即使输出为 1 的真值赋值。
电路可满足性问题："给定由 AND、OR、NOT 门组成的布尔组合电路，它是否可满足？"为形式化，需约定电路的标准编码：电路规模（size）= 元件数 + 导线数；可设计类似图的编码，把电路 $C$ 映射为长度为规模多项式的二进制串 $\langle C\rangle$。于是 $\text{CIRCUIT-SAT}=\{\langle C\rangle:C\text{ 是可满足的布尔组合电路}\}.$
应用：计算机辅助硬件优化——若某子电路恒输出 0，它就是多余的，可用省去所有逻辑门、直接输出常数 0 的更简单子电路替换。因此人们希望该问题有多项式时间算法。
朴素方法：检查所有输入赋值；若电路有 $k$ 个输入，需检查多达 $2^k$ 种赋值……（续见下一块：当电路规模是 $k$ 的多项式时该方法为超多项式时间；引理 34.5 CIRCUIT-SAT ∈ NP、引理 34.6 CIRCUIT-SAT 是 NP 难，以及 34.4、34.5 节）

第 34 章本块覆盖部分小结（章末「本章要点」「与量化交易的关联」「推荐习题」由下一块在章末完成）

已覆盖：P、NP、NPC 的非正式与形式定义；判定问题 vs. 最优化问题；编码与多项式相关编码（一元编码导致的伪多项式陷阱）；形式语言框架下"接受"与"判定"的区别（定理 34.2：P 也是多项式时间可接受语言类）；验证算法与证书（证书长度须为多项式）；co-NP 与四种可能的类关系；多项式时间归约 $\le_{\mathrm P}$、引理 34.3、NP 完全与 NP 难的定义、定理 34.4；电路可满足性问题的定义。
与量化相关的提示（供下一块章末汇总参考）：组合构建中的基数约束（持仓数上限）、整手约束、最小交易单位使问题成为整数规划（NP 难），这正是第 29 章思考题 29-3 与本章的交汇点；34.1-4 提到的 0-1 背包伪多项式算法在"预算有限、按整数手数选股"时实际可用，因为权重用小整数编码时规模可控。

本块结束于 PDF 第 1093 页（原书第 1072 页），34.3 节"电路可满足性"讨论中途。

算法	预处理时间	匹配时间
朴素算法（Naive）	0	\(O((n-m+1)m)\)
Rabin–Karp	\(\Theta(m)\)	\(O((n-m+1)m)\)（期望好得多）
有限自动机（Finite automaton）	\(O(m\lvert\Sigma\rvert)\)	\(\Theta(n)\)
Knuth–Morris–Pratt	\(\Theta(m)\)	\(\Theta(n)\)

\(a\)	\(b\)	\(\lfloor a/b\rfloor\)	\(d\)	\(x\)	\(y\)
99	78	1	3	−11	14
78	21	3	3	3	−11
21	15	1	3	−2	3
15	6	2	3	1	−2
6	3	2	3	0	1
3	0	—	3	1	0

\(i\)	1	2	3	4	5	6	7
\(P[i]\)	a	b	a	b	a	c	a
\(\pi[i]\)	0	0	1	2	3	0	1