π这个东西大家从小就知道约等于3.14但真正把它算到小数点后几百位、几千位甚至更高的时候事情就没那么浪漫了。我最初接触高精度计算π倒不是想挑战世界纪录而是因为在一次算法竞赛的题目里需要用到几百位的π作校验值当时用Python的math.pi一查发现只给了十几位瞬间就卡住了。后来认真研究了一轮才发现“高精度”这三个字背后的门道远比想象中深涉及大整数乘法、傅里叶变换、级数收敛分析、浮点误差控制甚至还有内存和CPU的取舍。这篇文章就把我折腾高精度π计算的完整过程讲透包括算法怎么选、大数怎么算、代码怎么写、坑怎么填。适合参加算法竞赛的同学、对数值计算感兴趣的开发者以及任何想验证自己“能算多少位π”的偏执狂。1. 高精度算π整体思路与算法选型1.1 常用算法的收敛速度与实现成本提到计算π很多人第一反应是“随机撒点”的蒙特卡洛法或者莱布尼茨级数。这两种方法教学演示没问题但拿来算高精度纯属自我折磨。蒙特卡洛法的误差收敛速度是O(1/√N)也就是说要精确到小数点后10位至少要10^20个随机点这谁能顶得住。莱布尼茨级数π/4 1 - 1/3 1/5 - 1/7 ...更惨每一项大概只能贡献0.2位有效数字算到第1000项都只能勉强到小数点后三位属于典型的“眼睛会了、电脑废了”。稍微靠谱一点的是马青公式π/4 4arctan(1/5) - arctan(1/239)马青公式是早期手工计算和简单程序计算π的常用方案每一项大概能贡献1.4位小数收敛速度线性但系数还算友好。用这个公式算几千位是可行的但要到百万位就力不从心了。再往上走有几条主流路线一是高斯-勒让德算法AGM利用算术几何平均迭代收敛速度是二次的每迭代一次有效位数翻一倍。这个算法在数学上很漂亮实现起来需要高精度的平方根、除法和加法。二是Chudnovsky算法这是目前绝大多数世界纪录保持者使用的公式它每计算一项就可以稳定提供14位左右的有效数字效率极高但递推式里有复杂的组合数运算。1.2 为什么实战中我优先选择Chudnovsky公式高斯-勒让德虽然收敛快但它依赖高精度浮点数运算每一步都要处理开方和除法中间很容易引入舍入误差而且对数值稳定性要求很高。Chudnovsky公式属于整数级数求和公式长这样1/π 12 Σ(-1)^k (6k)! (13591409 545140134k) / ((3k)! (k!)^3 640320^(3k3/2))看着唬人实际操作时可以把根号640320的部分提取到求和外面统一处理中间就只剩下整数乘除和累加。只要控制好大整数运算的精度结果天然稳定不需要反复做浮点修正。这就是为什么Chudnovsky在工程实现里反而更顺手。从项数估算上也能看出优势每项大约贡献14.18位十进制数要算10万位只需要7050项左右。对比马青公式要算7万多项差距非常明显。1.3 精度目标决定实现路径这里必须实话实说不同数量级的位数技术方案完全不同。算几千位直接用Python的decimal模块就行几秒出结果。算几十万位必须自己实现大整数乘法至少得上FFT否则时间会长到让人怀疑CPU坏了。算千万位以上需要多线程、内存池、针对FFT的优化甚至要写汇编级别的乘法基本是专业玩家的领域。我自己的实践是从1000位开始的逐步往百万位走。每次往上跳一个数量级都会遇到新问题。这篇文章主要讲千位到百万位之间的实现路线这个区间是大多数算法爱好者能体验到的“高精度乐趣区”。2. 高精度计算的地基大数运算与数值稳定性2.1 为什么不能直接用浮点数很多人一开始会想我不就多算几位小数吗把浮点数精度调高不就行了吗问题是现代CPU和编程语言里的浮点数遵循IEEE 754标准双精度double只有53位二进制有效位换算成十进制只有15到17位有效数字。你就算把long double搬出来也就二三十位想算到几千位纯属白日梦。解决方案有两种一是用专门的任意精度浮点库比如MPFR、Boost.Multiprecision它们在底层把数字拆成多个机器字来存储支持自定义精度二是直接用大整数运算把问题转化为整数计算最后再通过一次高精度除法变成小数输出。第二种方案优点在于误差可控中间过程不涉及浮点舍入所以我个人更推荐在高精度算π时走整数计算路线。Python里的int天生就是任意精度所以用Python做原型验证非常爽但性能一般。真要追求速度C配合GMP库或者自己写FFT才是王道。2.2 大整数乘法从竖式到FFT大整数运算里最关键的是乘法基数越大乘法效率越成为瓶颈。小学学的竖式乘法复杂度是O(n²)如果算10万位的π中间要做的乘法规模可能达到几十万位的整数乘法O(n²)完全扛不住。Karatsuba算法把复杂度降到O(n^1.585)思路是把大数拆成高低两半通过三次小规模乘法凑出结果。这个算法实现简单中等规模的数据下效果很好。再往上还有Toom-Cook算法把大数拆成三部分甚至更多复杂度进一步下降。但当数字位数到几十万甚至上百万时最靠谱的是FFT乘法。核心思想是把两个大整数看成两个多项式计算它们的卷积这正是多项式乘法的本质。卷积在时域是O(n²)但通过离散傅里叶变换DFT转换到频域两个多项式的乘法就变成了对应位置的点乘再逆变换回来即可。这个过程复杂度是O(n log n)比O(n²)快了几个数量级。这里可以打个比方竖式乘法就像两个人隔着一条街互相传球两人手里的每个球都要和对方每个球碰一次FFT就像把所有人拉到一个操场上统一报数报完数再各回各家效率高得多。2.3 除法与平方根的牛顿迭代处理高精度算π的过程中除法和开方也逃不掉。Chudnovsky公式最后一步需要做C / S其中C和S都是几十万位的大数如果直接做长除法效率低得吓人。常见优化是先求出除数的倒数再做乘法而求倒数用到的是牛顿迭代法设f(x) 1/x - D迭代公式为x_{n1} x_n(2 - D*x_n)每迭代一次精度翻倍这和高斯-勒让德算法的收敛特性如出一辙。开平方也一样求1/√a可以用牛顿迭代x_{n1} x_n(3 - a*x_n²) / 2算完后再乘一次a就是√a。这些操作在高精度计算里几乎是标配。实际写代码时如果你用的是成熟的GMP库直接调用mpn_divrem和mpn_sqrtrem就行不需要自己造轮子。但如果是为了学习和竞赛亲手实现一遍牛顿迭代会非常有帮助。3. 实操用Chudnovsky公式算到万位3.1 精简代码与核心递推逻辑下面这段代码是基于Chudnovsky公式的经典实现用Python的decimal模块处理任意精度浮点数逻辑比较直观适合理解公式到代码的映射过程。我直接用它可以算到1万位再往上就慢很多但不影响理解原理。from decimal import Decimal, getcontext def compute_pi(digits): # 比目标位数多留10位减少末尾舍入波动 getcontext().prec digits 10 # Chudnovsky 公式中的常数 C 426880 * Decimal(10005).sqrt() # 初始化递推项 L Decimal(13591409) K Decimal(6) M Decimal(1) X Decimal(1) S L # 累加和 # 每项约提供14.18位有效数字 term_count int(digits / 14.181647462725476) 2 for i in range(1, term_count 1): # 更新 M对应组合数部分的递推 M M * (K**3 - 16 * K) / (i**3) # 更新 L L 545140134 # 更新 X对应 640320^(3k) 项带上正负号 X -X * 262537412640768000 # 累加当前项 S M * L / X # K 每次增加12 K 12 pi C / S return pi if __name__ __main__: result compute_pi(10000) print(str(result)[:10012])这段代码核心就是把级数求和拆成了四个递推量M对应组合数部分L对应线性项X对应幂次项K对应递增的下标。每一步之间都是简单的乘除和加减但每一步涉及的大数位数都在膨胀所以归根结底还是大整数运算的比拼。3.2 运行结果与位数验证我在一台普通笔记本上跑了这段代码算1万位大约需要十几秒结果开头的数字是3.14159265358979323846264338327950288419716939937510…和公认值完全一致。注意因为decimal的精度只比目标位数多10位所以必须在打印结果做截断把尾部可能不稳的几位藏掉。这就是所谓“留出冗余精度”的工程细节。如果想在C里算几十万位就绝对不能依赖类似decimal的通用浮点库而是要把整个求和过程变成纯整数递推最后只做一次高精度除法。对中间步骤还要用FFT优化的乘法。3.3 从万位到百万位的工程优化从1万位往上走最先遇到的是时间瓶颈。Python的decimal是经过优化的C库但依然不够快。我测试时算5万位大约要一分钟十几万位就要二十多分钟了。这里有几个能明显提升性能的手段把中间的所有Decimal运算改成整型运算也就是大整数运算避免浮点除法的开销。大整数乘法从内置算法换成Karatsuba或FFT后者通常在位数大于几千时就开始有优势。项数控制在digits / 14.181647 3左右不要盲目多加项否则白白增加计算量。使用多线程并行计算级数求和不同项之间可以分组求和最后再合并。我实际用C配合自定义FFT乘法把计算速度提升了上百倍20万位大约能在几秒内完成。这个路线就是目前主流高精度π计算库的基本框架。4. 常见问题与排查技巧实录4.1 为什么结果最后几位一直不对这是高精度计算π最典型的坑几乎每个人都会踩。我最初用decimal算的时候只给了和目标位数相同的精度结果后面几位总是往上飘。原因很简单级数求和的中间结果需要更大的中间精度最后计算C / S时一旦S或C的精度不足商的末尾位就会失真。解决方法也很粗暴计算时精度至少比目标位数多留10到20位算完再截断。为什么不是留1位因为除法的试商和舍入可能在最高位产生连锁影响留位数太少不够安全。我实际测试下来留10位在1万位以内几乎不会出错。4.2 计算速度慢如何定位性能瓶颈高精度计算是个天然的性能分析标本。我在优化时习惯用性能剖析工具先看哪一行代码耗时最长。经验是绝大多数时间花在乘法上尤其是大数之间的乘法其次是高精度除法。如果乘法是O(n²)的实现位数一上去就原形毕露。建议把乘法单独抽出来做基准测试分别测几百位、几千位、几万位的耗时画个趋势线马上就能看出来是自己的乘法太慢还是别的环节在拖后腿。4.3 怎么确认算出来的π是准的验证结果这件事反而比计算本身还重要。我常用的方法有三种和已知结果对比从网络上找一个公认的π值常量文件比对前几百位或几千位最简单直接。利用BBP公式抽查BBP公式可以直接算出十六进制表示下某一位的值不需要前面所有位。比如我想验证第100万位的十六进制值用BBP公式针对性计算再和自己算出的二进制结果转换后比对即可。双算法交叉验证用Chudnovsky公式和AGM算法分别算一遍如果两者在目标位数内完全一致基本可以断定结果正确。最后一种方法我特别推荐因为如果两个独立算法在相同位数下结果吻合几乎排除了系统性的算法实现错误。5. 实测数据与环境建议5.1 不同位数在Python下的耗时实测我在自己的开发机上跑了一组测试机器配置是Intel i5-12400内存32GBPython 3.11使用decimal实现。结果仅供参考目标位数级数项数实测耗时1,000730.3秒10,00070615秒50,0003526约80秒100,0007051约300秒可以看出Python的decimal在万位级别已经明显吃力了10万位直接到5分钟级别。这个表现并不是说算法不行而是通用浮点环境下的实现冗余太多。同样的位数如果换成C配合GMP或者自写FFT通常能在几秒内跑完。5.2 环境与依赖建议如果你只是体验一下直接装Python 3.8以上版本就能跑上面的代码不需要额外依赖。如果打算挑战百万位级别推荐以下工具组合之一C GMP MPFR工业级方案稳定高效。C NTT数论变换用模素数替代浮点FFT避免浮点误差。Python gmpy2在Python生态里最接近GMP性能的选择。我个人建议先把Python的decimal版本跑通理解公式和递推关系再加优化。一上来就上C和FFT很容易被一堆底层细节淹没反而不利于理解高精度计算的核心思想。5.3 参考资源与延伸话题如果你对高精度计算π真正产生了兴趣后续还可以研究BBP公式。BBP公式是一个十六进制下的级数展开式它最神奇的地方是可以跳过前面的所有位直接算出任意指定位置附近的十六进制数字这对验证π结果的正确性非常有用。再进一步可以了解GMP底层的大整数乘法实现它涵盖了Karatsuba、Toom-Cook、FFT等多种算法根据数字大小自动选择最优策略值得精读。另外很多超频爱好者和硬件测评博主喜欢用SuperPi这类工具跑分它的底层算法就是上面这套东西高精度π计算也因此被拿来评估CPU的单线程性能和内存稳定性。我个人在实际操作中最大的体会是高精度计算π不是一个孤立的小玩具它是大整数运算、数值分析、算法优化三个方向的综合演练场。一个看起来只和数学常数有关的题目最后逼着你去理解FFT、牛顿迭代、复杂度分析和并行优化这种收获远比“算出了多少位π”本身更重要。如果你也正好在研究这块建议从小位数开始亲手把代码写一遍再逐步把位数往上推我相信你会在这个过程里踩到比我更多、也更精彩的坑。