突破 100 万位圆周率：从 Machin 级数到 Chudnovsky 算法与二进制拆分

在前端网页里高精度计算圆周率 $\pi$，很多人的第一印象是经典的梅钦类公式（Machin-like formula）：

$\frac{\pi}{4} = 4 \arctan\left(\frac{1}{5}\right) - \arctan\left(\frac{1}{239}\right)$

结合 JavaScript 原生 BigInt 放大定点整数，通过泰勒级数展开：

$\arctan(x) = x - \frac{x^3}{3} + \frac{x^5}{5} - \frac{x^7}{7} + \dots$

在计算 2,000 位以内的 $\pi$ 时，Machin 公式极度轻量，几十毫秒即可秒开。然而，一旦用户将目标精度提升至 100,000 位乃至 1,000,000 位（百万位），Machin 公式便遭遇了不可逾越的算法天花板：由于其每次迭代仅能推进约 1.4 位有效数字，且大整数单步相除的复杂度呈 $O(N^2)$ 增长，百万位计算将在浏览器中拖长至数小时甚至导致内存溢出。

为了打破这一物理限制，本站 圆周率计算器 引入了自适应分阶混合引擎：低位数（$<15,000$ 位）保留 Machin 公式秒开；高位数（$\ge 15,000$ 位）自动切换至 Chudnovsky 楚德诺夫斯基超高阶级数 + 二进制拆分（Binary Splitting）+ 牛顿-拉夫逊任意精度开方。

这篇文章我们将深度拆解这一算法跃迁背后的数学原理与纯 TypeScript 工程实现。

---

1. Machin 公式与 Chudnovsky 级数的收敛阶对比

要在算法上产生质的飞跃，关键在于提高“级数每增加一项所能获取的有效十进制位数”。

Machin 公式：O(N^2) 的渐进墙

Machin 反正切级数的项展开为 $x^{2k+1} / (2k+1)$。对于 $x = 1/5$，以 10 为底的对数显示：

$\log_{10}(5^2) = \log_{10}(25) \approx 1.397$

即级数每增加一项，仅能贡献约 1.397 位精度。计算 $N$ 位需要 $k \approx N / 1.397$ 项。因为常规递推中每次迭代都要对大整数做一次长除法，算法的总时间复杂度为：

$\text{Time}_{\text{Machin}}(N) = O\left(N \cdot N\right) = O(N^2)$

当 $N = 1,000,000$ 时，$N^2 = 10^{12}$ 次操作在单线程 JavaScript 引擎中是不可接受的。

Chudnovsky 级数：每项产生 14.181 位精确数字

1987 年，Chudnovsky 兄弟基于拉马努金（Ramanujan）椭圆模函数提出了如下超高阶快速收敛级数：

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

注意到基数 $640320^3 = 262537412640768000$。每一项的衰减因子为：

$\log_{10}(640320^3) = \log_{10}(2.625 \times 10^{17}) \approx 14.181647$

这意味着：级数每计算一项，精确的 $\pi$ 小数位就瞬间增加 14.181 位！计算 1,000,000 位只需要约 $\approx 70,510$ 项级数。

---

2. 为什么不能直接循环累加？二进制拆分（Binary Splitting）解密

虽然 Chudnovsky 级数仅需 7 万项，但如果直接用循环逐项计算大整数分子分母：

$S = \sum_{k=0}^{K-1} \frac{A(k)}{B(k)}$

直接做分数加法需要频繁求通分大分母，大整数长除法与乘法的运算规模随 $k$ 呈不均衡增长，整体复杂度依然维持在 $O(N^2)$。

为了将复杂度压进 $O(N \log N^3)$，必须采用**二进制拆分（Binary Splitting）**分治算法。

分解为矩阵多项式三元组

我们将 Chudnovsky 级数的第 $k$ 项表示为通分项比例：

$\frac{T(k)}{B(k)}$

对区间 $[a, b)$ 定义三元组 $(P(a,b), Q(a,b), T(a,b))$，使其满足合并定理：

对于区间 $[a, b)$（其中 $m = \lfloor(a+b)/2\rfloor$）：

1. 基础叶子节点（当 $b = a + 1$ 时）：

   • $P(a, a+1) = (6a-1)(2a-1)(6a-5)$ （当 $a=0$ 时为 1）
   • $Q(a, a+1) = a^3 \cdot \frac{640320^3}{24} = a^3 \cdot 10939058860032000$ （当 $a=0$ 时为 1）
   • $T(a, a+1) = P(a, a+1) \cdot (545140134 a + 13591409)$ （当 $a$ 为奇数时符号翻转）

2. 树状分治合并方程（分治左右区间 $L=[a, m)$ 与 $R=[m, b)$）：

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

