企业官网建设流程全解析

这个系列写到今天，终于轮到完结篇了。前面几篇里，我们玩过蒙特卡洛撒点、莱布尼茨级数，也折腾过韦达公式，把 π 算到一千位左右已经有点吃力的感觉。说实话，这些方法用来理解“π 为什么能被计算出来”足够有趣，但论性能全是花拳绣腿。这篇我会换一个画风，不玩“多撒几个点”的策略，直接从马青公式起步，再用 Chudnovsky 算法把 π 推到百万位量级，最后用 BBP 公式和本地文本做双向验证，确保这个百万位不是算着好玩的。整个过程我尽量用 Python 写完，保持可复现性。适合谁看？写过几行 Python、想真正体验一把高精度计算的读者，或者跟着系列一路走过来、好奇最终篇怎么收尾的朋友。

1. 从蒙特卡洛到马青公式：为什么朴实方法注定打不到百万位

1.1 先说说为什么我会“抛弃”前几篇的方法

蒙特卡洛算 π 的思路很直白：往一个正方形里随机撒点，数一数落在内切圆里的比例，再乘以 4。这个方法的误差大约正比于 1/√N，也就是撒 10000 个点，大概只能稳定到 2 位小数；撒到一百万，也就 3 到 4 位。想攒出 10 位正确小数，理论上是 10^20 这个量级的采样，完全没戏。

莱布尼茨级数 π/4 = 1 - 1/3 + 1/5 - 1/7 + ... 更惨，误差约等于 1/(2N+1)。你要 10 位精度，大约要取 50 亿项，这还是在浮点误差不捣乱的前提下。就算你把级数加速方法全用上，离百万位的目标仍然差好几个数量级。

所以真正用来冲精度的公式，必须是“每一项能多给好几个有效数字”那种。这类公式的核心套路是反三角函数展开：arctan(1/x) 的泰勒级数里，每往后一项，数值量级会缩小 x² 倍。只要 x 够大，收敛就快得肉眼可见。

1.2 马青公式：一个两三百年前的实用派方案

马青公式是 1706 年提出的：

π/4 = 4·arctan(1/5) - arctan(1/239)

为什么它成立？我当时用正切加法定理手推了一遍：设 α = arctan(1/5)，先用二倍角公式得到 tan(2α)=5/12，再用一次二倍角得到 tan(4α)=120/119，最后算 tan(4α - arctan(1/239))，分子分母恰好约干净，结果等于 1。既然正切值是 1，角度就是 π/4。整个过程每一项都来自普通整数比例，没有半点魔法，所以公式本身可以直接写进代码里。

Python 里实现 arctan(1/x) 的级数很简单：

from decimal import Decimal, getcontext def arctan_1_over_x(x, terms): # 在外部已经设置好 getcontext().prec x_inv = Decimal(1) / Decimal(x) x_pow = x_inv total = Decimal(0) for k in range(terms): term = x_pow / (2 * k + 1) if k % 2 == 0: total += term else: total -= term x_pow *= x_inv * x_inv return total def machin_pi(digits): getcontext().prec = digits + 20 terms = int(digits / 1.3) + 50 # 因为 arctan(1/5) 每一项约增加 1.4 位数字 pi = 4 * (4 * arctan_1_over_x(5, terms) - arctan_1_over_x(239, terms)) return format(pi, 'f')[:digits + 2]

这段代码算 1000 位绰绰有余。arctan(1/5) 每项贡献大约 log10(25)=1.3979 位数字，arctan(1/239) 收敛更快，每项约 4.76 位。所以 1000 位只需要几百项。问题出在“需要多少项”之外：每次循环的x_pow都在和同样是千位级精度的x_inv * x_inv相乘，一万位时 Decimal 的低速乘除法就开始拖后腿，到百万位时几乎不可接受。

这也就是为什么终极篇必须换到 Chudnovsky 赛道的根本原因。

2. Chudnovsky算法：真正的主角是整数运算

2.1 公式的样子，以及我为什么不直接用 Decimal 硬算

