【从 163 到一千万位：用 Chudnovsky 算法计算 π 的完整实践】

原创已于 2026-08-10 22:33:58 修改 · 388 阅读

6 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

收录于

Python

数学之美

仿真验证

于 2026-08-10 22:17:37 首次发布

从 163 到一千万位：用 Chudnovsky 算法计算 π 的完整实践

摘要：本文完整记录一次"用 Chudnovsky 算法把 π 算到 10,000,000 位"的实践：从算法根由（163 与模形式的魔法、Ramanujan 级数、二分分裂）到 60 行 gmpy2 实现（8 秒、86MB、三重独立验证），再到与 mpmath / decimal / Machin 公式 / y-cruncher 的基准对决。文末附全部源码与可复现步骤。期间还踩到一个差点毁掉基准数据的 mpmath 缓存陷阱，一并记录。

0. 结果速览

指标	数值
位数	10,000,000（“3.” + 1000 万位小数）
算法	Chudnovsky 公式 + 二分分裂 + GMP 整数平方根（gmpy2/mpz）
耗时	≈ 8 s（7.99–8.23 s，单线程；二分分裂占大头）
峰值内存	≈ 86 MB
环境	Windows 11 / Python 3.11.9 / gmpy2 2.3.1 / Intel i7-1255U (10C/12T)
验证	mpmath 交叉验证 ✅、angio.net 官方 1M 参考逐字节一致 ✅、50M 二进制参考前 10M 位 SHA-256 相同 ✅
SHA-256	`46059c61a4de67d6c916fa958168789da324a03ee8a85c30e9ca292c3712eb25`
统计结论	单数字 χ²=2.78 (p=0.97)、数字对 χ²=91.7 (p=0.69)、三位组 χ²=1071 (p=0.057)，均不拒绝均匀分布假设

一句话：在 2023 年的普通笔记本上，单核 Python 8 秒算完的 1000 万位，已经超过了 1997 年日本超级计算机（日立 SR2201，29 小时算 515 亿位）的全程平均速度。这就是算法与硬件双重复利的缩影。

1. 算法根由：π 计算的百年军备竞赛

1.1 从割圆术到 Machin 公式

π 的计算史是一部"收敛率"进化史：

多边形逼近（阿基米德，前 3 世纪；刘徽割圆术，3 世纪；祖冲之密率 355/113，5 世纪）：每翻一倍边数只多出约 0.5 位精度，想多要 100 位需要天文数字的计算量。
arctan 级数（Madhava 14 世纪，Gregory-Leibniz 17 世纪）： $\arctan x = x - \frac{x^3}{3} + \frac{x^5}{5} - \cdots$ 。取 $x = 1$ 就是莱布尼茨级数，每项只贡献约 0.3 位；但取 $x = 1/5$ 时每项贡献 $\log_{10}25 \approx 1.4$ 位。
Machin 公式（1706）： $\frac{\pi}{4} = 4\arctan\frac{1}{5} - \arctan\frac{1}{239}$ 这是"小参数 arctan 组合"的鼻祖，把 π 算到 100 位，统治了此后 200 多年（Rutherford、Shanks 都是它的变体；Shanks 1873 年算到 707 位，其中 180 位是错的，直到 1946 年才被发现）。
电子计算机时代：ENIAC 1949 年 70 小时算 2037 位；1961 年 IBM 7090 算 100,265 位；1973 年 CDC 7600 算 100 万位（23 小时 18 分）——用的还是 Machin 类公式，复杂度 $O(n^2)$ ，100 万位基本到顶。

1.2 Ramanujan 与 Chudnovsky：来自模形式的馈赠

1914 年，Srinivasa Ramanujan 在论文 Modular equations and approximations to π 中给出了一族计算 $1/\pi$ 的超几何级数，其中最著名的是：

$\frac{1}{\pi} = \frac{2\sqrt{2}}{9801}\sum_{k=0}^{\infty}\frac{(4k)!\,(1103+26390k)}{(k!)^4\,396^{4k}}$

每项贡献约 8 位——比 Machin 快了近 6 个数量级的收敛速度。

1989 年，David 与 Gregory Chudnovsky 兄弟在 Ramanujan 理论的基础上给出了收敛更快的变体（即本文主角，每项 14.18 位），并用它加上一台用邮购零件在自己公寓里攒出来的并行机 “mzero”，先后打破 4.8 亿位（1989）、10.1 亿位（1989）、22.6 亿位（1991）的纪录，把 CRAY、日立等超级计算机实验室甩在身后——这是计算数学史上最著名的"车库逆袭"之一。

1.3 163 的魔法：j-不变量与"几乎整数"

Chudnovsky 级数里那个神秘的 $640320^3$ 从哪来？答案是模形式。

考虑 $j$ -不变量，它是上半平面上的模函数，有 Fourier 展开：

$j(\tau) = q^{-1} + 744 + 196884q + \cdots,\qquad q = e^{2\pi i \tau}$

取 $\tau = \frac{1+\sqrt{-163}}{2}$ ，则 $-e^{-\pi\sqrt{163}}$ 。关键事实：

163 是最大的 Heegner 数。Heegner 数（1, 2, 3, 7, 11, 19, 43, 67, 163）恰好是类数为 1 的虚二次域对应的判别式（Heegner 1952 年证明完备性）。
类数 1 意味着 $j(\tau)$ 是代数整数——对判别式 $- 163$ 而言，它精确地等于一个普通整数：