通过这种树状分解，大整数的乘法在树的底部是极小整数乘法，在树的顶部是大整数乘法。结合现代 JavaScript 引擎（如 V8）内部对 BigInt 乘法采用的 Karatsuba / Toom-Cook 算法，整个分治合并的代价被降至：

$\text{Time}_{\text{BinarySplitting}}(N) = O(N \log N \log \log N) \approx O(N \log N^3)$

---

3. 手写 BigInt 任意精度 Newton-Raphson 整数开方

分治合并求得最终的 $Q(0, K)$ 与 $T(0, K)$ 后，根据 Chudnovsky 公式：

$\pi = \frac{426880 \sqrt{10005} \cdot Q(0, K)}{T(0, K)}$

表达式中包含关键项 $\sqrt{10005}$。在定点整数放大了 $10^{N + \text{guard}}$ 倍后，必须在 BigInt 域下精确求解 $\lfloor \sqrt{10005 \cdot 10^{2(N+\text{guard})}} \rfloor$。

常规 Math.sqrt() 仅能提供 IEEE 754 双精度 53 位尾数（约 15 位十进制数），无法处理百万位大数。我们手写了基于 牛顿-拉夫逊（Newton-Raphson）迭代法 的任意精度整型开方算法：

/**
 * 使用牛顿迭代法求解 floor(sqrt(n))，n 为 BigInt
 * 迭代公式: x_{k+1} = (x_k + n / x_k) >> 1
 */
function bigintSqrt(n: bigint): bigint {
	if (n < 0n) throw new Error('Square root of negative number');
	if (n === 0n) return 0n;
	if (n < 4n) return 1n;

	// 使用位长度估算初始猜测值，使收敛步数最少
	const bitLen = n.toString(2).length;
	let x = 1n << BigInt(Math.ceil(bitLen / 2));

	while (true) {
		const nextX = (x + n / x) >> 1n;
		if (nextX >= x) return x;
		x = nextX;
	}
}

牛顿迭代法具有二次收敛性（Quadratic Convergence）——即每轮迭代有效精确位数翻倍。即便针对 $2,000,000$ 位的巨大整型，仅需约 20 余次迭代即可收敛到末位精确值。

---

4. 算法性能实测与 1.5 万位交叉点（Crossover）

既然 Chudnovsky 算法复杂度如此优秀，为什么不直接全量使用 Chudnovsky，而保留 Machin 公式？

工程实测表明：Chudnovsky 算法拥有较高的常数初始化开销。

下表为在 Core i7 / Apple M2 浏览器环境下的实际性能对比数据：

小数位数 $N$	Machin 耗时	Chudnovsky 耗时	优胜算法
200 位	0.8 ms	12.4 ms	Machin 秒开
2,000 位	14.2 ms	45.1 ms	Machin 秒开
10,000 位	320 ms	380 ms	Machin 微弱领先
15,000 位	850 ms	840 ms	交叉点（Crossover）
50,000 位	12,400 ms	2,150 ms	Chudnovsky 领先 5.7 倍
100,000 位	58,000 ms	4,800 ms	Chudnovsky 领先 12 倍
1,000,000 位	$> 1.5$ 小时 (冻结)	仅需约 8.5 秒	Chudnovsky 压倒性优势

因此，我们在引擎控制层设置了 $N = 15,000$ 的分界点：

• $N < 15,000$：执行逻辑简单、无额外内存分配的 Machin 反正切递推；
• $N \ge 15,000$：自动切入 Chudnovsky + Binary Splitting 架构。

---

5. 防主线程冻结与 14 位安全冗余（Guard Digits）

在浏览器单线程环境中计算百万位，如果同步执行分治树合并，主线程仍会短暂卡顿。为此，计算器建立了 4 阶段异步分块流水线：

1. Stage 1 (0%–10%)：计算二进制拆分树范围及参数分配；
2. Stage 2 (10%–70%)：递归二进制拆分合并（每合并一定层级抛出 await new Promise(r => setTimeout(r, 0)) 释放主线程，刷新 UI 进度条）；
3. Stage 3 (70%–95%)：执行 bigintSqrt(10005) 牛顿开方与最终大整数相除；
4. Stage 4 (95%–100%)：结果转换为十进制字符串并执行 1,000 位自动换行格式化。

同时，为了防止最后一步相除时的舍入误差向左传递，算法额外引入了 $\text{guard} = 14$ 位安全冗余（const extraDigits = 14），最终输出时做截断裁切，确保输出的百万位小数 100% 数学准确。

欢迎前往本站 圆周率 π 计算器 亲体验百万位极限计算与 CPU 性能测速。