Chudnovsky 算法长这样：

π = 426880√10005 / ( Σ_{k=0}^∞ (-1)^k (6k)!(13591409+545140134k) / ((3k)!(k!)^3·640320^{3k}) )

这个公式每计算一项，大约能贡献 14.181647462725477 位十进制数字。也就是说，想要 100 万位，只要算大约 70495 项就够了。单看项数，比马青公式还少一个量级，但直接写成 Decimal 循环蹲坑的话，分母里那个 640320^{3k} 会膨胀成一个天文数字，每算一项都要做一次真正的“大数除法”，循环七万次，光运算时间就会让人崩溃。

正确姿势是二元分割（binary splitting）。它的核心思想：与其一项一项用小数除法，不如把整个区间 [0, n) 的求和先揉成一个有理数 T/Q，最后一次做高精度除法。揉的过程中只用整数乘除法，完全没有浮点数污染。最终的常数、开方、除法加起来，才是真正烧精度的时间点。

2.2 二元分割的 Python 实现

我直接给出完整脚本。这个实现受了 mpmath 里 Chudnovsky 模块的启发，改成独立可跑的版本：

import sys from decimal import Decimal, getcontext # 640320^3 // 24，提前算好，避免重复大整数运算 C3_OVER_24 = 640320**3 // 24 def chudnovsky_bs(a, b): """返回区间 [a, b) 的 (P, Q, T)，使系列部分和 = T / Q。""" if b - a == 1: if a == 0: P = Q = 1 else: P = (6*a - 5) * (2*a - 1) * (6*a - 1) Q = a * a * a * C3_OVER_24 T = P * (13591409 + 545140134 * a) if a & 1: T = -T return P, Q, T m = (a + b) // 2 Pam, Qam, Tam = chudnovsky_bs(a, m) Pmb, Qmb, Tmb = chudnovsky_bs(m, b) P = Pam * Pmb Q = Qam * Qmb T = Qmb * Tam + Pam * Tmb return P, Q, T def chudnovsky_pi(target_digits): # target_digits 表示小数点后要保留的正确位数 getcontext().prec = target_digits + 50 terms = int(target_digits / 14.181647462725477) + 2 P, Q, T = chudnovsky_bs(0, terms) sqrt_10005 = Decimal(10005).sqrt() pi = Decimal(426880) * sqrt_10005 * Q / T pi_str = format(pi, 'f') return pi_str[:target_digits + 2] # 加 2 是因为整数位 '3' 和小数点 if __name__ == "__main__": sys.setrecursionlimit(1000000) print(chudnovsky_pi(1000)[:120])

注意第 20 到 22 行的合并逻辑：

• P 是分子部分的乘积；
• Q 是整套分母因子的乘积；
• T 是“加了符号的累计项”，合并时需要对跨越中点两边的项做交叉乘加。

这一手其实就是把柯西求和的分配律用代码写了一遍。合并之后 T/Q 正好是区间内所有项的精确有理和。最后一步Decimal(426880) * sqrt_10005 * Q / T，把这个有理数转到 Decimal 再除，精度由getcontext().prec控制。

提示：不要在这个阶段用math.sqrt(10005)。math.sqrt返回的是 double，只有 15 位有效数字，放在百万位计算里等于在源头撒了一把沙。

2.3 实测下来，位数和时间的关系

在我自己的普通笔记本上，粗测的体感时间是这样的，机器性能差异很大，只当一个量级参考：

百万位能在分钟级跑出来，已经是二元分割加 Python 大整数的功劳了。如果换回马青公式那套 Decimal 循环，这个时间至少要翻一个数量级。

3. 校验闭环：BBP公式和已知文本，怎么信任一个百万位的答案

3.1 第一层校验：前 100 位必须对

算出一个 100 万位的字符串，第一件事永远是拿前 100 位对一遍，因为哪怕算法里的符号写反一个，前 20 位也会错得离谱。标准参考值是：