$j\left(\frac{1+\sqrt{-163}}{2}\right) = -640320^3 = -262537412640768000$

于是由 $j(\tau) \approx q^{-1} + 744$ （ $q$ 小到 $10^{-18}$ 量级）：

$e^{\pi\sqrt{163}} \approx 640320^3 + 744 \approx 262537412640768743.99999999999925\ldots$

这就是著名的 Ramanujan 常数（“几乎整数”）。1975 年 Martin Gardner 在 Scientific American 的愚人节专栏里煞有介事地宣布 $e^{\pi\sqrt{163}}$ 是整数，骗过了无数读者。

Chudnovsky 公式里的 $640320^3$ 正是这个 $j$ -值——级数本身来自 Ramanujan 的 $1/\pi$ 模方程理论。数学的"无用之美"（模形式）就这样直接变成了 2025 年 314 万亿位纪录的引擎。

1.4 每项 14.18 位是怎么算出来的

设级数项 $a_k = \dfrac{(-1)^k (6k)!\,(13591409+545140134k)}{(3k)!\,(k!)^3\,C^k}$ ， $C=640320^3$ 。用 Stirling 近似：

$\frac{(6k)!}{(3k)!\,(k!)^3} \sim \text{const}\cdot k^{-3/2}\cdot 1728^k$

因此相邻项之比：

$\frac{a_{k+1}}{a_k} \longrightarrow \frac{1728}{640320^3} \approx 6.58\times 10^{-15}$

每项贡献的十进制位数：

$-\log_{10}\frac{1728}{640320^3} = \log_{10}\frac{640320^3}{1728} = 14.181647462725477\ldots$

对比：Machin 每项约 1.4 位，Ramanujan 9801 系列约 8 位，Chudnovsky 14.18 位。

1.5 二分分裂：把 $O(n^2)$ 压成 $O(M(n)\log n)$

级数收敛快还不够——逐项累加 70 万个 1000 万位的分数，代价是 $O(n^2)$ 的。二分分裂（binary splitting） 是破局点：

把部分和表示为有理数 $S = T (0, N) / Q (0, N)$ ，其中：

$P(a,b)=\prod_{k=a}^{b-1}P_k,\quad Q(a,b)=\prod_{k=a}^{b-1}Q_k,\quad T(a,b)=\sum_{k=a}^{b-1}P(a,k)\,T_k\,Q(k+1,b)$

分治合并：

$P(a,b)=P(a,m)P(m,b),\quad Q(a,b)=Q(a,m)Q(m,b),\quad T(a,b)=Q(m,b)T(a,m)+P(a,m)T(m,b)$

每层合并是若干次大数乘法，共 $\log n$ 层 → 总复杂度 $O(M(n)\log n)$ ，其中 $M (n)$ 是 $n$ 位大数乘法的复杂度（GMP 对超大数使用 Schönhage-Strassen FFT 乘法）。

一个关键工程事实：CPython 内置 int 的乘法只到 Toom-Cook（无 FFT），而 GMP 有 FFT——所以同样的二分分裂，用 CPython int 算 1000 万位需要数小时，用 gmpy2（GMP 绑定）只要 8 秒。这也是本实现选择 gmpy2 的根本原因。

1.6 为什么现在纪录都用 Chudnovsky 而不是 AGM

1976 年 Brent–Salamin 发现的高斯-勒让德（AGM）算法是二次收敛（位数翻倍式增长），曾是 1980 年代超级计算机的主力。但它在 FFT 时代反而劣势：每次迭代都要做大数平方根，而 Chudnovsky 全程只需要一次平方根。自 y-cruncher 2009 年问世以来，所有 π 世界纪录（10T→62.8T→105T→314T）全部基于 Chudnovsky + 二分分裂。

2. 数学推导：公式整理与整数化

Chudnovsky 公式的标准形式：

$\frac{1}{\pi} = 12\sum_{k=0}^{\infty}\frac{(-1)^k (6k)!\,(13591409 + 545140134k)}{(3k)!\,(k!)^3\, 640320^{3k+3/2}}$

整理成便于二分分裂的形式。令 $C = 640320^3 = 262537412640768000$ ：

$\pi = \frac{426880\,\sqrt{10005}}{S},\qquad S = \sum_{k=0}^{\infty}\frac{(-1)^k (6k)!\,(13591409+545140134k)}{(3k)!\,(k!)^3\,C^k}$

（推导要点： $640320^{3/2} = 640320\cdot 8\sqrt{10005}$ ，与分母的 12 合并得 $426880\sqrt{10005}$ 。）

整数化（全程无浮点）：要得到 $\pi$ 的 $D$ 位十进制数字，计算

$\lfloor \pi\cdot 10^{D}\rfloor = \left\lfloor \frac{426880 \cdot \mathrm{isqrt}(10005\cdot 10^{2D})\cdot Q(0,N)}{T(0,N)} \right\rfloor$

其中 $\mathrm{isqrt}$ 是 GMP 的整数平方根（误差 < 1 ulp，由保护位吸收）。取 $D = 10{,}000{,}000 + 1000$ ，算完截断前 1000 万位即可；级数项数 $\lceil D/14.181647462725477\rceil + 2 = 705{,}209$ 。

3. 工程实现

3.1 代码结构

核心 bs(a,b) 实现二分分裂，两个工程要点：

