第 31 章 数论、字符串匹配与计算几何

本章合并原书第 31、32、33 章。这三章各自是一个独立领域：数论算法支撑着密码学（交易系统与交易所 API 的签名、TLS、密钥），字符串匹配是文本检索和流式事件识别的基础，计算几何处理点、线段、多边形。它们与因子研究、风险模型的距离比前几章远，所以本章只讲清"每个领域解决什么问题、核心算法是什么、为什么对、多快"，证明保留关键思路，细节注明回原书的位置。与量化关系最近的几处——随机数生成器的周期、Miller–Rabin 的贝叶斯分析与因子多重检验的同构、凸包与有效前沿、帕累托分层、期权报价的凸性修正——在量化实战中展开。

学习目标

读完本章，你应当能够：

1. 用扩展欧几里得算法求 gcd 与模逆元，用中国剩余定理在模数之间转换，用反复平方法在 \(O(\beta)\) 次模乘内计算模幂，并解释 RSA 加密与签名为什么正确。
2. 说清费马测试为何会被 Carmichael 数欺骗、Miller–Rabin 如何修补，以及"每轮误判率 \(\le1/2\)"与"通过检验后为素数的后验概率"之间的区别——并把它类比到因子挖掘的多重检验。
3. 写出 Rabin–Karp 的滚动哈希与 KMP 的前缀函数，说明二者的复杂度与适用场景，并能把 KMP 改写为逐 tick 处理的流式检测器。
4. 用叉积判断转向与线段相交，写出 Graham 扫描求凸包，理解扫描线和分治最近点对的思想。
5. 在量化场景中用凸包近似有效前沿、用极大层给策略做帕累托分层、用下凸包消除期权报价中的蝶式套利。

读前导读

这一章在解决什么问题

这一章是三个独立小领域的合集，与因子研究、风险模型的距离都比较远。对量化从业者，可以把它当作"工具箱说明书"来读，知道每样工具能干什么、什么时候该想起它，不必掌握全部证明。

数论研究整数的整除和"取余数"运算。它对量化的意义主要在工程层面：交易所和券商 API 的签名认证、TLS 加密连接都建立在 RSA 这类算法上；蒙特卡洛模拟用的伪随机数生成器，其周期长短也由数论决定。此外，素性测试中"通过检验后真是素数的概率"这个贝叶斯问题，和因子挖掘里"通过显著性检验的因子有多少是真的"是同一个公式。这一点你在 CFA 里学过基础率（base rate）和条件概率，读起来会很熟悉。

字符串匹配是在一长串字符里找某个模式出现的位置。把行情离散成"涨/跌/平"的符号串后，K 线形态识别就变成了字符串匹配；KMP 算法的价值在于它可以逐 tick 流式运行，只需记住一个整数状态。

计算几何处理平面上的点和线段。与金融最直接的联系是"凸"：有效前沿是凹曲线，看涨期权价格关于执行价必须是凸函数（否则存在蝶式套利），多目标比较策略时的帕累托前沿也是一个几何概念。31.10.4 节把这三件事用代码串了起来。

需要先想起来的数学

1. 取余与同余。 \(a\bmod n\) 是 \(a\) 除以 \(n\) 的余数，例如 \(17\bmod5=2\)。\(a\equiv b\pmod n\) 读作"\(a\) 与 \(b\) 模 \(n\) 同余"，意思是两者除以 \(n\) 余数相同，例如 \(17\equiv2\pmod5\)。日常例子是钟表：14 点就是下午 2 点，即 \(14\equiv2\pmod{12}\)。取余运算可以和加法、乘法交换顺序：先乘后取余，与先取余再乘再取余，结果相同。这是本章所有算法能在大整数上快速运行的原因。

2. 条件概率与贝叶斯公式。 \(\Pr\{A\mid B\}=\frac{\Pr\{B\mid A\}\Pr\{A\}}{\Pr\{B\}}\)。要点是不要把 \(\Pr\{\text{通过}\mid\text{合数}\}\) 和 \(\Pr\{\text{合数}\mid\text{通过}\}\) 混为一谈。见 第 00 册第 07 章 概率中的分析工具。

3. 二维行列式与凸集。 \(2\times2\) 行列式 \(\det\begin{pmatrix}a&b\\c&d\end{pmatrix}=ad-bc\)，几何上是两个向量张成的平行四边形的有向面积，见 第 00 册第 06 章 线性代数速成。"凸"的意思是：集合中任取两点，连线上的点也在集合内；凸函数的图像上任取两点，连线在图像上方。

怎么读这一章

建议的读法是：先通读三个部分开头的说明文字，了解每个领域在解决什么问题；然后直接读 31.10 量化实战，遇到不懂的算法再回到前面查。

值得细读的有：31.8 末尾的"贝叶斯修正"（与因子多重检验直接相关）、32.2 Rabin–Karp 的滚动哈希（思想与滚动窗口统计相同）、33.1 叉积、33.3 凸包和 33.5 极大层。

第一次可以跳过：31.3 群结构、31.4 模线性方程、31.5 中国剩余定理的细节、定理 31.38 的证明、31.9 Pollard rho、32.3 有限自动机、33.2 扫描线、33.4 最近点对。它们是各领域内部的经典内容，但离量化工作较远。

第一部分 数论算法（原书第 31 章）

31.1 输入规模与代价模型

数论算法里的"大输入"是指大整数，不是很多个数。输入规模按二进制位数计：输入为整数 \(a_1,\dots,a_k\) 的算法，若运行时间是 \(\lg a_1,\dots,\lg a_k\) 的多项式，才称为多项式时间算法。所以"从 2 试除到 \(\sqrt n\)"虽然只做 \(\sqrt n\) 次除法，却是 \(\Theta(2^{\beta/2})\)——位数 \(\beta\) 的指数。

大整数的算术不再是单位时间，改用位操作（bit operations）计数：普通方法乘两个 \(\beta\) 位整数要 \(\Theta(\beta^2)\)，除法取余同样 \(\Theta(\beta^2)\)；分治乘法（第 30 章思考题 30-1）为 \(\Theta(\beta^{\lg3})\)，Schönhage–Strassen 为 \(\Theta(\beta\lg\beta\lg\lg\beta)\)。本章同时按"算术运算次数"和"位操作数"分析。

31.2 整除与最大公约数

基本概念。\(d\mid a\)（\(d\) 整除 \(a\)）指存在整数 \(k\) 使 \(a=kd\)。\(a>1\) 且只有约数 1 和 \(a\) 的数称为素数（prime），否则为合数（composite）；1 既不是素数也不是合数。带余除法定理（division theorem）：对任意整数 \(a\) 和正整数 \(n\)，存在唯一的 \(q,r\) 使 \(a=qn+r\)、\(0\le r<n\)，记 \(r=a\bmod n\)。\(a\equiv b\pmod n\) 表示 \(a\bmod n=b\bmod n\)，模 \(n\) 的等价类记 \([a]_n\)，全体等价类构成 \(\mathbb Z_n=\{0,1,\dots,n-1\}\)。唯一分解定理：每个合数唯一地写成素数幂之积，如 \(6000=2^4\cdot3\cdot5^3\)。

最大公约数 \(\gcd(a,b)\) 是同时整除 \(a,b\) 的最大整数，约定 \(\gcd(0,0)=0\)。核心性质：

定理 31.2 \(a,b\) 不全为 0 时，\(\gcd(a,b)\) 是线性组合集合 \(\{ax+by:x,y\in\mathbb Z\}\) 中最小的正元素。

证明思路：设 \(s\) 是最小的正线性组合。\(a\bmod s=a-\lfloor a/s\rfloor s\) 也是线性组合且小于 \(s\)，只能为 0，所以 \(s\mid a\)；同理 \(s\mid b\)，故 \(s\le\gcd(a,b)\)。反过来 \(\gcd(a,b)\) 整除 \(a,b\)，所以整除 \(s\)，\(\gcd(a,b)\le s\)。\(\square\)

由此可推出：\(\gcd(a,b)=1\)（互素，relatively prime）当且仅当存在 \(x,y\) 使 \(ax+by=1\)；素数 \(p\mid ab\) 则 \(p\mid a\) 或 \(p\mid b\)。