3.14159265358979323846264338327950288419716939937510 58209749445923078164062862089986280348253421170679

把chudnovsky_pi(100)的结果和这段对比，基本能排除八成的低级错误。但注意，Decimal 的截断和四舍五入可能导致最后 1 位不同，所以对比前不要对末尾最后 2 位吹毛求疵，后面我会讲为什么需要冗余位。

3.2 第二层校验：BBP公式可以隔空取数字

如果只信任前 100 位，那中间几万位是不是算错了？最稳的做法是用 BBP 公式。它的亮点是：不需要算出 π 的全部前序数字，就能直接提取第 n 位十六进制小数。公式是：

π = Σ_{k=0}^∞ (1/16^k) [ 4/(8k+1) - 2/(8k+4) - 1/(8k+5) - 1/(8k+6) ]

要取第 n 位，核心是算模幂：把每一项中的 16^{n-k} 拆成“整数倍 + 余数”，整数倍部分不影响小数，余数部分除以分母，才是我们需要的。下面是一个直接可用的提取脚本：

def pow_mod(base, exp, mod): result = 1 % mod base %= mod while exp > 0: if exp & 1: result = (result * base) % mod base = (base * base) % mod exp >>= 1 return result def bbp_hex_digit(n): # 返回 π 的十六进制小数点后第 n 位，n 从 0 开始 getcontext().prec = 80 total = Decimal(0) for k in range(n + 1): m = 8 * k + 1 total += Decimal(4 * pow_mod(16, n - k, m)) / Decimal(m) m = 8 * k + 4 total -= Decimal(2 * pow_mod(16, n - k, m)) / Decimal(m) m = 8 * k + 5 total -= Decimal(pow_mod(16, n - k, m)) / Decimal(m) m = 8 * k + 6 total -= Decimal(pow_mod(16, n - k, m)) / Decimal(m) for k in range(n + 1, n + 30): factor = Decimal(1) / (Decimal(16) ** (k - n)) total += factor * ( Decimal(4) / Decimal(8 * k + 1) - Decimal(2) / Decimal(8 * k + 4) - Decimal(1) / Decimal(8 * k + 5) - Decimal(1) / Decimal(8 * k + 6) ) frac = total - int(total) return int(frac * 16)

快速验证：因为 π 的十六进制开头是3.243F6A8885A308D313198A2E...，所以小数点后第 0 位是十六进制的 2，第 1 位是 4。跑一下bbp_hex_digit(0)应该得到 2，bbp_hex_digit(1)应该得到 4。能对上，说明提取逻辑没写歪。然后你可以拿它抽查你自己计算结果的第 5000、第 50000、第 500000 位附近，再把十六进制结果和 Chudnovsky 算出的十进制字符串统一进制转换后对比。两边独立跑了不同的数学流程，能对上，基本可以放心。

提示：我这里的 BBP 脚本不是性能版本。n 到十万以上时，模幂循环也要跑很久，我只用它“抽查关键位置”，不会拿它逐位验证一百万位。

3.3 第三层校验：和权威文本做切片对比

最笨也最可靠的方法，是下载一份别人已经算好的 π 文本文件，比如 100 万位或更多，然后把我的输出和它做切片对比。

做法是：把 Chudnovsky 结果写进pi_mine.txt，用现成的 diff 类工具和权威文本比对。如果整个文件太大，不想全量比对，可以取几个固定偏移量，比如第 1、10000、100000、999000 位附近各截 50 个字符，两边一致就算通过。我自己的经验是至少取三段：开头一段、中间一段、接近末尾一段。因为文件末尾最容易出问题，那里往往是截断和进位边界。

4. 工程化问题：内存、精度和时间，藏在“算出来”背后的账

4.1 为什么不要在循环里做 Decimal 除法

我在第一次改造代码时，天真地把每一项直接算成 Decimal 再累加。结果是 1 万位之后速度陡降，最后实在忍不了。根源是：Decimal 的除法不是免费的，尤其是除数为超长整数时，内部要走复杂的修正算法，比整数乘法贵很多。二元分割把七万项揉成一次除法，是质的飞跃。