叶节点块化（CHUNK=64）：不递归到单个项（70 万个 Python 函数调用太贵），而是每叶直接累算 64 项。实测把 10M 位耗时从 21.70s → 8.23s（2.6 倍加速），且结果 SHA-256 完全一致（可复现性验证）。
峰值内存采样：psutil 后台线程每 200ms 采样 RSS，10M 位峰值约 86MB。

3.2 验证策略（三重独立 + 可复现）

验证	范围	结果
mpmath 1.4.1（独立 Chudnovsky 实现）	前 100,000 位	✅ 一致
angio.net 官方 `pi1000000.txt`	前 1,000,000 位	✅ 逐字节一致
angio.net `pi50.4.bin`（50M 位 BCD 打包）	前 10,000,000 位	✅ SHA-256 完全相同
两次独立运行（不同分裂实现 v1/v2）	全量	✅ SHA-256 相同

单靠"跑得快"是不够的，高精度计算必须用独立实现交叉验证——这也是 y-cruncher 世界纪录必须附带"验证文件"（用另一套公式独立复核）的原因。

4. 统计分析：数字规律

对 1000 万位做频率、卡方、游程、子串分析（analysis.py）：

4.1 单数字频率（χ²=2.784, df=9, p=0.972）

数字	次数	占比	数字	次数	占比
0	999,440	9.9944%	5	1,000,466	10.0047%
1	999,333	9.9933%	6	999,337	9.9934%
2	1,000,306	10.0031%	7	1,000,207	10.0021%
3	999,964	9.9996%	8	999,814	9.9981%
4	1,001,093	10.0109%	9	1,000,040	10.0004%

最大偏差仅 0.011 个百分点（3σ 界为 ±0.0285 个百分点）——均匀得"过分"。

4.2 高阶结构

数字对：χ²=91.7 (df=99, p=0.69)，最高 08（100,816 次）/ 最低 32（99,314 次），期望 100,000。
三位组：χ²=1071 (df=999, p=0.057)，最高 242（10,296 次）/ 最低 067（9,632 次）。p 值接近 5% 但只是 1.6σ 波动，不显著。
游程：最长连续 7 位（0000000 起于第 3,794,572 位），10M 位内无 8 连；Feynman 点 999999 在第 762 位 ✅（与已知值一致）。
特殊子串：今天日期 20260810 未出现（期望仅 0.1 次），但其 7 位前缀 2026081 恰好出现在第 9,382,034 位；0123456789 未出现（已知首次出现在第 17,387,594,880 位，远超 10M）。

结论：所有检验均不拒绝"独立均匀随机"假设。π 的十进制展开在已检验尺度下符合正态数（normal number）行为——"没有规律"本身就是答案；那些"巧合"（Feynman 点、日期串）都是随机性的正常产物。

5. 基准对标

5.1 同机实测（i7-1255U，同一进程内顺序测量，含十进制转换、不含写盘）

实现	1k	10k	100k	1M	10M
本实现（gmpy2 手写 Chudnovsky，块=64）	0.000s	0.001s	0.029s	0.435s	7.99s（125 万位/秒）
mpmath 1.4.1（gmpy2 后端）	0.001s	0.008s	0.039s	0.481s	8.53s（117 万位/秒）
decimal 模块（C 实现，无 FFT）	0.001s	0.058s	—	—	50k=0.42s（10M 预计小时级）
Machin 公式（1706，decimal）	0.003s（930 项）	—	—	—	5k=0.033s（4,631 项）

扩展性： $T(10\mathrm{M})/T(1\mathrm{M}) = 18.4$ （理论超线性 $n^{1.1}$ 量级，与 $O(M(n)\log n)$ 一致）。

要点：

本实现与 mpmath 基本持平（同为 Chudnovsky + gmpy2，实现细节差异约 ±7%）；块化优化后反超 mpmath 约 7%。
1000 位所需级数项数：Machin ≈930 项 vs Chudnovsky ≈72 项——收敛率差距的直观体现。
decimal 模块没有 FFT，50k 位就到 0.42s，1000 万位是"小时级"任务。

5.2 与 y-cruncher 对比（官方基准表，2024-07 更新，含十进制转换）

CPU	25M 位耗时	折合 10M 位（按 n^1.1 外推）
Core i3 8121U（低功耗笔记本）	1.951 s	≈0.75 s
Core i7 11800H（笔记本）	0.490 s	≈0.19 s
Ryzen 9 7940HS（笔记本）	0.410 s	≈0.16 s
Ryzen 9 7950X（桌面）	0.287 s	≈0.11 s
Ultra 9 285K（桌面）	0.206 s	≈0.08 s

y-cruncher（多线程 + AVX-512 + 手写汇编）比我们的 Python 实现快约 15–100 倍——这符合预期，也说明 Python 生态的合理定位是"验证与教学"，追求极限速度需要 C/汇编路线。

5.3 历史坐标（全程均值速度）

年份	纪录	机器	速度
1949	2,037 位	ENIAC（70h）	0.008 位/秒
1961	100,265 位	IBM 7090（8.7h）	3 位/秒
1973	1,000,000 位	CDC 7600（23h18m）	12 位/秒
1997	515 亿位	日立 SR2201（29h）	49 万位/秒
2002	1.24 万亿位	日立 SR8000（602h）	57 万位/秒
2009	2.7 万亿位	单台桌面机 i7 920（Bellard，131 天）	24 万位/秒
2019	31.4 万亿位	Google Cloud（121 天）	300 万位/秒
2024	105 万亿位	2× AMD Epyc 9754（75 天）	1620 万位/秒
2025-11	314 万亿位	2× AMD Epyc 9965（110 天）	3300 万位/秒