欧几里得算法（Euclid's algorithm）依赖 GCD 递归定理（定理 31.9）：\(\gcd(a,b)=\gcd(b,a\bmod b)\)（\(b>0\)）。证明：\(a\bmod b\) 是 \(a,b\) 的线性组合，\(a\) 是 \(b\) 与 \(a\bmod b\) 的线性组合，两边的公约数集合相同。

def euclid(a, b):
    return a if b == 0 else euclid(b, a % b)

例：\(\mathrm{EUCLID}(30,21)\to(21,9)\to(9,3)\to(3,0)\)，结果 3。

运行时间与斐波那契数。引理 31.10：若 \(a>b\ge1\) 且递归 \(k\ge1\) 次，则 \(a\ge F_{k+2}\)、\(b\ge F_{k+1}\)（归纳，关键是 \(b+(a\bmod b)\le a\)）。由此得 Lamé 定理（定理 31.11）：若 \(b<F_{k+1}\)，递归少于 \(k\) 次。上界是紧的——相邻斐波那契数 \(\mathrm{EUCLID}(F_{k+1},F_k)\) 恰好递归 \(k-1\) 次，是最坏输入。因 \(F_k\approx\phi^k/\sqrt5\)（\(\phi\) 为黄金分割比），递归次数 \(O(\lg b)\)；对 \(\beta\) 位输入，\(O(\beta)\) 次算术运算，细致分析（思考题 31-2）为 \(O(\beta^2)\) 次位操作。