另外，Decimal 的prec一旦设置，并不会自动给中间计算加保护位。很多人只设置prec = target_digits，最后几位立刻开始波动。我现在的习惯是至少加 50 个保护位。如果你想要 100 万位，可以把prec设成target_digits + 100，牺牲一点速度换稳定。

4.2 内存的账要心里有数

二元分割的代价是内存。看代码你会发现整个递归过程中，P、Q、T 三个大整数在每一层都存在。百万位时，Q 和 T 的十进制位数不是 100 万，而是大约 120 万到 130 万，因为 640320^{3k} 这一项导致分子分母本身膨胀得比 π 还长。Python 的int在每个二进制 limb 上还有对象开销，实际跑百万位时峰值内存可能到几百 MB。如果觉得吃力，可以用迭代版的二叉堆合并来压缩峰值，但代码可读性会差很多。

4.3 什么时候应该从 Python 切到 GMP 或 C++

Python 的int底层虽然也是 C 实现的，但它的乘法算法在大整数场景下不如 GMP 激进。想继续往千万位冲，基本三条路：

• 用gmpy2在 Python 里继续做，它直接把 GMP 的mpz、mpfr暴露给 Python，能把瓶颈顶上去不少；
• 换 C++ 加 boost::multiprecision::cpp_int，自己实现二元分割，性能优势和自由度都更大；
• 追求世界纪录级别的，直接去用 y-cruncher，它内部用了 FFT 乘法等大杀器，Python 这个系列到这里就该体面退场了。

我自己的建议是：Python 版能跑通百万位，已经足够理解 π 计算的完整链条。真要继续冲千万位，第一件事不是调 Python 代码，而是先把数据结构和乘法库换了。

5. 收尾：我这次实操中最后悔没早点知道的事

5.1 三个让我白跑几小时的错误

第一是math.sqrt。第一版我把Decimal(10005).sqrt()写成了math.sqrt(10005)，结果前 30 位一点问题没有，越靠后错得越离谱。这个错误的隐蔽性在于：它不会让你“完全算错”，而是让你“后面全错”，尤其是验证前 100 位时几乎不可能发现。后来我把校验点放到第 5000 位附近，才暴露出问号。

第二是对末尾截断的执念。π 的十进制展开没有循环，任何有限项逼近都会让最后 1 到 2 位出现进位抖动。我第一次切字符串时直接取前N + 2个字符，导致第 N 位偶尔少 1。现在的做法是：用两个不同保护位数（比如 +50 和 +100）各算一遍，比较结果的前 N 位；稳定不变的那些位才是可信的。

第三是输出环节。百万位字符串直接print到终端，不仅慢，搞不好终端直接卡死。正确做法是写到文件，然后对文件做校验。你以为算法是瓶颈，最后发现 I/O 和输出格式才是，这种乌龙我遇到不止一次。

5.2 给想复现的人的最后一点建议

如果你是第一次跑这套流程，强烈建议从 1000 位开始，把马青公式、Chudnovsky 公式、BBP 验证三段全部串通，再跳到 10 万位，最后才冲 100 万位。每一档之间改的东西只有目标位数和prec，如果某一档结果对不上，能迅速定位是算法问题还是环境问题。直接一上来跑百万位，出了问题反而不好查。

另外，把公式里的 640320、13591409、545140134 这些常数当作文物一样对待，不要随手“优化”成浮点数。它们存在的原因和共轭数、模方程有关，我们在这个系列里不需要深挖数论，但至少要知道一句话：这些数字不是拍脑袋来的，少一个因子，结果就会往错误方向狂奔。

我这个系列从撒点到级数，再到今天的二元分割，把 π 从小数点后几十位一路推到百万位，整个过程最大的体会是：算 π 难的不是“列出公式”，而是别让精度、截断、进位和内存这些工程细节在半路把你绊倒。如果你也打算复现，希望这篇能让你少走几圈我走过的弯路。