5.4 插曲：一个差点毁掉基准数据的缓存陷阱

第一版 benchmark_compare.py 中，mpmath 的 10M 位测出了 1.414s——比单独运行（8.4s）快 6 倍，而且位数还是对的。排查后真相是：

mpmath 的 constant_memo 装饰器把缓存存在闭包内部的原始函数上（f.memo_prec / f.memo_val），而不是包装函数 g 上。我执行的 L.pi_fixed.memo_prec = -1 只是给 g 设了个无效属性——重置从未生效。于是第一次调用（交叉验证阶段）算完 10M 位后，基准循环阶段全部命中缓存，1.414s 只是 str() 十进制转换的耗时。

正确重置要绕道闭包：

g = libelefun.pi_fixed
f = g.__closure__[0].cell_contents   # 原始 pi_fixed
f.memo_prec = -1
f.memo_val = None

教训：基准测试前先确认"被测代码真的在干活"，缓存命中测出来的"性能"毫无意义。

6. 完整源码

6.1 计算脚本 `pi_chudnovsky.py`

# -*- coding: utf-8 -*-
"""
Chudnovsky 算法 + 二分分裂 (Binary Splitting) 计算 pi 的十进制数字 (gmpy2/mpz)。

核心公式:
    1/pi = 12 * sum_{k=0..oo} (-1)^k (6k)! (13591409 + 545140134k) / ((3k)! (k!)^3 * 640320^(3k+3/2))
    化简: pi = 426880 * sqrt(10005) / S
          S  = sum_k (-1)^k (6k)! (13591409+545140134k) / ((3k)! (k!)^3 * C^k),  C = 640320^3
    二分分裂: P(a,b), Q(a,b), T(a,b) 递推, S = T(0,N)/Q(0,N)
    每项约贡献 log10(640320^3/1728) = 14.181647462725477 位十进制精度
整数化 (避免浮点):
    floor(pi * 10^D) = floor(426880 * isqrt(10005 * 10^(2D)) * Q / T)
    (isqrt 为 GMP 全精度整数平方根, 误差 < 1 ulp, 由 extra 保护位吸收)
"""
import sys, time, hashlib, threading, argparse


def compute_pi_scaled(digits, extra=1000, verbose=True):
    """返回 s = floor(pi * 10^(digits+extra)) 的十进制字符串, 以及分阶段耗时统计."""
    import gmpy2
    from gmpy2 import mpz, isqrt

    C3_24 = mpz(10939058860032000)          # 640320^3 / 24
    sys.setrecursionlimit(100000)

    D = digits + extra
    N = int(D / 14.181647462725477) + 2     # 需要的项数
    CHUNK = 64                              # 叶节点块大小: 每叶直接累算 CHUNK 项

    def bs(a, b):
        if b - a <= CHUNK:
            Pk = []; Qk = []; Tk = []
            for k in range(a, b):
                if k == 0:
                    p = q = mpz(1)
                else:
                    p = (6 * k - 5) * (2 * k - 1) * (6 * k - 1)
                    q = C3_24 * k * k * k
                t = p * (13591409 + 545140134 * k)
                if k & 1:
                    t = -t
                Pk.append(p); Qk.append(q); Tk.append(t)
            Pab = mpz(1)
            for p in Pk:
                Pab *= p
            Qab = mpz(1)
            for q in Qk:
                Qab *= q
            m = len(Pk)
            Qsuf = [mpz(1)] * (m + 1)
            for i in range(m - 1, -1, -1):
                Qsuf[i] = Qsuf[i + 1] * Qk[i]
            Tab = mpz(0); Ppre = mpz(1)
            for i in range(m):
                Tab += Ppre * Tk[i] * Qsuf[i + 1]
                Ppre *= Pk[i]
            return Pab, Qab, Tab
        m = (a + b) // 2
        Pam, Qam, Tam = bs(a, m)
        Pmb, Qmb, Tmb = bs(m, b)
        return Pam * Pmb, Qam * Qmb, Qmb * Tam + Pam * Tmb

    # 峰值内存采样线程
    peak = [0]
    stop = threading.Event()
    def sampler():
        try:
            import psutil
            p = psutil.Process()
            while not stop.is_set():
                m = p.memory_info().rss
                if m > peak[0]:
                    peak[0] = m
                time.sleep(0.2)
        except Exception:
            pass
    th = threading.Thread(target=sampler, daemon=True)
    th.start()

    t0 = time.perf_counter()
    P, Q, T = bs(0, N)
    t1 = time.perf_counter()
    tenD = mpz(10) ** D
    sqrtC = isqrt(mpz(10005) * tenD * tenD)   # floor(sqrt(10005)*10^D)
    t1b = time.perf_counter()
    num = mpz(426880) * Q * sqrtC
    t1c = time.perf_counter()
    pi_scaled = num // T
    t2 = time.perf_counter()
    s = str(pi_scaled)
    t3 = time.perf_counter()
    stop.set()
    th.join(timeout=2)

    if verbose:
        print(f"[bs   ] terms={N}  {t1 - t0:.2f}s", flush=True)
        print(f"[sqrt ] {t1b - t1:.2f}s  [mul] {t1c - t1b:.2f}s  [div] {t2 - t1c:.2f}s  [str] {t3 - t2:.2f}s", flush=True)
        print(f"[total] {t3 - t0:.2f}s  peakRSS={peak[0] / 1e6:.0f}MB", flush=True)

    assert s[0] == '3' and len(s) == D + 1
    stats = {"bs": t1 - t0, "sqrt": t1b - t1, "mul": t1c - t1b, "div": t2 - t1c,
             "str": t3 - t2, "total": t3 - t0, "terms": N, "peak_rss": peak[0]}
    return s, stats


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("digits", type=int)
    ap.add_argument("out", nargs="?", default=None)
    ap.add_argument("--extra", type=int, default=1000)
    ap.add_argument("--verify", action="store_true",
                    help="用 mpmath 交叉验证前 min(100000, digits) 位")
    args = ap.parse_args()

    s, stats = compute_pi_scaled(args.digits, args.extra)
    digits_str = s[1:1 + args.digits]
    content = "3." + digits_str
    assert len(digits_str) == args.digits

    if args.verify:
        import mpmath as mp
        n = min(args.digits, 100000)
        mp.mp.dps = n + 20
        ref = str(mp.pi)[2:2 + n]
        if digits_str[:n] == ref:
            print(f"[verify] mpmath 交叉验证前 {n} 位: 一致 OK", flush=True)
        else:
            for i, (a, b) in enumerate(zip(digits_str[:n], ref)):
                if a != b:
                    print(f"[verify] 不一致 @ 位置 {i}: ours={a} ref={b}", flush=True)
                    break
            sys.exit(1)

    sha = hashlib.sha256(content.encode("ascii")).hexdigest()
    print(f"[sha256] {sha}", flush=True)
    print(f"[head] {content[:101]}", flush=True)
    print(f"[tail] ...{content[-101:]}", flush=True)

    if args.out:
        t = time.perf_counter()
        with open(args.out, "w", encoding="ascii") as f:
            f.write(content)
        print(f"[save] {args.out}  {time.perf_counter() - t:.2f}s", flush=True)


