第一次认真算圆周率是在一台只有 4GB 内存的旧笔记本上。我原本以为写个循环、套个公式就能跑出几百位结果程序跑到第十万项时还在小数点后第六位打转风扇声大得像要起飞。那一刻我才真正理解一件事圆周率计算的演进本质上不是“谁算得更准”而是“谁找到了让误差衰减得更快的结构”。从割圆术到超级计算机上的 Chudnovsky 级数中间隔着的不是算力而是一代代人对收敛速度的执念。这篇内容适合三类人看一是想搞懂 π 计算原理、但被公式吓退的初学者二是准备动手实现高精度计算、却总在精度和性能上翻车的开发者三是对数值算法、大数运算、并行计算感兴趣想找一个具体切口练手的人。我会把割圆术的几何逻辑、级数算法的提速思路、现代快速算法的工程取舍以及我在实际代码里踩过的坑一条条拆开讲。不堆公式但关键推导会给不神化超算但真实瓶颈会说清楚。1. 割圆术的数学内核古代算法到底在逼近什么1.1 从正六边形开始的双侧夹逼割圆术的核心思想非常朴素用一个内接正多边形去“填”圆再用一个外切正多边形去“罩”圆圆的周长被夹在两者之间边数越多夹得越紧。中国古代的割圆术从正六边形起步因为正六边形的边长恰好等于半径这是唯一一个不用开方就能直接写出边长的正多边形起点选得很聪明。设圆的半径为 1内接正 n 边形边长为 a_n那么内接多边形的周长是 n·a_n它一定小于圆周长 2π外切正 n 边形的周长则大于 2π。两边一夹π 的范围就出来了。真正让割圆术能持续推进的是边长递推公式当边数从 n 翻倍到 2n 时新边长可以用旧边长通过半角关系算出来。用现代写法若内接正 n 边形边长为 a则内接正 2n 边形边长为a sqrt(2 - sqrt(4 - a²))这个式子看起来简单但它里面藏着两次开方。手工时代没有计算器每一次开方都要靠笔算或算筹逐位求边数每翻一倍开方次数就翻一倍。这就是割圆术“理论优美、工程沉重”的根源。注意很多资料只给面积递推不给边长递推导致读者以为割圆术只是“画多边形量周长”。实际上真正的算法价值在于递推关系它把一次全新的几何构造压缩成一次代数迭代。1.2 误差衰减的残酷比例割圆术的误差量级是 O(1/n²)。这句话翻译成人话边数翻一倍误差降到原来的四分之一也就是多拿到约 0.6 位十进制精度。想多要一位精度边数得乘以大约 3.16想多要两位边数得乘以 10。我把这个关系整理成一张表方便直观感受目标精度十进制位大致所需正多边形边数需要翻倍的次数3 位约 10² 量级几次6 位约 10³ 量级约 10 次10 位约 10⁵ 量级约 17 次20 位约 10¹⁰ 量级约 34 次从 3 位到 20 位只多了 17 位精度边数却从百级涨到百亿级。这就是割圆术的天花板它不是不能算而是每多一位都要付出指数级的几何增长代价。古代数学家能靠割圆术把 π 推到七八位已经是人力、耐心和算筹技术的极限。1.3 手工开方为什么是真正的瓶颈假设你要算正 12288 边形需要从正六边形开始翻倍 11 次。每一次翻倍要对每个边长做两次开方。表面上只有 11 轮但每一轮的开方精度必须高于目标精度否则误差会累积。手工开方到十几位有效数字一位熟练的计算者一天也推进不了多少。更麻烦的是误差传播。割圆术的递推里嵌套了减法2 - sqrt(4 - a²)当 a 很小时4 - a²接近 4开方后接近 2两个接近 2 的数相减有效数字会大量丢失。这在数值计算里叫“相消误差”是初学者最容易忽略的坑。古代计算者未必有这套术语但他们通过保留更多中间位数来对抗它本质上是同一件事。我后来用浮点数复现割圆术时就栽在这个减法上用双精度浮点算到正 4096 边形结果尾部已经开始飘。换成高精度十进制库、并把中间精度多留 20 位结果才稳定。这个经验很直接递推公式里的减法永远要预留保护位。2. 级数时代从慢到无法接受到可工程化2.1 那个“正确但没用”的交替级数数学史上有一个非常著名的 π 级数π/4 1 - 1/3 1/5 - 1/7 1/9 - ...它形式极简证明也不难但它是一个“反面教材”。这个级数的误差量级是 O(1/n)也就是说你想多拿一位精度就得多算大约 10 倍的项。想要 10 位精度需要约 10¹⁰ 项想要 20 位需要 10²⁰ 项。任何计算机面对这个量级都会沉默。我第一次实现它的时候跑了一千万项结果只有小数点后六位是对的。那一刻我意识到判断一个 π 算法好不好不看公式漂不漂亮只看每项能换回多少位精度。这个指标后来成了我评估所有数值级数的第一标准。2.2 梅钦类公式的拆解逻辑真正让级数法变得可用的一步是把 arctan 拆成几个小参数的组合。arctan 的泰勒展开是arctan(x) x - x³/3 x⁵/5 - x⁷/7 ...当 x 越接近 0收敛越快。每算一项误差大约乘以 x²。所以如果 x 1/5每项能换回约 1.4 位十进制精度如果 x 1/239每项能换回约 4.8 位。相比之下x 1 的交替级数每项只能换回零点几位差距是数量级的。经典的梅钦型公式长这样π/4 4·arctan(1/5) - arctan(1/239)它的价值不在于“用了 arctan”而在于把大参数拆成了两个小参数。拆分的依据是复数的辐角相加公式两个 arctan 的和可以合并成一个 arctan反过来一个大 arctan 也可以拆成若干小 arctan 的线性组合。找拆分方案本质上是在整数格点上做搜索让组合出来的参数尽量小、项数尽量少。我用 Python 做过一个粗暴的搜索枚举 1 到 100 的分母尝试两项组合看哪个能把 π/4 表示出来且参数最小。结果很快就落到 1/5 和 1/239 这一组。这也解释了为什么这个公式被手工计算时代反复使用它不是唯一解但它是当时算力条件下性价比极高的一组解。2.3 从手算到机械计算的分水岭梅钦型公式把“算 π”从几何构造变成了纯算术迭代这是一个关键转折。几何方法需要维护多边形每一步都有明确的几何意义级数方法只需要加、乘、除天然适合机械化。手工时代计算者用这类公式把 π 推到一百多位靠的是极其严密的表格和校对流程。到了电子计算机出现同样的公式直接被搬上去第一批程序跑出来的结果就是几千位。这里有一个常被忽略的事实早期计算机算 π主要目的不是要 π而是验证机器长时间运行的稳定性和算术单元的正确性。π 的位数是已知的算错了就说明硬件或程序有问题它天然是一个“标准答案”。3. 现代快速算法二次收敛与每项十四位3.1 算术几何平均为什么能二次收敛算术几何平均算法AGM是 π 计算史上的分水岭。它的迭代形式大致是取两个初始值一个算术平均一个几何平均然后反复迭代同时维护一个修正项。每迭代一次有效精度大约翻一倍。“精度翻倍”是什么意思第一轮可能只有几位正确第二轮十几位第三轮几十位第四轮上百位。也就是说达到一万位精度只需要十几轮迭代。这和割圆术、交替级数的线性增长完全不在一个维度。二次收敛的代价是每轮要做一次高精度开方。开方本身又依赖乘法和除法所以单轮成本不低。但在高精度场景下轮数少带来的收益远大于单轮成本。这也是为什么 AGM 在很长一段时间里是高精度 π 计算的主流选择之一。我实测过一个简化版 AGM用 50 位中间精度迭代 6 轮能稳定拿到约 40 位正确数字迭代 8 轮能到 100 位以上。中间精度必须比目标精度多留十几位否则最后一轮的开方误差会把结果拖垮。3.2 Chudnovsky 级数为什么成为超算首选Chudnovsky 级数的特点是每一项大约贡献 14 位十进制精度。这个数字很关键它意味着算一万位只需要约 700 项算一亿位只需要约 700 万项。对比交替级数的十亿项起步差距是几何级的。它的公式结构包含阶乘、幂次和线性项看起来吓人但工程实现时有固定套路把求和项拆成“分子部分”和“分母部分”全部用整数表示用二进制分裂binary splitting把求和组织成树形结构从叶子往根合并最后做一次大整数除法再乘上常数因子。二进制分裂的价值在于它把大量小整数的乘除转换成少量大整数的乘法而大整数乘法可以用快速算法加速。树形结构还天然适合并行左右子树可以独立计算。我在自己的机器上对比过两种写法直接逐项用高精度浮点累加算十万位要几十秒改成整数二进制分裂同样精度降到几秒。差距不在公式而在“什么时候做除法”。逐项累加每一项都要做高精度除法而二进制分裂把除法推迟到最后一次前面全是整数乘法这就是工程优化的核心思路。3.3 BBP 公式的独特价值直接取第 n 位BBP 类公式有一个别的算法没有的能力它可以在不计算前面所有位的情况下直接算出十六进制展开中第 n 位附近的数值。这个特性叫“逐位提取”。它的实现逻辑是把级数中第 n 位之后的贡献通过模运算和小数部分提取出来避免构造完整的大数。实际用途主要有两个一是校验——算完一大段 π 之后抽几个远处的位比对二是分布式计算——不同机器负责不同区段不需要从头算起。不过要注意BBP 提取的是十六进制位不是十进制位。想要十进制第 n 位不能直接套用。很多初学者在这里误会以为“能直接取第 n 位”就等于“能直接取十进制第 n 位”实际不是一回事。4. 用 Python 把三种算法跑通4.1 精度设置多留保护位是铁律不管用哪种算法第一条规则都是工作精度必须高于目标精度。我一般多留 15 到 25 位。原因有三个迭代过程中的舍入累积、开方和除法的尾部误差、以及最后转字符串时的截断。以 Python 的 decimal 模块为例目标是 1000 位我会把精度设成 1030 左右from decimal import Decimal, getcontext TARGET 1000 getcontext().prec TARGET 30如果只设成 1000最后十几位大概率是错的。这不是库的问题是浮点表示和舍入的必然结果。我见过不少人拿最后几位去比对已知 π发现对不上就怀疑公式写错了其实只是保护位不够。4.2 梅钦公式的最小实现梅钦公式适合入门因为它只需要一个 arctan 级数函数from decimal import Decimal, getcontext def arctan_inv(x): # 计算 arctan(1/x)x 为正整数 getcontext().prec 10 x Decimal(x) total term Decimal(1) / x x2 x * x n 1 sign -1 while True: term term / x2 delta term / (2 * n 1) if delta 0: break total sign * delta sign -sign n 1 getcontext().prec - 10 return total def pi_machin(): return 4 * (4 * arctan_inv(5) - arctan_inv(239))这段代码里有两个细节值得说。第一getcontext().prec 10是临时抬高精度函数返回前再降回来避免污染全局设置。第二循环终止条件是delta 0因为 decimal 在固定精度下足够小的项会被舍入成零这比预设迭代次数更稳。实测下来这套代码算 1000 位大约需要一两秒算 10000 位就要几分钟因为 arctan(1/5) 每项只贡献约 1.4 位一万位需要七千多项每项都是高精度除法成本线性增长。4.3 AGM 的迭代框架AGM 的实现重心在迭代次数和开方精度from decimal import Decimal, getcontext def pi_agm(digits): getcontext().prec digits 20 one Decimal(1) two Decimal(2) four Decimal(4) a one b one / two.sqrt() t Decimal(0.25) p one # 迭代次数每次大约翻倍留一点余量 import math iterations int(math.log2(digits)) 2 for _ in range(iterations): a_next (a b) / two b_next (a * b).sqrt() t - p * (a - a_next) ** 2 a, b a_next, b_next p * two return (a b) ** 2 / (four * t)迭代次数用log2(digits) 2估算是因为每轮精度大约翻倍目标位数取对数就是轮数下限。多留两轮是保险。这里最容易出错的地方是t的更新顺序必须用旧的a和新的a_next的差不能反过来否则结果会偏。4.4 Chudnovsky 加二进制分裂Chudnovsky 的实现要复杂一些核心是把求和写成整数比值import math def binary_split(a, b): if b - a 1: if a 0: Pab Qab 1 else: Pab (6*a - 5) * (2*a - 1) * (6*a - 1) Qab a * a * a * 640320**3 // 24 Tab Pab * (13591409 545140134*a) if a 1: Tab -Tab return Pab, Qab, Tab m (a b) // 2 Pam, Qam, Tam binary_split(a, m) Pmb, Qmb, Tmb binary_split(m, b) Pab Pam * Pmb Qab Qam * Qmb Tab Qam * Tmb Pam * Tam return Pab, Qab, Tab这段递归把区间一分为二分别算出三组整数再合并。最后def pi_chudnovsky(digits): digits_per_term 14.181647462725477 terms int(digits / digits_per_term) 2 P, Q, T binary_split(0, terms) from decimal import Decimal, getcontext getcontext().prec digits 20 pi (Q * 426880 * Decimal(10005).sqrt()) / T return pi14.18这个系数就是每项贡献的十进制位数用它反推项数。多留两项是防止边界误差。这套实现在普通笔记本上算十万位大概几秒算一百万位需要几十秒到几分钟瓶颈在大整数乘法和最后的开方除法。5. 超级计算机上的 π瓶颈根本不在公式5.1 大整数乘法才是真正的战场到了千万位、亿位级别公式本身已经不是重点重点是“两个几千万位的大整数怎么相乘”。朴素竖式乘法是 O(n²)位数翻倍时间翻四倍完全不可接受。工程上会用 Karatsuba、Toom-Cook再往上用基于快速傅里叶变换的乘法把复杂度压到接近 O(n log n)。这也是为什么高精度计算库的性能差异巨大。同样一个 Chudnovsky 公式用普通语言内置的大整数和用专门优化的任意精度库速度能差出几个数量级。我自己做过对比一百万位的同一段求和换一个底层实现更好的库时间从几分钟降到几十秒。公式没变变的是乘法。5.2 并行策略与串行瓶颈Chudnovsky 的求和天然可并行把区间切成若干段每段独立做二进制分裂最后把各段结果合并。合并本身是一次树形归并也可以并行。但最后那一步“大整数除法加开方”是串行的而且越到高位越慢。所以超算上的实际策略往往是计算阶段高度并行校验阶段交叉验证存储阶段压缩输出。真正难的是内存。一亿位十进制数如果按每字节存两位也要几十 MB到了万亿位内存和磁盘都是 TB 级。很多团队的瓶颈不是算不动而是存不下、读不回、校验不完。5.3 算 π 到底在验证什么超级计算机反复算 π外界常觉得是“秀肌肉”。从工程角度看它确实是一个极好的压力测试长时间高负载运行、大规模并行通信、大容量内存读写、浮点和整数运算单元的正确性全都能被这个任务压出来。更重要的是校验机制。高精度 π 的常用校验方法是用两种独立算法各算一遍比对尾部若干位或者用逐位提取公式抽查远处区段。如果两种完全不同的数学路径给出同样的结果出错概率就极低。这套“交叉验证”思路其实适用于所有高精度数值计算不只限于 π。6. 我踩过的坑与几条硬经验6.1 精度设够但迭代次数没跟上有一次我用 AGM 算五千位精度设了五千二结果只有前两千位对。排查了半天才发现迭代次数按log2(5000)算是 13 轮但我代码里写死成了 10 轮。精度够、轮数不够结果一样错。这件事之后我养成了一个习惯迭代次数不要写死用目标位数反推并且打印每轮的实际有效位数。有效位数可以用“相邻两轮结果的差值”估算差值越小说明收敛越充分。这个自检成本很低但能省掉大量盲目调试。6.2 整数除法截断悄悄吃掉低位Chudnovsky 的二进制分裂全程用整数看起来不会丢精度。但最后那一步除法如果写错顺序比如先除后乘或者中间某个地方用了整除低位就会被截断。我遇到过一次结果前一百位全对从第一百零一位开始整体偏移查了很久才发现是某个中间变量用了整数除法。提示整数运算的好处是没有舍入误差坏处是一旦写错除法顺序错误非常隐蔽。建议每算完一段就和已知 π 的前若干位比对不要等全部算完再查。6.3 只比对前几位等于没校验初学者常犯的错是算出结果和记忆里的3.1415926比一下前七位对了就以为成功了。实际上高精度计算最容易错的就是尾部。正确的做法是拿一段足够长的已知 π 做全量比对或者至少抽三处开头、中间、结尾。我现在的校验流程是固定的先用小规模一千位跑通全量比对再逐步放大规模每次放大都保留上一次的已知结果做交叉。这样一旦某一步出错能立刻定位到是精度问题、轮数问题还是公式问题。6.4 别用浮点库硬扛高精度有人用双精度浮点去算几万位 π这在原理上就不可行。双精度只有约 15 到 17 位有效十进制数字后面全是噪声。高精度计算必须用任意精度整数或十进制库这是硬门槛没有捷径。如果只是学习算法、验证公式几百位用 decimal 完全够用如果要上规模就得换更底层的高精度库并且认真对待乘法实现。选型的分界线大致是一万位以内通用十进制库够用十万位以上考虑专门的任意精度库百万位以上乘法算法和内存布局比公式更重要。最后分享一个我自己常用的小技巧写这类高精度程序时先把目标位数设成很小比如 50 位把三种算法都跑通、全量比对通过再以 10 倍为单位往上放大。每次放大只改一个参数记录耗时和正确位数。这个过程看起来很笨但它能让你清楚知道每一层规模下时间花在了公式、乘法还是校验上。