扩展欧几里得算法同时求出系数 \(x,y\) 使 \(d=\gcd(a,b)=ax+by\)（式 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)
    d, x, y = extended_euclid(b, a % b)
    return (d, y, x - (a // b) * y)

推导：递归返回 \(d=bx'+(a\bmod b)y'\)，代入 \(a\bmod b=a-\lfloor a/b\rfloor b\) 得 \(d=ay'+b(x'-\lfloor a/b\rfloor y')\)。原书图 31.1 的例子：

即 \(\gcd(99,78)=3=99\cdot(-11)+78\cdot14\)。

推导拆解：扩展欧几里得最常见的用途是求"模逆元"，即在模 \(n\) 的世界里做除法。以 \(5^{-1}\bmod11\) 为例：我们要找 \(x\) 使 \(5x\equiv1\pmod{11}\)。

先用欧几里得算法：\(11=2\cdot5+1\)，\(5=5\cdot1+0\)，所以 \(\gcd(5,11)=1\)。再把余数 1 倒推回去：\(1=11-2\cdot5\)，即 \(5\cdot(-2)+11\cdot1=1\)。两边对 11 取余，\(11\cdot1\) 这一项消失，得 \(5\cdot(-2)\equiv1\pmod{11}\)。\(-2\) 加上 11 变成正数 9，所以 \(5^{-1}\equiv9\pmod{11}\)。验算：\(5\times9=45=4\times11+1\)。

表格从下往上读就是这个"倒推"过程：每一层用下一层返回的 \((x',y')\)，按 \(x=y'\)、\(y=x'-\lfloor a/b\rfloor y'\) 算出本层的系数。RSA 生成私钥 \(d=e^{-1}\bmod\phi(n)\) 用的正是这一步。

31.3 模运算的群结构

原书用群论把"每步取模"的直觉严格化。要点如下，细节见原书 31.3 节。

• 群（group）\((S,\oplus)\)：封闭、有单位元、结合、每个元素有逆元；满足交换律的称阿贝尔群。
• 模 \(n\) 加法群 \((\mathbb Z_n,+_n)\)，大小 \(n\)。模 \(n\) 乘法群 \((\mathbb Z_n^*,\cdot_n)\)，\(\mathbb Z_n^*\) 是与 \(n\) 互素的剩余类，例如 \(\mathbb Z_{15}^*=\{1,2,4,7,8,11,13,14\}\)。乘法逆元由扩展欧几里得给出：\(\gcd(a,n)=1\) 时 \(ax+ny=1\)，于是 \(x\) 就是 \(a^{-1}\bmod n\)。例：\(5^{-1}\equiv9\pmod{11}\)，\(7^{-1}\equiv13\pmod{15}\)。
• 欧拉 \(\phi\) 函数 \(\phi(n)=|\mathbb Z_n^*|=n\prod_{p\mid n}(1-1/p)\)，如 \(\phi(45)=45\cdot\frac23\cdot\frac45=24\)；\(p\) 为素数时 \(\phi(p)=p-1\)。
• 拉格朗日定理（定理 31.15）：有限群的子群大小整除群的大小；所以真子群至多有一半元素（推论 31.16，Miller–Rabin 的分析要用）。
• 元素 \(a\) 的阶（order）\(\mathrm{ord}(a)\) 是使 \(a^t=e\) 的最小正整数 \(t\)，等于它生成的子群 \(\langle a\rangle\) 的大小；幂序列以 \(\mathrm{ord}(a)\) 为周期。由拉格朗日定理，\(\mathrm{ord}(a)\) 整除群的大小。

31.4 模线性方程

求解 \(ax\equiv b\pmod n\)。设 \(d=\gcd(a,n)\)，由扩展欧几里得得 \(ax'+ny'=d\)。

• 有解当且仅当 \(d\mid b\)（推论 31.21）；
• 有解时恰有 \(d\) 个模 \(n\) 的解：\(x_0=x'(b/d)\bmod n\)，\(x_i=x_0+i(n/d)\)，\(i=0,\dots,d-1\)（定理 31.23、31.24）；
• \(\gcd(a,n)=1\) 时解唯一；特别地 \(ax\equiv1\) 的解就是逆元（推论 31.26）。

例：\(14x\equiv30\pmod{100}\)。\(\mathrm{EXTENDED\text{-}EUCLID}(14,100)=(2,-7,1)\)，\(2\mid30\)，\(x_0=(-7)(15)\bmod100=95\)，两个解为 95 与 45。复杂度 \(O(\lg n+\gcd(a,n))\) 次算术运算。

31.5 中国剩余定理

约公元 100 年，孙子问：被 3、5、7 除分别余 2、3、2 的数是多少？答案 23（全部解为 \(23+105k\)）。

定理 31.27（中国剩余定理） 设 \(n=n_1n_2\cdots n_k\)，\(n_i\) 两两互素。则 \(a\leftrightarrow(a\bmod n_1,\dots,a\bmod n_k)\) 是 \(\mathbb Z_n\) 与 \(\mathbb Z_{n_1}\times\cdots\times\mathbb Z_{n_k}\) 之间的双射，且加、减、乘可在各分量上独立进行。

反向构造：令 \(m_i=n/n_i\)，\(c_i=m_i(m_i^{-1}\bmod n_i)\)，则

例：\(a\equiv2\pmod5\)，\(a\equiv3\pmod{13}\)。\(13^{-1}\equiv2\pmod5\)，\(5^{-1}\equiv8\pmod{13}\)，\(c_1=26\)，\(c_2=40\)，\(a\equiv52+120\equiv42\pmod{65}\)。

两种用途：结构上，\(\mathbb Z_n\) 等同于若干小模数的乘积；算法上，大模数下的运算可以拆到小模数下并行做，最后再拼回（RSA 解密的常用加速，以及多模数 FFT/哈希）。

31.6 元素的幂与反复平方

欧拉定理（定理 31.30）：\(a\in\mathbb Z_n^*\) 时 \(a^{\phi(n)}\equiv1\pmod n\)。费马小定理（定理 31.31）：\(p\) 为素数时 \(a^{p-1}\equiv1\pmod p\)。二者都是"元素的阶整除群的大小"的直接推论。

若 \(\mathrm{ord}_n(g)=\phi(n)\)，称 \(g\) 为模 \(n\) 的原根（primitive root），此时 \(\mathbb Z_n^*\) 是循环群，每个元素都是 \(g\) 的幂，指数称为离散对数。例：3 是模 7 的原根（\(3^i\bmod7\)：1,3,2,6,4,5），2 不是（1,2,4 循环）。\(\mathbb Z_n^*\) 有原根当且仅当 \(n=2,4,p^e,2p^e\)（\(p\) 为奇素数）。

定理 31.34 \(p\) 为奇素数时，\(x^2\equiv1\pmod{p^e}\) 只有 \(x\equiv\pm1\) 两个解。由此（推论 31.35）：若存在模 \(n\) 的"1 的非平凡平方根"，\(n\) 必为合数。例：\(6^2=36\equiv1\pmod{35}\)，所以 35 是合数。这是 Miller–Rabin 的理论基础。

反复平方法（repeated squaring）计算 \(a^b\bmod n\)：从高到低扫描 \(b\) 的二进制位，每位先平方，位为 1 再乘 \(a\)。

def modular_exponentiation(a, b, n):
    d = 1
    for bit in bin(b)[2:]:          # 从最高位到最低位
        d = d * d % n
        if bit == '1':
            d = d * a % n
    return d

循环不变式：处理完 \(b\) 的前缀 \(c\) 后，\(d=a^c\bmod n\)。原书图 31.4：\(a=7\)，\(b=560=(1000110000)_2\)，\(n=561\)，\(c\) 依次为 1,2,4,8,17,35,70,140,280,560，\(d\) 依次为 7,49,157,526,160,241,298,166,67,1。\(\beta\) 位输入需 \(O(\beta)\) 次模乘、\(O(\beta^3)\) 次位操作。同样的"反复平方"思想可用于矩阵幂：思考题 31-3 用 \(\begin{pmatrix}0&1\\1&1\end{pmatrix}^n\) 在 \(O(\lg n)\) 次矩阵乘法内求斐波那契数；量化中马尔可夫链的 \(n\) 步转移矩阵、AR 模型的多步预测都可以这样算。

31.7 RSA 公钥密码系统

公钥密码系统让每个参与者有一对密钥：公钥 \(P\) 公开，私钥 \(S\) 保密，二者互为逆变换。它有两种用法：

• 加密：Bob 用 Alice 的公钥算密文 \(C=P_A(M)\)，只有 Alice 能用 \(S_A\) 解密；
• 数字签名：Alice 用私钥算签名 \(\sigma=S_A(M')\)，任何人用公钥验证 \(P_A(\sigma)=M'\)。签名同时认证了签名者身份和消息内容，消息改动一位签名即失效。

RSA 密钥生成：随机选两个大素数 \(p,q\)（如各 1024 位），\(n=pq\)；选与 \(\phi(n)=(p-1)(q-1)\) 互素的小奇数 \(e\)；用扩展欧几里得求 \(d=e^{-1}\bmod\phi(n)\)。公钥 \((e,n)\)，私钥 \((d,n)\)。变换为

正确性（定理 31.36）：\(ed=1+k(p-1)(q-1)\)。若 \(p\nmid M\)，由费马小定理 \(M^{ed}=M\cdot(M^{p-1})^{k(q-1)}\equiv M\pmod p\)；\(p\mid M\) 时显然。模 \(q\) 同理，再由中国剩余定理得 \(M^{ed}\equiv M\pmod n\)。

安全性依赖大整数分解的困难：能分解 \(n\) 就能算出 \(d\)。实用中 RSA 很慢，所以采用混合模式：用 RSA 加密一个对称密码的短密钥，再用对称密码加密长消息；签名时先用抗碰撞哈希函数 \(h\) 算出消息的短"指纹"，只对 \(h(M)\) 签名。公钥的归属由受信任机构签发的证书（certificate）保证。

对量化工程师，这一节的直接价值在于理解交易所/券商 API 的认证：RSA/ECDSA 签名、TLS 握手、证书链都是这一套；"签名 = 私钥作用于消息哈希"这句话能帮你排查大部分验签失败（消息拼接顺序、编码、时间戳不一致导致哈希不同）。

31.8 素性测试

素数的密度。素数定理（定理 31.37）：\(\pi(n)\sim n/\ln n\)。\(n=10^9\) 时 \(\pi(n)=50{,}847{,}534\)，\(n/\ln n\approx48{,}254{,}942\)。所以随机取一个数是素数的概率约 \(1/\ln n\)，找一个 1024 位素数平均要试约 \(\ln2^{1024}\approx710\) 个随机数（只试奇数减半）。

费马测试（伪素数测试）：若 \(2^{n-1}\not\equiv1\pmod n\)，\(n\) 必为合数；否则"猜"它是素数。只会犯一种错：把以 2 为基的伪素数（如 341、561、645、1105）判为素数。10,000 以内只错 22 个；随机 512 位数上出错概率小于 \(10^{-20}\)。但换几个基也无法根除错误：Carmichael 数（561、1105、1729……）对所有 \(a\in\mathbb Z_n^*\) 都满足 \(a^{n-1}\equiv1\)。

Miller–Rabin 测试做了两点改进：随机选多个基 \(a\)；并且在计算 \(a^{n-1}\) 时检查 1 的非平凡平方根。写 \(n-1=2^tu\)（\(u\) 奇），先算 \(x_0=a^u\)，再连续平方 \(t\) 次。若某次平方得到 1 而平方前不是 \(\pm1\)，就找到了非平凡平方根，\(n\) 必为合数；若最后 \(a^{n-1}\ne1\)，费马测试失败，也是合数。

def witness(a, n):
    t, u = 0, n - 1
    while u % 2 == 0:
        t += 1; u //= 2
    x = modular_exponentiation(a, u, n)
    for _ in range(t):
        x_new = x * x % n
        if x_new == 1 and x != 1 and x != n - 1:
            return True                 # 找到 1 的非平凡平方根
        x = x_new
    return x != 1                       # a^{n-1} != 1：费马测试失败

例：\(n=561\)，\(n-1=2^4\cdot35\)，取 \(a=7\)。序列 \(x_0..x_4=241,298,166,67,1\)：最后一次从 67 平方得 1，67 是非平凡平方根，判合数。

定理 31.38 若 \(n\) 为奇合数，则 \(n\) 为合数的证据（witness）至少 \((n-1)/2\) 个。证明的核心是把所有"非证据"装进 \(\mathbb Z_n^*\) 的一个真子群，由拉格朗日定理其大小不超过一半。非 Carmichael 数时取 \(B=\{b:b^{n-1}\equiv1\}\)；Carmichael 数时先证它不是素数幂、写成两个互素奇数之积 \(n_1n_2\)，再用中国剩余定理构造一个不在 \(B=\{x:x^{2^ju}\equiv\pm1\}\) 中的元素。

定理 31.39 Miller–Rabin\((n,s)\) 的出错概率至多 \(2^{-s}\)，且出错只可能是"把合数判为素数"。

贝叶斯修正。\(2^{-s}\) 是"给定 \(n\) 是合数，通过检验的概率"，不是"通过检验后 \(n\) 是合数的概率"。随机取 \(\beta\) 位的 \(n\)，先验 \(\Pr\{\text{素数}\}\approx1/\ln n\)，于是

推导拆解：这个公式就是贝叶斯公式。记先验 \(\pi=\Pr\{\text{素数}\}\approx1/\ln n\)。素数一定通过检验（\(\Pr\{\text{通过}\mid\text{素数}\}=1\)），合数以至多 \(2^{-s}\) 的概率通过。于是

\[\Pr\{\text{素数}\mid\text{通过}\}=\frac{\pi\cdot1}{\pi\cdot1+(1-\pi)2^{-s}}=\frac{1}{1+2^{-s}\frac{1-\pi}{\pi}}.\]

把 \(\pi=1/\ln n\) 代入，\(\frac{1-\pi}{\pi}=\ln n-1\)，就得到原文的式子。它说明：先验越低（\(\ln n\) 越大，随机数越可能是合数），需要的检验轮数越多，后验才能过半。

这和 CFA 里"罕见病检测"的例题结构完全相同：检测的准确率很高，但如果患病率很低，阳性结果中仍可能大多是误报。上一段末尾的提醒（实际误判率远小于 \(2^{-s}\)）意味着这里给出的是后验的保守下界。

章末注记提到，2002 年 Agrawal–Kayal–Saxena 给出了确定性多项式时间素性测试（AKS），但实践中随机化测试仍更快。

31.9 整数分解：Pollard rho

判素容易，分解难——RSA 的安全性就建立在这个差距上。Pollard rho 启发式能快速找到大数的小素因子：迭代 \(x_{i+1}=(x_i^2-1)\bmod n\)，保存下标为 2 的幂的那些 \(x\) 作为 \(y\)，每步计算 \(\gcd(y-x_i,n)\)，得到非平凡值就输出因子。

为什么有效：设 \(p\) 是 \(n\) 的一个小因子，序列模 \(p\) 的值 \(x_i\bmod p\) 满足同一个递推，且只有 \(p\) 种取值。把它近似看作随机函数，由生日悖论，期望 \(\Theta(\sqrt p)\) 步就会出现重复，进入一个圈——画出来像希腊字母 ρ，故名。一旦 \(x_i\equiv x_j\pmod p\)，\(p\) 就整除 \(x_i-x_j\)，gcd 把它揭示出来。所以找到因子 \(p\) 期望只需 \(\Theta(\sqrt p)\) 次运算，而且只用常数个存储单元。原书图 31.7：\(n=1387=19\cdot73\)，\(x_1=2\)，序列 2,3,8,63,1194,1186,177,…，在 \(x_7=177\) 时 \(\gcd(63-177,1387)=19\)。完全分解 \(\beta\) 位合数期望约 \(n^{1/4}=2^{\beta/4}\) 次运算，仍是指数级。保存 \(x_1,x_2,x_4,\dots\) 的技巧本身是一种判圈算法（习题 31.9-2；Floyd/Brent 算法），可以用来检测任何确定性迭代（伪随机数生成器、状态机）是否进入循环。

第二部分 字符串匹配（原书第 32 章）

32.1 问题与算法概览

文本 \(T[1..n]\)、模式 \(P[1..m]\) 都是有限字母表 \(\Sigma\) 上的字符串。若 \(T[s+1..s+m]=P[1..m]\)，称 \(P\) 以偏移 \(s\) 出现，\(s\) 为有效偏移（valid shift）。字符串匹配问题要找出全部有效偏移。原书图 32.2 总结四种算法：

记号：\(w\sqsubset x\) 表示 \(w\) 是 \(x\) 的前缀，\(w\sqsupset x\) 表示后缀；\(P_k=P[1..k]\)。重叠后缀引理（引理 32.1）：\(x,y\) 都是 \(z\) 的后缀时，短的必是长的后缀。

朴素算法把模式当模板在文本上逐位滑动，最坏 \(\Theta((n-m+1)m)\)（如 \(T=a^n\)、\(P=a^m\)）。它低效的根源是丢弃了已比较过的文本信息。不过在随机文本上它期望只做约 \(2(n-m+1)\) 次比较（习题 32.1-3），并不差。

32.2 Rabin–Karp：滚动哈希

把长度为 \(m\) 的字符串看作 \(d\) 进制数（\(d=|\Sigma|\)）。模式的值 \(p\) 用 Horner 法则 \(\Theta(m)\) 算出；文本窗口的值可以滚动更新：

金融直觉：滚动哈希与滚动窗口求和是同一个技巧。算 20 日移动和时，不必每天把 20 个数重新加一遍，只要"加上新的一天、减去滑出窗口的那一天"，每步 \(O(1)\)。Rabin–Karp 对"窗口的数值"做同样的事：减去最高位（滑出的字符），整体左移一位（乘 \(d\)），加上新的最低位（新字符）。

取模 \(q\) 的作用类似于给每个窗口算一个"校验码"：校验码不同，内容一定不同；校验码相同，内容多半相同，但要再逐字核对一次（伪命中）。这和对账时先比汇总金额、金额一致再逐笔核对是一个思路。

数值太大怎么办？对一个素数 \(q\) 取模：

def rabin_karp(T, P, d=256, q=1_000_000_007):
    n, m = len(T), len(P); h = pow(d, m - 1, q)
    p = t = 0
    for i in range(m):
        p = (d * p + ord(P[i])) % q; t = (d * t + ord(T[i])) % q
    out = []
    for s in range(n - m + 1):
        if p == t and T[s:s + m] == P:            # 哈希相等再核对
            out.append(s)
        if s < n - m:
            t = (d * (t - ord(T[s]) * h) + ord(T[s + m])) % q
    return out

若把模 \(q\) 看作随机映射，伪命中期望 \(O(n/q)\) 个，期望匹配时间 \(O(n)+O(m(v+n/q))\)（\(v\) 为有效偏移数），\(v=O(1)\)、\(q\ge m\) 时为 \(O(n+m)\)。Rabin–Karp 的优势是容易推广：多个模式（把模式的哈希放进哈希表，习题 32.2-2）、二维模式（习题 32.2-3）。习题 32.2-4 的多项式指纹更是一个通用思想：比较两个长文件是否相同，只需比较 \(\sum a_ix^i\bmod q\) 在随机点 \(x\) 处的值，碰撞概率不超过 \(n/q\)——数据同步中的校验、rsync 的滚动校验和都是这个原理。

32.3 有限自动机

有限自动机（finite automaton）\(M=(Q,q_0,A,\Sigma,\delta)\) 从初始状态出发，每读一个字符按转移函数 \(\delta\) 换一个状态，处于接受状态集 \(A\) 时"接受"。

为模式 \(P\) 构造的匹配自动机：状态 \(q\in\{0,\dots,m\}\) 表示"\(P\) 的、同时是已读文本后缀的最长前缀的长度"。用后缀函数 \(\sigma(x)=\max\{k:P_k\sqsupset x\}\) 定义转移

缺点是转移表有 \((m+1)|\Sigma|\) 项，字母表大时预处理昂贵（朴素构造 \(O(m^3|\Sigma|)\)，利用前缀函数可降到 \(O(m|\Sigma|)\)）。

32.4 KMP 算法

KMP 绕开转移表，只用一个长度 \(m\) 的辅助数组——前缀函数（prefix function）：

白话解释：用表中的 \(P=ababaca\) 举例。假设已经匹配了前 5 个字符 \(ababa\)，第 6 个字符失配。朴素算法会把模式右移一位、从头重新比较。KMP 注意到：刚匹配上的 \(ababa\) 的结尾 \(aba\) 恰好也是模式的开头 \(aba\)，所以可以直接把模式右移 2 位，让开头的 \(aba\) 对准文本中刚读过的 \(aba\)，并且已知这 3 个字符不用再比，从第 4 个字符接着比。\(\pi[5]=3\) 记录的就是这个"3"。

关键在于，文本中刚读过的内容恰好就是模式的前 \(q\) 个字符，所以"失配后能保留多少"只取决于模式本身，可以事先算好，匹配时文本指针永远不回退。这正是它能逐 tick 流式运行的原因。

\(P=ababaca\) 的前缀函数（原书图 32.11）：

def compute_prefix_function(P):          # 0 下标版：pi[q] 对应原书 π[q+1]
    m = len(P); pi = [0] * m; k = 0
    for q in range(1, m):
        while k > 0 and P[k] != P[q]:
            k = pi[k - 1]
        if P[k] == P[q]:
            k += 1
        pi[q] = k
    return pi

def kmp_matcher(T, P):
    pi = compute_prefix_function(P); m = len(P); q = 0; out = []
    for i, c in enumerate(T):
        while q > 0 and P[q] != c:
            q = pi[q - 1]                # 失配：回退
        if P[q] == c:
            q += 1                       # 匹配：前进
        if q == m:
            out.append(i - m + 1)        # 找到一次出现
            q = pi[q - 1]                # 继续找下一次（允许重叠）
    return out

两个过程结构相同：一个把 \(T\) 与 \(P\) 匹配，一个把 \(P\) 与自身匹配。

运行时间（聚合分析）：以 COMPUTE-PREFIX-FUNCTION 为例，\(k\) 每轮至多加 1，总增量不超过 \(m-1\)；while 循环每执行一次都使 \(k\) 严格减小，且 \(k\) 永不为负。所以 while 的总次数不超过总增量，预处理 \(\Theta(m)\)；同理匹配 \(\Theta(n)\)。这是第 17 章聚合分析的典型应用。

正确性：引理 32.5 说明反复迭代 \(\pi\)（\(\pi[q],\pi[\pi[q]],\dots\)）恰好枚举出所有"是 \(P_q\) 真后缀的前缀"；推论 32.7 据此给出 \(\pi[q]\) 的递推。KMP 的 while 循环就是按从长到短的顺序尝试这些候选，因此它在每一步得到的 \(q\) 与有限自动机的状态相同，只是把 \(\delta\) 按需"现算"。

第三部分 计算几何（原书第 33 章）

33.1 叉积：几何算法的基本原语

叉积（cross product）

它可以只用加、减、乘和比较回答三个基本问题，完全避免除法和三角函数（后两者代价高，而且在线段接近平行时对舍入误差极其敏感）：

1. 转向：走 \(p_0\to p_1\to p_2\)，看 \((p_2-p_0)\times(p_1-p_0)\)：负则在 \(p_1\) 处左转，正则右转，零则共线。（本章代码用等价的 \((p_1-p_0)\times(p_2-p_0)\)，正为左转。）
2. 线段相交：两线段相交当且仅当"互相跨越对方所在直线"，或"某个端点落在另一线段上"（边界情况）。四次叉积给出四个端点相对另一条线段的方向 \(d_1,\dots,d_4\)：\(d_1,d_2\) 异号且 \(d_3,d_4\) 异号即互相跨越；某个 \(d_k=0\) 时，再用坐标包围盒检查该端点是否在线段上。\(O(1)\) 时间。
3. 极角排序：比较两点相对某原点的极角，只要看叉积的符号，不必算角度。

习题中的几个小工具很实用：射线法判断点是否在多边形内（33.1-7），鞋带公式 \(\frac12\left|\sum(x_iy_{i+1}-x_{i+1}y_i)\right|\) 求多边形面积（33.1-8）。

33.2 扫描线：判断是否有线段相交

\(n\) 条线段中是否存在相交的一对？两两检查要 \(\Theta(n^2)\)。扫描（sweeping）技术用一条假想的竖直扫描线从左向右扫过平面，把 \(x\) 坐标当作时间。它维护两样东西：

• 扫描线状态：与扫描线相交的线段按上下顺序排成的全预序，用红黑树存储，比较用叉积；
• 事件点调度：所有端点按 \(x\) 坐标排序（同 \(x\) 时左端点先于右端点，同类按 \(y\)）。

遇到左端点就插入线段，并检查它与上下相邻线段是否相交；遇到右端点，先检查它的上下邻居是否相交（删除后二者会变成相邻），再删除。正确性的关键是：最左的交点处，相交的两条线段在某个事件点上一定会成为相邻。总时间 \(O(n\lg n)\)。原书同时说明：把"返回"改成"打印并继续"并不能输出所有交点（习题 33.2-3），要输出全部 \(k\) 个交点需把交点也作为动态事件（Bentley–Ottmann，\(O((n+k)\lg n)\)，习题 33.2-7）。

一维的扫描线在量化系统里很常见：多源事件按时间戳合并、统计订单生命周期的重叠数、检测时间窗口冲突。按"同一时刻先开始后结束（或反之）"的规则排序事件，正是本节端点排序规则的一维版本。

33.3 凸包

点集 \(Q\) 的凸包（convex hull）\(CH(Q)\) 是包含 \(Q\) 中所有点的最小凸多边形——把点想成木板上的钉子，凸包就是箍住所有钉子的橡皮筋。

Graham 扫描（\(O(n\lg n)\)）：

1. 取 \(y\) 最小（并列取最左）的点 \(p_0\)，它一定是凸包顶点；
2. 其余点按相对 \(p_0\) 的极角逆时针排序（叉积比较），同极角只保留最远者；
3. 用栈扫描：依次压入各点，压入前若栈顶两点与新点构成非左转（右转或共线），就弹出栈顶，直到左转为止。

结束时栈中自底向顶就是逆时针顺序的凸包顶点。正确性（定理 33.1）的循环不变式：处理第 \(i\) 个点前，栈恰为前 \(i-1\) 个点的凸包。被弹出的点落在某个三角形 \(p_0p_rp_i\) 内部或边上，不可能是凸包顶点。运行时间：排序 \(O(n\lg n)\)；每个点入栈一次、出栈至多一次，栈操作总共 \(O(n)\)——又一次聚合分析。

Jarvis 步进（gift wrapping，\(O(nh)\)，\(h\) 为凸包顶点数）：从最低点出发，每次选相对当前点极角最小的点作为下一个顶点，像包礼物一样绕一圈。\(h\) 很小时（\(h=o(\lg n)\)）比 Graham 快。

其他方法：增量法、分治法 \(O(n\lg n)\)；Kirkpatrick–Seidel 的剪枝搜索 \(O(n\lg h)\) 渐近最优。下界：在比较模型下按顺序输出凸包需 \(\Omega(n\lg n)\)，因为可以把排序归约到"求抛物线 \(y=x^2\) 上的点的凸包"（习题 33.3-2）。凸包常作为其他问题的第一步，例如平面最远点对一定是凸包顶点（旋转卡壳）。

33.4 最近点对

在 \(n\) 个点中找欧氏距离最近的两点。分治：按 \(x\) 坐标把点集平分为左右两半，递归求出两侧最近距离，取 \(\delta=\min(\delta_L,\delta_R)\)。跨越分界线的更近点对只能落在以分界线为中心、宽 \(2\delta\) 的竖直带内；带内按 \(y\) 排序后，每个点只需与后面 7 个点比较。原因是 \(\delta\times2\delta\) 的矩形中，左右两半各至多容纳 4 个两两距离不小于 \(\delta\) 的点。关键实现技巧是预排序：只在开始时按 \(x\)、\(y\) 各排序一次，之后每层用线性时间把有序数组拆成有序子数组（归并的逆操作），否则会变成 \(O(n\lg^2n)\)。总时间 \(T(n)=2T(n/2)+O(n)=O(n\lg n)\)。

高维的近邻问题（相似股票、相似行情片段检索）用 KD 树、近似最近邻（ANN）等方法，平面算法不能直接套用，但"分治 + 只检查边界附近"的思想相通。

33.5 思考题中的两个重要概念

凸层（思考题 33-1）：剥掉凸包，剩余点再求凸包，如此反复。

极大层（思考题 33-2）：若 \(x\ge x'\) 且 \(y\ge y'\)，称 \((x,y)\) 支配（dominates）\((x',y')\)。不被任何点支配的点称为极大点（maximal），它们构成第 1 极大层；删去后剩下点的极大点构成第 2 层，依此类推。极大点就是帕累托前沿（Pareto frontier）。原书给出 \(O(n\lg n)\) 算法：按 \(x\) 从大到小扫描，维护每层当前最左点的 \(y\) 值（它们严格递减），新点用二分查找确定归入哪一层。

31.10 量化实战

31.10.1 数论工具箱的验证

先把第一部分的算法写成可运行代码，复现原书的数字，并完成 RSA 玩具例（原书习题 31.7-1：\(p=11,q=29,e=3\)）。

import random

def extended_euclid(a, b):
    if b == 0:
        return (a, 1, 0)
    d, x, y = extended_euclid(b, a % b)
    return (d, y, x - (a // b) * y)

def mod_inverse(a, n):
    d, x, _ = extended_euclid(a, n)
    if d != 1:
        raise ValueError("not invertible")
    return x % n

def modular_exponentiation(a, b, n):
    d = 1
    for bit in bin(b)[2:]:
        d = d * d % n
        if bit == '1':
            d = d * a % n
    return d

def witness(a, n):
    t, u = 0, n - 1
    while u % 2 == 0:
        t += 1; u //= 2
    x = modular_exponentiation(a, u, n)
    for _ in range(t):
        x_new = x * x % n
        if x_new == 1 and x != 1 and x != n - 1:
            return True
        x = x_new
    return x != 1

def miller_rabin(n, s=50, rng=random.Random(0)):
    if n in (2, 3): return True
    if n < 2 or n % 2 == 0: return False
    return not any(witness(rng.randint(1, n - 1), n) for _ in range(s))

def crt(a_list, n_list):
    n = 1
    for ni in n_list: n *= ni
    x = 0
    for ai, ni in zip(a_list, n_list):
        mi = n // ni
        x += ai * mi * mod_inverse(mi % ni, ni)
    return x % n

print("EXTENDED-EUCLID(99,78) =", extended_euclid(99, 78))
print("5^{-1} mod 11 =", mod_inverse(5, 11), "; 7^{-1} mod 15 =", mod_inverse(7, 15))
print("CRT([2,3],[5,13]) =", crt([2, 3], [5, 13]), "; 孙子问题 CRT([2,3,2],[3,5,7]) =", crt([2, 3, 2], [3, 5, 7]))
print("7^560 mod 561 =", modular_exponentiation(7, 560, 561), "（561=3·11·17 是 Carmichael 数）")
print("费马测试 2^560 mod 561 =", modular_exponentiation(2, 560, 561), "-> 被骗；Miller-Rabin:", miller_rabin(561))
print("witness(7,561) =", witness(7, 561))
carm = [561, 1105, 1729, 2465, 2821, 6601, 8911]
print("Carmichael 数被 Miller-Rabin 判为素数的个数:", sum(miller_rabin(c) for c in carm))
p = 2**61 - 1
print("2^61-1 是素数?", miller_rabin(p), "; 2^61+1 是素数?", miller_rabin(2**61 + 1))

p_, q_, e = 11, 29, 3
n = p_ * q_; phi = (p_ - 1) * (q_ - 1)
d = mod_inverse(e, phi)
M = 100
C = modular_exponentiation(M, e, n)
print(f"n={n}, φ(n)={phi}, d={d}, 密文 C={C}, 解密得 {modular_exponentiation(C, d, n)}")
sig = modular_exponentiation(42, d, n)        # 用私钥对"消息摘要"42 签名
print(f"签名 σ={sig}, 用公钥验证 σ^e mod n = {modular_exponentiation(sig, e, n)}")

输出：

EXTENDED-EUCLID(99,78) = (3, -11, 14)
5^{-1} mod 11 = 9 ; 7^{-1} mod 15 = 13
CRT([2,3],[5,13]) = 42 ; 孙子问题 CRT([2,3,2],[3,5,7]) = 23
7^560 mod 561 = 1 （561=3·11·17 是 Carmichael 数）
费马测试 2^560 mod 561 = 1 -> 被骗；Miller-Rabin: False
witness(7,561) = True
Carmichael 数被 Miller-Rabin 判为素数的个数: 0
2^61-1 是素数? True ; 2^61+1 是素数? False
n=319, φ(n)=280, d=187, 密文 C=254, 解密得 100
签名 σ=180, 用公钥验证 σ^e mod n = 42

这只是教学用的玩具。实盘系统里的签名、加密必须用经过审计的密码库（如 Python 的 cryptography），不要自己实现 RSA：真实实现还需要填充方案（OAEP、PSS）、常数时间运算以防侧信道等，本书不涉及。Python 内置的 pow(a, b, n) 就是反复平方法，pow(a, -1, n) 直接给出模逆元。

31.10.2 随机数生成器的周期与"检验的后验"

乘法同余生成器的周期。经典的乘法同余生成器（Lehmer 生成器）\(x_{i+1}=ax_i\bmod m\)，\(m\) 为素数时，状态在 \(\mathbb Z_m^*\) 中运动，周期恰好是 \(\mathrm{ord}_m(a)\)。由拉格朗日定理 \(\mathrm{ord}_m(a)\mid m-1\)；要验证 \(a\) 是原根（满周期 \(m-1\)），只需对 \(m-1\) 的每个素因子 \(q\) 检查 \(a^{(m-1)/q}\not\equiv1\)。

贝叶斯修正与因子挖掘。31.8 节的后验公式对任何"先验很低、检验有一定误报率"的筛选都成立。在大量候选因子中，真正有效的比例（先验）很低；即使单次检验的显著性水平是 5%，通过检验的因子也可能大多是假阳性。

import math

def factorize(n):
    f, d = {}, 2
    while d * d <= n:
        while n % d == 0:
            f[d] = f.get(d, 0) + 1; n //= d
        d += 1
    if n > 1: f[n] = f.get(n, 0) + 1
    return f

def order(a, m):
    """m 为素数时，ord_m(a) 整除 m-1：从 m-1 出发逐个去掉素因子"""
    t = m - 1
    for q in factorize(m - 1):
        while t % q == 0 and pow(a, t // q, m) == 1:
            t //= q
    return t

m = 2**31 - 1                          # Mersenne 素数
print("m-1 的分解:", factorize(m - 1))
for a in (16807, 48271, 65539 % m, 2):
    print(f"a={a:6d}: 周期 ord_m(a) = {order(a, m):,}" + ("（原根，满周期）" if order(a, m) == m - 1 else ""))

def posterior(prior, power, alpha):
    return prior * power / (prior * power + (1 - prior) * alpha)

beta = 1024
prior = 1 / (beta * math.log(2))       # 随机 1024 位整数为素数的先验约 1/ln n
for s in (1, 5, 9, 20):
    print(f"Miller-Rabin s={s:2d}: P(素数 | 通过) >= {posterior(prior, 1.0, 2.0**-s):.4f}")
print("因子挖掘：先验 1% 为真，检验功效 80%")
for alpha in (0.05, 0.01, 0.001):
    print(f"  显著性 {alpha:<6}: P(真因子 | 通过) = {posterior(0.01, 0.8, alpha):.3f}")

输出：

m-1 的分解: {2: 1, 3: 2, 7: 1, 11: 1, 31: 1, 151: 1, 331: 1}
a= 16807: 周期 ord_m(a) = 2,147,483,646（原根，满周期）
a= 48271: 周期 ord_m(a) = 2,147,483,646（原根，满周期）
a= 65539: 周期 ord_m(a) = 1,073,741,823
a=     2: 周期 ord_m(a) = 31
Miller-Rabin s= 1: P(素数 | 通过) >= 0.0028
Miller-Rabin s= 5: P(素数 | 通过) >= 0.0432
Miller-Rabin s= 9: P(素数 | 通过) >= 0.4194
Miller-Rabin s=20: P(素数 | 通过) >= 0.9993
因子挖掘：先验 1% 为真，检验功效 80%
  显著性 0.05  : P(真因子 | 通过) = 0.139
  显著性 0.01  : P(真因子 | 通过) = 0.447
  显著性 0.001 : P(真因子 | 通过) = 0.890

解读：

• 16807（MINSTD）和 48271 都是模 \(2^{31}-1\) 的原根，周期约 21 亿；65539 只有一半周期；2 的周期只有 31，完全不能用。约 21 亿的周期对现代蒙特卡洛来说太短（一次大规模路径模拟就可能用完），而且低维格点结构明显。实务中应使用 PCG64（numpy 默认）、Mersenne Twister 或 Philox 等生成器，并通过 SeedSequence.spawn 给并行任务分配独立的流。
• Miller–Rabin 的数字和原书一致：\(\lg(\ln n-1)\approx9.5\) 轮之前后验不到 1/2。
• 因子挖掘的数字是同一个公式：若 100 个候选里只有 1 个真因子，用 5% 显著性筛选，通过者中只有约 14% 是真的；把门槛提到 0.1%，才接近九成。这就是多重检验下要提高 \(t\) 值门槛（例如 \(t>3\)）、或做 FDR 控制的理由，详见第 03 册第 10b 章（多重检验与策略过拟合）。

31.10.3 行情符号化与流式形态识别

把价格序列离散成符号串（涨 U、跌 D、平 F，或更细的 SAX 表示），"K 线形态识别"就成了字符串匹配。KMP 的一个工程优点是可以流式运行：只保存当前状态 \(q\)，每来一个符号摊还 \(O(1)\) 时间，不需要回看历史数据。

import numpy as np, time
# compute_prefix_function, kmp_matcher, rabin_karp 定义同上

def naive(T, P):
    m = len(P)
    return [s for s in range(len(T) - m + 1) if T[s:s + m] == P]

print("π(ababaca) =", compute_prefix_function("ababaca"))
print("KMP 在 abababacaba 中找 ababaca:", kmp_matcher("abababacaba", "ababaca"))

rng = np.random.default_rng(3)
r = rng.standard_normal(1_000_000) * 0.01
sym = np.where(r > 0.005, 'U', np.where(r < -0.005, 'D', 'F'))
T = ''.join(sym)
P = "DDDU"                                         # 连跌三天后反弹
for f in (naive, rabin_karp, kmp_matcher):
    t0 = time.perf_counter(); occ = f(T, P); dt = time.perf_counter() - t0
    print(f"{f.__name__:12s} 找到 {len(occ)} 次, 用时 {dt:.2f}s")

occ = np.array(kmp_matcher(T, P))                  # 形态之后的下一日收益是否显著？
nxt = occ + len(P)
nxt = nxt[nxt < len(r)]
fwd = r[nxt]
tstat = fwd.mean() / fwd.std(ddof=1) * np.sqrt(len(fwd))
print(f"形态后下一日平均收益 {fwd.mean()*1e4:.2f} bp, t = {tstat:.2f}, 样本 {len(fwd)}")

class StreamingKMP:                                # 流式版本：每来一个符号调用一次 feed
    def __init__(self, P):
        self.P, self.pi, self.q = P, compute_prefix_function(P), 0
    def feed(self, c):
        P, pi, q = self.P, self.pi, self.q
        while q > 0 and P[q] != c:
            q = pi[q - 1]
        if P[q] == c:
            q += 1
        hit = (q == len(P))
        if hit:
            q = pi[q - 1]
        self.q = q
        return hit

det = StreamingKMP(P)
hits = [i - len(P) + 1 for i, c in enumerate(T) if det.feed(c)]
print("流式检测与批量 KMP 结果一致:", hits == kmp_matcher(T, P))

输出：

π(ababaca) = [0, 0, 1, 2, 3, 0, 1]
KMP 在 abababacaba 中找 ababaca: [2]
naive        找到 8972 次, 用时 0.03s
rabin_karp   找到 8972 次, 用时 0.08s
kmp_matcher  找到 8972 次, 用时 0.03s
形态后下一日平均收益 1.04 bp, t = 0.99, 样本 8972
流式检测与批量 KMP 结果一致: True

三点说明：

• 在 Python 里，朴素算法的每次比较由 C 实现的切片比较完成，模式又很短，所以它并不慢；Rabin–Karp 的逐字符 Python 运算反而最慢。渐近复杂度要在同一实现层面比较才有意义；生产中直接用 str.find、正则引擎或 C++ 实现。KMP 真正的价值是流式：在逐 tick 的行情回调里，它只需一个整数状态。
• 本例中模拟收益是独立同分布的，形态之后的收益不显著（\(t=0.99\)），这是对的。真实数据上找到"显著"形态时，要记住上一小节的后验公式：尝试过的形态越多，单个形态的 \(t\) 值门槛就应越高。算法高效不代表信号有效。
• 多个模式同时匹配（几百个关键词扫描公告、新闻）用 Aho–Corasick 自动机，它是 KMP 的多模式推广；近重复新闻检测用滚动哈希做 shingling，再配合 MinHash。

31.10.4 凸包、有效前沿与套利修正

凸包近似有效前沿。在 (波动率 \(\sigma\), 期望收益 \(\mu\)) 平面上，两个组合 \(A\)、\(B\) 按 \(\lambda,1-\lambda\) 混合后，收益是线性插值，而波动率满足 \(\sigma_{\lambda A+(1-\lambda)B}\le\lambda\sigma_A+(1-\lambda)\sigma_B\)（标准差的三角不等式）。所以混合组合位于连接 \(A\)、\(B\) 的线段上方或左方：可行集的左上边界是凹的，凸包上左链上任意一点都被某个真实可达的组合弱支配。因此，对一批候选组合（随机生成，或者历史上的实际持仓）求凸包，其上左链就是有效前沿的一个内侧近似。

金融直觉：这段话是 CFA 组合管理里"分散化效应"的几何表述。两资产组合的标准差 \(\sigma_p=\sqrt{\lambda^2\sigma_A^2+(1-\lambda)^2\sigma_B^2+2\lambda(1-\lambda)\rho\sigma_A\sigma_B}\)，只有相关系数 \(\rho=1\) 时才等于 \(\lambda\sigma_A+(1-\lambda)\sigma_B\)；\(\rho<1\) 时更小。所以在 \((\sigma,\mu)\) 图上，两个组合的混合会"向左鼓出"，有效前沿因此是一条向左上凸出的曲线。凸包的上左链用直线段连接已有的组合，相当于假设 \(\rho=1\)、不享受分散化，所以它总在真实前沿的右侧（内侧），是一个保守的近似。

期权凸性的直觉也一样：蝶式组合（买 \(K-h\)、卖 2 份 \(K\)、买 \(K+h\)）到期收益永远非负，所以今天价格不能为负，即 \(C(K-h)-2C(K)+C(K+h)\ge0\)，这正是"\(C\) 关于 \(K\) 是凸函数"的离散形式。取下凸包，就是用不违反这个条件、又尽量贴近原报价的折线替换报价曲线。

帕累托分层。比较策略时往往有多个目标（收益高、回撤小），它们只构成偏序，没有唯一的"最好"。极大层算法一次给出所有策略的帕累托层级，第 1 层就是不被任何策略支配的那些。

期权报价的凸性。同一到期日的看涨期权价格对执行价必须是凸函数，否则买 1 份 \(K-h\)、卖 2 份 \(K\)、买 1 份 \(K+h\) 的蝶式组合价格为负，构成套利。带噪声的报价常违反这一点；取报价点的下凸包（最大凸下界）就能得到一组最接近原报价、又满足凸性的价格。

import numpy as np
from functools import cmp_to_key
from bisect import bisect_left
from scipy.optimize import minimize
from scipy.stats import norm

def cross(o, a, b):
    """(a-o) × (b-o)：>0 表示 o→a→b 左转（逆时针），<0 右转，=0 共线"""
    return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])

def graham_scan(pts):
    pts = [tuple(p) for p in pts]
    p0 = min(pts, key=lambda p: (p[1], p[0]))               # y 最小，并列取最左
    rest = [p for p in pts if p != p0]
    def by_angle(a, b):                                      # 用叉积比较极角
        c = cross(p0, a, b)
        if c > 0: return -1
        if c < 0: return 1
        da = (a[0]-p0[0])**2 + (a[1]-p0[1])**2; db = (b[0]-p0[0])**2 + (b[1]-p0[1])**2
        return -1 if da > db else 1                          # 同极角：远的排前
    rest.sort(key=cmp_to_key(by_angle))
    uniq = []                                                # 同极角只保留最远者
    for p in rest:
        if uniq and cross(p0, uniq[-1], p) == 0:
            continue
        uniq.append(p)
    S = [p0, uniq[0], uniq[1]]
    for p in uniq[2:]:
        while len(S) >= 2 and cross(S[-2], S[-1], p) <= 0:  # 非左转就弹栈
            S.pop()
        S.append(p)
    return S                                                 # 逆时针顺序的凸包顶点

# ---- 1) 随机组合在 (σ, μ) 平面的凸包 vs 真实有效前沿 ----
rng = np.random.default_rng(11)
N = 6
mu = rng.uniform(0.02, 0.15, N)
A = rng.standard_normal((N, N)); vol = rng.uniform(0.10, 0.35, N)
S = A @ A.T + N * np.eye(N); d = 1 / np.sqrt(np.diag(S)); C = S * np.outer(d, d)   # 相关矩阵
Sigma = C * np.outer(vol, vol)
W = rng.dirichlet(np.full(N, 0.3), 20000)                    # 2 万个随机多头组合
pts = np.column_stack([np.sqrt(np.einsum('ij,jk,ik->i', W, Sigma, W)), W @ mu])
H = graham_scan(pts)
top = max(range(len(H)), key=lambda i: H[i][1]); left = min(range(len(H)), key=lambda i: H[i][0])
chain = [H[i % len(H)] for i in range(top, left + (len(H) if left < top else 0) + 1)]  # 逆时针从最高点走到最左点
print(f"凸包顶点 {len(H)} 个，其中上左链（有效部分）{len(chain)} 个")

def frontier_sigma(target):                                  # 多头约束下给定收益的最小波动
    cons = [{'type': 'eq', 'fun': lambda w: w.sum() - 1}, {'type': 'eq', 'fun': lambda w: w @ mu - target}]
    res = minimize(lambda w: w @ Sigma @ w, np.full(N, 1 / N), bounds=[(0, 1)] * N, constraints=cons, method='SLSQP')
    return np.sqrt(res.fun)
cs = np.array(chain)[::-1]                                   # 按 σ 升序
for target in np.linspace(cs[:, 1].min(), cs[:, 1].max(), 4):
    s_hull = np.interp(target, cs[:, 1], cs[:, 0])
    print(f"目标收益 {target:.3f}: 凸包链 σ = {s_hull:.4f}, 优化前沿 σ = {frontier_sigma(target):.4f}")

# ---- 2) 帕累托极大层（原书思考题 33-2）：策略按 (收益, -回撤) 分层 ----
def maximal_layers(points):
    """O(n lg n)：按 x 从大到小扫描，ys[i] 为第 i 层当前最左点的 y（严格递减）"""
    order = sorted(range(len(points)), key=lambda i: -points[i][0])
    ys, layer = [], [0] * len(points)
    for i in order:
        y = points[i][1]
        j = bisect_left([-v for v in ys], -y)                # 找第一个 ys[j] < y 的层
        if j == len(ys): ys.append(y)
        else: ys[j] = y
        layer[i] = j + 1
    return layer
strat = np.column_stack([rng.normal(0.08, 0.05, 200), -np.abs(rng.normal(0.15, 0.08, 200))])
lay = maximal_layers([tuple(p) for p in strat])
brute = [i for i in range(200) if not any((strat[j] >= strat[i]).all() and (strat[j] > strat[i]).any() for j in range(200))]
print(f"200 个策略：第 1 层（帕累托前沿）{lay.count(1)} 个，与暴力法一致: {sorted(brute) == [i for i in range(200) if lay[i] == 1]}；共 {max(lay)} 层")

# ---- 3) 期权报价的下凸包：消除蝶式套利 ----
K = np.arange(60, 142, 2.0); S0, r, T, sig = 100, 0.02, 0.25, 0.3
d1 = (np.log(S0 / K) + (r + sig**2 / 2) * T) / (sig * np.sqrt(T))
call = S0 * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d1 - sig * np.sqrt(T))
quote = np.maximum(call + rng.normal(0, 0.05, len(K)), 0)    # 带噪声的报价
def lower_hull(x, y):                                       # 单调链：只保留下凸包
    H = []
    for p in zip(x, y):
        while len(H) >= 2 and cross(H[-2], H[-1], p) <= 0:
            H.pop()
        H.append(p)
    return np.array(H)
LH = lower_hull(K, quote)
fixed = np.interp(K, LH[:, 0], LH[:, 1])
bfly = lambda c: np.sum(c[:-2] - 2 * c[1:-1] + c[2:] < -1e-12)   # 等间距蝶式价格为负即套利
print(f"蝶式套利个数：原始报价 {bfly(quote)}，下凸包修正后 {bfly(fixed)}；最大修正幅度 {np.max(quote - fixed):.3f}")

输出：

凸包顶点 26 个，其中上左链（有效部分）10 个
目标收益 0.093: 凸包链 σ = 0.0719, 优化前沿 σ = 0.0715
目标收益 0.109: 凸包链 σ = 0.0769, 优化前沿 σ = 0.0761
目标收益 0.125: 凸包链 σ = 0.0865, 优化前沿 σ = 0.0862
目标收益 0.141: 凸包链 σ = 0.1041, 优化前沿 σ = 0.1040
200 个策略：第 1 层（帕累托前沿）4 个，与暴力法一致: True；共 26 层
蝶式套利个数：原始报价 13，下凸包修正后 0；最大修正幅度 0.125

解读：

• 2 万个随机组合的凸包只有 26 个顶点，上左链 10 个；在每个目标收益上，凸包链的波动率都略高于优化得到的前沿（内侧近似），差距在 0.1 个百分点以内。资产数一多，随机抽样很难覆盖前沿附近的"角落"组合，差距会变大——凸包适合对已有的有限组合集合（历史持仓、候选策略组合）做筛选，求真正的前沿还是要解二次规划（第 04 册第 16a、16b 章）。
• 极大层算法的结果与 \(O(n^2)\) 暴力法一致。第 1 层只有 4 个策略，这正说明多目标下"最优"往往不唯一；后续可以在第 1 层内部再用夏普比率、容量等指标挑选。
• 下凸包消除了全部 13 处蝶式套利。它只保证凸性，完整的无套利修正还要求斜率在 \([-e^{-rT},0]\) 内（看涨价格关于执行价单调递减、且下降速度有界）以及跨到期日的日历价差约束，通常写成带约束的最小二乘问题求解。

本章小结

数论部分：gcd 是最小正线性组合，欧几里得算法递归 \(O(\lg b)\) 次，扩展版给出 Bézout 系数从而求出模逆元；\(ax\equiv b\pmod n\) 有解当且仅当 \(\gcd(a,n)\mid b\)；中国剩余定理把大模数拆成两两互素的小模数；拉格朗日定理推出欧拉定理与费马小定理；反复平方法 \(O(\beta)\) 次模乘计算模幂；RSA 的正确性来自费马小定理加中国剩余定理，安全性来自分解困难；Miller–Rabin 通过寻找 1 的非平凡平方根修补了费马测试，每轮误判率不超过 1/2，但解读结果时要做贝叶斯修正。字符串匹配部分：Rabin–Karp 用滚动哈希做快速过滤，KMP 用前缀函数在 \(\Theta(n+m)\) 内完成匹配并可流式运行。计算几何部分：叉积只用加减乘就能判断转向和相交；扫描线把二维问题变成按事件排序的一维过程；Graham 扫描 \(O(n\lg n)\) 求凸包；分治最近点对 \(O(n\lg n)\)；极大层就是帕累托分层。

练习

基础

1. 手算 \(\mathrm{EXTENDED\text{-}EUCLID}(899,493)\)，写出每层的 \((d,x,y)\)。（原书 31.2-2。答案 \(d=29\)。）
2. 求解 \(35x\equiv10\pmod{50}\) 的全部解。（原书 31.4-1。答案 \(6,16,26,36,46\)。）
3. 求被 9、8、7 除分别余 1、2、3 的最小正整数。（原书 31.5-2。答案 \(10\)，模 \(504\)。）
4. 对 \(p=11,q=29,e=3\) 完整走一遍 RSA：求 \(d\)，加密 \(M=100\)，再解密。（原书 31.7-1。答案 \(d=187\)，密文 254。）
5. 计算模式 \(ababbabbabbababbabb\) 的前缀函数。（原书 32.4-1。用 31.10.3 节的代码核对。）
6. 在 \(T=3141592653589793\) 中用 Rabin–Karp（模 \(q=11\)）找 \(P=26\)，会遇到多少个伪命中？（原书 32.2-1。）
7. 说明为什么 Graham 扫描检查的是"非左转"而不只是"右转"。（提示：共线的中间点不是凸多边形的顶点。）

进阶

8. 写出从右到左扫描指数二进制位的模幂算法，并说明它与原书算法的循环不变式有何不同。（原书 31.6-2。）
9. 证明：若 \(x\) 是模 \(n\) 的 1 的非平凡平方根，则 \(\gcd(x-1,n)\) 与 \(\gcd(x+1,n)\) 都是 \(n\) 的非平凡约数。据此说明 Miller–Rabin 在什么情况下能顺便给出一个因子。（原书 31.8-3。）
10. 用 \(T'T'\) 中查找 \(T\) 的方法，在线性时间内判断 \(T\) 是否是 \(T'\) 的循环移位。（原书 32.4-7。）
11. 设计 \(O(n^2\lg n)\) 算法判断 \(n\) 个点中是否有三点共线，并解释为什么对每个点做一次极角排序就够了。（原书 33.1-4。）
12. 在 31.10.4 节中，把随机组合的数量从 2 万改为 2000 和 20 万，观察凸包链与优化前沿的差距如何变化；再把资产数从 6 改为 30，解释差距为何显著变大。
13. 31.10.2 节中，若同时检验 1000 个候选因子、其中 10 个为真，用 Benjamini–Hochberg 方法控制 FDR 在 10%，模拟比较它与"\(p<0.05\) 直接筛选"的假发现比例。

原书推荐习题：31.2-2，31.4-1，31.5-2，31.6-2，31.6-3，31.7-1，31.8-3，31.9-2，思考题 31-1（二进制 gcd）、31-3（矩阵快速幂求斐波那契）；32.2-1，32.2-4，32.3-1，32.4-1，32.4-3，32.4-7，32.4-8；33.1-3，33.1-7，33.1-8，33.2-3，33.3-2，33.4-2，思考题 33-2（极大层）。

原书对照