if __name__ == "__main__":
    main()

6.2 分析比对脚本 `benchmark_compare.py`

（多实现基准 + 正确性交叉验证 + Machin 历史对照；含 mpmath 缓存陷阱的修复）

# -*- coding: utf-8 -*-
"""
benchmark_compare.py — 多实现 π 计算基准比对 + 正确性交叉验证

参与对比的实现:
  1. 本实现   : Chudnovsky 公式 + 二分分裂 (gmpy2/mpz), 块=64        [pi_chudnovsky.py]
  2. mpmath   : mpmath 1.4.1 官方 Chudnovsky (gmpy2 后端), 冷启动     [需重置闭包内缓存]
  3. decimal  : Python decimal 模块 (C 实现, 无 FFT) 的 Chudnovsky
  4. Machin   : π/4 = 4·arctan(1/5) − arctan(1/239) (1706 年公式, 历史对照)

输出: stdout 比对表 + benchmark_results.csv / benchmark_results.md
"""
import os, sys, time, csv

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import pi_chudnovsky as pc          # 本实现
from decimal import Decimal, getcontext


# ---------- 1. 本实现 ----------
def ours(digits):
    s, stats = pc.compute_pi_scaled(digits, extra=100, verbose=False)
    return s[1:1 + digits], stats["total"]


# ---------- 2. mpmath (强制冷启动) ----------
# 注意: mpmath 的 constant_memo 把缓存存在闭包内部的原始函数上 (f.memo_prec / f.memo_val),
# 对包装函数 g 赋值无效。必须通过 __closure__ 拿到内部函数才能真正重置缓存。
from mpmath.libmp import libelefun as _L
_MPMATH_PI_INNER = _L.pi_fixed.__closure__[0].cell_contents

def _cold_reset_mpmath_pi():
    _MPMATH_PI_INNER.memo_prec = -1
    _MPMATH_PI_INNER.memo_val = None

def mpmath_pi(digits):
    import mpmath as mp
    _cold_reset_mpmath_pi()                 # 真·冷启动
    mp.mp.dps = digits + 5
    t0 = time.perf_counter(); x = +mp.pi; t1 = time.perf_counter()
    s = str(x); t2 = time.perf_counter()
    return s[2:2 + digits], (t1 - t0, t2 - t1)     # (位数串, (计算s, 转换s))


# ---------- 3. decimal 模块 Chudnovsky ----------
def pi_decimal(digits):
    getcontext().prec = digits + 10
    C = Decimal(640320) ** 3
    C3_24 = C // 24
    def bs(a, b):
        if b - a == 1:
            if a == 0:
                Pab = Qab = Decimal(1)
            else:
                Pab = (6 * a - 5) * (2 * a - 1) * (6 * a - 1)
                Qab = C3_24 * a * a * a
            Tab = Pab * (13591409 + 545140134 * a)
            if a & 1:
                Tab = -Tab
            return Pab, Qab, Tab
        m = (a + b) // 2
        Pam, Qam, Tam = bs(a, m)
        Pmb, Qmb, Tmb = bs(m, b)
        return Pam * Pmb, Qam * Qmb, Qmb * Tam + Pam * Tmb
    N = int((digits + 10) / 14.181647462725477) + 2
    P, Q, T = bs(0, N)
    pi = (Decimal(426880) * Decimal(10005).sqrt() * Q) / T
    return str(pi)[2:2 + digits]


# ---------- 4. Machin 公式 (历史对照) ----------
def pi_machin(digits):
    """π/4 = 4·arctan(1/5) − arctan(1/239), arctan 用泰勒级数."""
    getcontext().prec = digits + 10
    def atan_inv(x):
        """arctan(1/x) = Σ (-1)^k / ((2k+1) x^(2k+1))"""
        term = Decimal(1) / x
        x2 = Decimal(x) * x
        s = term
        n = 1
        while True:
            term /= x2
            d = term / (2 * n + 1)
            s = s - d if n & 1 else s + d
            if abs(d) < Decimal(1).scaleb(-digits - 5):
                break
            n += 1
        return s, n + 1
    a5, n5 = atan_inv(5)
    a239, n239 = atan_inv(239)
    return str(16 * a5 - 4 * a239)[2:2 + digits], n5 + n239


# ---------- 交叉验证 ----------
print("=" * 78)
print("交叉验证 (正确性)")
print("=" * 78)
ours_10m, _ = ours(10_000_000)
mp_10m, _ = mpmath_pi(10_000_000)
dec_10k = pi_decimal(10_000)
mach_1k, mach_terms = pi_machin(1_000)
checks = [
    ("本实现 vs mpmath (10M 位)", ours_10m == mp_10m),
    ("本实现 vs decimal (10k 位)", ours_10m[:10_000] == dec_10k),
    ("本实现 vs Machin 公式 (1k 位)", ours_10m[:1_000] == mach_1k),
    ("mpmath vs decimal (10k 位)", mp_10m[:10_000] == dec_10k),
]
for name, ok in checks:
    print(f"  [{'PASS' if ok else 'FAIL'}] {name}")
assert all(ok for _, ok in checks), "交叉验证失败!"

# ---------- 基准 ----------
print()
print("=" * 78)
print("基准比对 (本机, 含十进制转换, 不含写盘)")
print("=" * 78)
SIZES_OURS_MP = [1_000, 10_000, 100_000, 1_000_000, 10_000_000]
SIZES_DEC = [1_000, 10_000, 50_000]
SIZES_MACHIN = [1_000, 5_000]

rows = []
def add(impl, digits, secs, extra=""):
    rate = digits / secs
    rows.append({"实现": impl, "位数": f"{digits:,}", "总耗时(s)": round(secs, 3),
                 "万位/秒": round(rate / 1e4, 1), "备注": extra})
    print(f"  {impl:<16} {digits:>10,} 位  {secs:>8.3f}s  {rate/1e4:>10.1f} 万位/秒  {extra}")

for d in SIZES_OURS_MP:
    _, t = ours(d)
    add("本实现(gmpy2)", d, t)
for d in SIZES_OURS_MP:
    _, (c, st) = mpmath_pi(d)
    add("mpmath", d, c + st)
for d in SIZES_DEC:
    t0 = time.perf_counter(); pi_decimal(d); t = time.perf_counter() - t0
    add("decimal", d, t)
for d in SIZES_MACHIN:
    t0 = time.perf_counter(); _, terms = pi_machin(d); t = time.perf_counter() - t0
    add("Machin(1706)", d, t, f"共 {terms:,} 项级数")

# Machin 项数对照 (1k 位)
_, t1k = pi_machin(1_000)
print(f"\n  [对照] 1000 位所需级数项数: Machin≈{t1k:,} 项 vs Chudnovsky≈{int(1000/14.181647462725477)+2} 项"
      f" (每项位数: Machin≈1.4, Chudnovsky≈14.18)")

# 扩展性
t10 = next(r["总耗时(s)"] for r in rows if r["实现"] == "本实现(gmpy2)" and r["位数"] == "10,000,000")
t1 = next(r["总耗时(s)"] for r in rows if r["实现"] == "本实现(gmpy2)" and r["位数"] == "1,000,000")
print(f"  扩展性: T(10M)/T(1M) = {t10/t1:.1f} (理论 ~ n^1.1)")

# ---------- 落盘 ----------
with open("benchmark_results.csv", "w", newline="", encoding="utf-8-sig") as f:
    w = csv.DictWriter(f, fieldnames=["实现", "位数", "总耗时(s)", "万位/秒", "备注"])
    w.writeheader()
    w.writerows(rows)
with open("benchmark_results.md", "w", encoding="utf-8") as f:
    f.write("| 实现 | 位数 | 总耗时(s) | 万位/秒 | 备注 |\n|---|---|---|---|---|\n")
    for r in rows:
        f.write(f"| {r['实现']} | {r['位数']} | {r['总耗时(s)']} | {r['万位/秒']} | {r['备注']} |\n")
print("\n已输出: benchmark_results.csv / benchmark_results.md")

6.3 统计分析脚本 `analysis.py`

（完整源码见项目目录；核心逻辑：Counter 统计频率/数字对/三位组，mpmath.gammainc 算卡方 p 值，单遍扫描游程，str.find/count 定位特殊子串，10×1M 分块稳定性。）

# -*- coding: utf-8 -*-
"""对 pi_10m_digits.txt 做数字规律统计分析: 频率、卡方检验、数字对、三位组、游程、特殊子串、分块统计."""
import json, sys
from collections import Counter

SRC = sys.argv[1] if len(sys.argv) > 1 else "pi_10m_digits.txt"
OUT = sys.argv[2] if len(sys.argv) > 2 else "stats.json"

raw = open(SRC, encoding="ascii").read().strip()
assert raw.startswith("3.")
D = raw[2:]                      # 去掉 "3." 的小数位字符串
n = len(D)
print(f"digits loaded: {n}", flush=True)

def chi2_pvalue(chi2, df):
    """用 mpmath 上不完全伽马函数算卡方 p 值."""
    import mpmath as mp
    return float(mp.gammainc(df / 2.0, chi2 / 2.0, regularized=True))

# ---------- 1. 单数字频率 ----------
freq = Counter(D)
expected = n / 10.0
chi2_digit = sum((freq[d] - expected) ** 2 / expected for d in "0123456789")
p_digit = chi2_pvalue(chi2_digit, 9)

# ---------- 2. 相邻数字对 ----------
pairs = Counter(D[i:i+2] for i in range(n - 1))
exp_pair = (n - 1) / 100.0
chi2_pair = sum((c - exp_pair) ** 2 / exp_pair for c in pairs.values())
chi2_pair += (100 - len(pairs)) * exp_pair
p_pair = chi2_pvalue(chi2_pair, 99)
top_pairs = pairs.most_common(10)
bot_pairs = sorted(pairs.items(), key=lambda kv: kv[1])[:10]

# ---------- 3. 三位组 ----------
triples = Counter(D[i:i+3] for i in range(n - 2))
exp_tri = (n - 2) / 1000.0
chi2_tri = sum((c - exp_tri) ** 2 / exp_tri for c in triples.values())
chi2_tri += (1000 - len(triples)) * exp_tri
p_tri = chi2_pvalue(chi2_tri, 999)
top_tri = triples.most_common(10)
bot_tri = sorted(triples.items(), key=lambda kv: kv[1])[:10]

# ---------- 4. 游程分析 ----------
max_run = {d: 0 for d in "0123456789"}
max_run_pos = {d: -1 for d in "0123456789"}
run_hist = Counter()   # 长度 -> 游程数
cur, cur_len, cur_start = D[0], 1, 0
for i in range(1, n):
    c = D[i]
    if c == cur:
        cur_len += 1
    else:
        run_hist[cur_len] += 1
        if cur_len > max_run[cur]:
            max_run[cur] = cur_len
            max_run_pos[cur] = cur_start + 1   # 1-based 小数位位置
        cur, cur_len, cur_start = c, 1, i
run_hist[cur_len] += 1
if cur_len > max_run[cur]:
    max_run[cur] = cur_len
    max_run_pos[cur] = cur_start + 1
overall_max_run = max(max_run.values())
overall_max_pos = max_run_pos[max(max_run, key=lambda d: max_run[d])]

# ---------- 5. 特殊子串 ----------
targets = [
    "20260810", "2026081", "202608", "2026",          # 今天的日期
    "314159", "271828", "161803", "141421", "577215", "693147",
    "123456", "1234567", "12345678", "0123456789", "9876543210",
    "000000", "111111", "222222", "333333", "444444",
    "555555", "666666", "777777", "888888", "999999",
    "9999999", "88888888", "13579", "24680",
]
substr = {}
for t in targets:
    first = D.find(t)
    substr[t] = {"first_pos": (first + 1) if first >= 0 else -1,  # 1-based
                 "count": D.count(t)}

# ---------- 6. 分块统计 (每 1M 位一块) ----------
blocks = []
for b in range(10):
    chunk = D[b * 1_000_000:(b + 1) * 1_000_000]
    cf = Counter(chunk)
    chi2_b = sum((cf[d] - 100_000.0) ** 2 / 100_000.0 for d in "0123456789")
    blocks.append({"block": b + 1, "range": f"{(b*1_000_000+1):,}-{(b+1)*1_000_000:,}",
                   "counts": {d: cf[d] for d in "0123456789"}, "chi2": chi2_b})

# ---------- 汇总 ----------
summary = {
    "digits": n,
    "freq": {d: freq[d] for d in "0123456789"},
    "freq_pct": {d: round(freq[d] / n * 100, 4) for d in "0123456789"},
    "chi2_digit": round(chi2_digit, 4), "p_digit": p_digit, "df_digit": 9,
    "chi2_pair": round(chi2_pair, 2), "p_pair": p_pair, "df_pair": 99, "exp_pair": exp_pair,
    "top_pairs": [[k, v] for k, v in top_pairs],
    "bot_pairs": [[k, v] for k, v in bot_pairs],
    "chi2_tri": round(chi2_tri, 2), "p_tri": p_tri, "df_tri": 999, "exp_tri": exp_tri,
    "top_tri": [[k, v] for k, v in top_tri],
    "bot_tri": [[k, v] for k, v in bot_tri],
    "max_run": max_run, "max_run_pos": max_run_pos,
    "overall_max_run": overall_max_run, "overall_max_pos": overall_max_pos,
    "run_hist": {str(k): v for k, v in sorted(run_hist.items())},
    "substr": substr,
    "blocks": blocks,
    "first_1000": D[:1000],
    "last_1000": D[-1000:],
}
json.dump(summary, open(OUT, "w", encoding="utf-8"), ensure_ascii=False, indent=1)

# ---------- 控制台摘要 ----------
print("\n=== 单数字频率 ===")
for d in "0123456789":
    print(f"  {d}: {freq[d]:,}  ({freq[d]/n*100:.4f}%)")
print(f"  chi2={chi2_digit:.4f} (df=9)  p={p_digit:.4f}")
print("\n=== 数字对 (top/bottom 5) ===")
print("  top:", [(k, v) for k, v in top_pairs[:5]])
print("  bot:", [(k, v) for k, v in bot_pairs[:5]])
print(f"  chi2={chi2_pair:.2f} (df=99)  p={p_pair:.4f}")
print("\n=== 三位组 (top/bottom 5) ===")
print("  top:", [(k, v) for k, v in top_tri[:5]])
print("  bot:", [(k, v) for k, v in bot_tri[:5]])
print(f"  chi2={chi2_tri:.2f} (df=999)  p={p_tri:.4f}")
print("\n=== 最长游程 ===")
for d in "0123456789":
    if max_run[d] >= 5:
        print(f"  {d*max_run[d]}  x{max_run[d]}  @ 位置 {max_run_pos[d]:,}")
print(f"  总体最长: {str(overall_max_run)} @ {overall_max_pos:,}")
print("\n=== 特殊子串 (首次出现位置 1-based / 总次数) ===")
for t in targets:
    fp, c = substr[t]["first_pos"], substr[t]["count"]
    print(f"  {t}: first@{fp if fp>0 else '无'}  count={c:,}")
print("\n=== 分块 chi2 ===")
for b in blocks:
    print(f"  块{b['block']}: chi2={b['chi2']:.3f}")
print("\nDONE ->", OUT)

6.4 其他文件

charts.py：4 张统计图（数字频率柱状图、数字对热力图、分块热力图、游程分布）→ charts/*.png
make_xlsx.py：8 工作表数据交付件 pi_stats.xlsx（概览/数字频率/数字对/三位组/游程/特殊子串/分块统计/基准对标）
stats.json：全量统计（机器可读）
pi_10m_digits.txt：最终结果（“3.” + 10,000,000 位）
benchmark_results.csv / benchmark_results.md：比对结果落盘

7. 复现与延伸

pip install gmpy2 mpmath psutil openpyxl matplotlib

# 计算 10M 位 (含 mpmath 交叉验证)
python pi_chudnovsky.py 10000000 pi_10m_digits.txt --verify

# 统计分析
python analysis.py pi_10m_digits.txt stats.json

# 多实现基准比对 + 交叉验证
python benchmark_compare.py

# 图表与 Excel
python charts.py && python make_xlsx.py

进一步加速的方向：多线程 FFT（GMP 的 mpn_fft 是单线程的，可并行化分治层）；换成 y-cruncher 的路线（AVX-512 + 手写汇编 + 磁盘交换，支持万亿位）；或 Bellard 2009 的 Fermat 数 FFT 方案。若只想要"现成最快"，直接用 y-cruncher——它官方基准里 10M 位在主流 CPU 上仅需 0.1–0.8 秒。

附件下载

本文涉及的完整源码和数据文件可通过以下链接下载：

文件	大小	说明
pi_chudnovsky.py	4.2 KB	主计算脚本（Chudnovsky算法+二分分裂）
benchmark_compare.py	5.8 KB	基准测试与交叉验证脚本
analysis.py	5.3 KB	统计分析脚本
pi_10m_digits.txt	9.5 MB	计算得到的1000万位π值（前1000位预览见下文）
stats.json	45 KB	完整的统计分析结果
pi_stats.xlsx	1.2 MB	Excel格式的统计报告（8个工作表）

使用说明：

下载所有文件到同一目录
确保已安装依赖：pip install gmpy2 mpmath psutil openpyxl matplotlib
运行 python pi_chudnovsky.py 1000000 计算100万位（10M位需要约8秒和86MB内存）

π的前1000位预览（完整1000万位见pi_10m_digits.txt）：

3.14159265358979323846264338327950288419716939937510
58209749445923078164062862089986280348253421170679
82148086513282306647093844609550582231725359408128
48111745028410270193852110555964462294895493038196
44288109756659334461284756482337867831652712019091
45648566923460348610454326648213393607260249141273
72458700660631558817488152092096282925409171536436
78925903600113305305488204665213841469519415116094
33057270365759591953092186117381932611793105118548
07446237996274956735188575272489122793818301194912
98336733624406566430860213949463952247371907021798
60943702770539217176293176752384674818467669405132
00056812714526356082778577134275778960917363717872
14684409012249534301465495853710507922796892589235
42019956112129021960864034418159813629774771309960
51870721134999999837297804995105973173281609631859
50244594553469083026425223082533446850352619311881
71010003137838752886587533208381420617177669147303
59825349042875546873115956286388235378759375195778
18577805321712268066130019278766111959092164201989

注意：由于CSDN附件大小限制（通常≤50MB），9.5MB的pi_10m_digits.txt可以正常上传。如需完整1000万位以外的更大数据文件，建议使用网盘分享。

8. 参考

D. V. Chudnovsky, G. V. Chudnovsky, The computation of classical constants, PNAS 86 (1989).
S. Ramanujan, Modular equations and approximations to π, Quart. J. Math. 45 (1914).
B. Haible, T. Papanikolaou, Fast multiprecision evaluation of series of rational numbers (1997).
R. P. Brent, Fast multiple-precision evaluation of elementary functions, JACM 23 (1976).
y-cruncher 官方纪录与基准: https://www.numberworld.org/y-cruncher/、http://www.numberworld.org/digits/Pi/
F. Bellard, 2009 年 2.7 万亿位纪录: https://bellard.org/pi/pi2700e9/
参考数字文件: https://www.angio.net/pi/digits.html（pi1000000.txt, pi50.4.bin）
Wikipedia: Chudnovsky algorithm, Chronology of computation of pi, Heegner number, Ramanujan constant