算圆周率这种事听起来特别像大学期末编程作业。但真把需求摆到台面上,目标是“小数点后10000位”性质就完全变了——它不再是你用Math.PI糊弄一下、截图交差的小练习而是实打实的高精度计算问题。你在网上搜一圈就会发现做这件事的人大致分三类一类是刚学完Python循环想看看自己能算到多少位一类是被面试题问到了比如“不用库怎么算pi到任意精度”还有一类是纯粹想验证某台机器或者某个语言的数值计算能力拿10000位当benchmark用。不管你是哪种这篇文章都打算把这件事彻底讲透从算法选型到浮点数陷阱从代码实现到末位校验全部走一遍。我会用Python为主顺便提一下如果是用C/C或者JavaScript思路要怎么变。1. 方案选型先想清楚再写循环1.1 为什么不能直接double乘到底先说一个最常见的坑。很多人的第一反应是圆周率不是有现成公式吗pi 4 * arctan(1)然后对着arctan的泰勒展开一项一项算不就行了理论上是这样但只要你真的用Python的float去跑算到小数点后十几位就开始不对劲到四五十位的时候已经完全不可信了。原因很简单IEEE 754的双精度浮点数只有53位有效二进制位换算成十进制大约是15到17位有效数字。你连第50位都算不出来更别提10000位。所以这个问题的第一课是如果你需要的小数位数超过15位一切浮点类型都要扔到一边必须用“高精度数”来做运算。Python里的decimal.Decimal、fractions.Fraction或者干脆用整数去模拟小数点都是路子。1.2 自己算还是用现成的高精度库很多人会问那直接用mpmath这种库不就完了一行mpmath.mp.dps 20000; print(mpmath.pi)搞定。这确实是最快的办法。但如果你正在准备面试、正在做算法作业或者只是想把计算过程彻底搞明白那用库反而会丢掉核心的乐趣。我自己倾向于把这当成两种梯度如果用库那么重点在于理解库背后用了什么公式、为什么快。如果自己写那么重点在于精度控制、存储设计和迭代收敛。我在这篇文章里会把两种路线都演示一遍先用Python标准库decimal配合一个收敛很快的公式拿到结果然后手写一个用整数运算实现的高精度计算版本最后对比校验。1.3 公式选哪个决定了你要写多少循环计算pi的公式非常多但并不是所有公式都适合算到10000位。Leibniz公式收敛慢到令人绝望算到几千项才几个有效位基本只能用来讲概念。Machin公式pi/4 4*arctan(1/5) - arctan(1/239)收敛快手工推导友好适合教学。每算一项大约能多出1.4位十进制精度。Chudnovsky算法收敛极快每项能多出约14位精度是目前各种破纪录计算的标准选项。打个比方如果Leibniz公式是走路Machin公式是骑自行车那Chudnovsky就是高铁。所以这篇文章里用decimal实现时我会直接用Chudnovsky算法自己手写整数版本时则用Machin公式——原因后面会说手写整数高精度除法和开方是有代价的Chudnovsky涉及大整数平方根初学者容易栽进去。2. Chudnovsky算法用Python标准库算出10000位2.1 核心原理Chudnovsky算法的公式长这样1/pi 12 * sum_{k0..∞} (-1)^k * (6k)! * (13591409 545140134k) / ((3k)! * (k!)^3 * 640320^(3k1.5))别被这串东西吓到。它本质上是一个收敛极快的级数。你每多算一项正确数字的位数大概能多出14位。这意味着算到10000位理论上只需要迭代到10000 / 14次左右也就是大约714项。实际计算时考虑到末尾截断误差你多迭代几十项就够了。这个公式能快速收敛的原因和模形式的深层数学性质有关不是这篇博客的重点。我们只需要把它当作一个“黑箱但透明”的迭代公式来用知道它快知道怎么实现知道每一项的计算方式即可。2.2 关键参数与预处理当我们要用decimal.Decimal来算时需要提前做几件事设置全局精度。如果要输出10000位不能只把精度设为10000因为计算过程中中间值有舍入误差一般要留出额外的保护位。我通常会设置到10000 50多的50位留给中间过程的舍入误差。把公式里的常数预先算出来不要放在循环里反复用高精度除法。比如640320**1.5、545140134、13591409这些循环外先算好。尽量用整数运算保存迭代过程中的递推关系而不是每次都从头计算阶乘。比如我们可以维护当前k对应的(6k)!、(3k)!、(k!)^3通过递推公式从k推进到k1这样能省掉大量乘法。2.3 代码实现下面这段代码是我实际测试过可以直接跑的版本from decimal import Decimal, getcontext # 计算位数 DIGITS 10000 # 多算一些保护位减少中间舍入误差 getcontext().prec DIGITS 50 def chudnovsky_pi(n_terms): 利用 Chudnovsky 级数计算圆周率 n_terms: 迭代项数每项约贡献 14 位精度 C 426880 * Decimal(10005).sqrt() M Decimal(1) L Decimal(13591409) X Decimal(1) K Decimal(6) S Decimal(13591409) for k in range(1, n_terms): M M * (K**3 - 16*K) / (Decimal(k)**3) L Decimal(545140134) X * Decimal(-262537412640768000) S Decimal(M * L) / X K Decimal(12) pi C / S return pi if __name__ __main__: pi chudnovsky_pi(800) print(str(pi)[:DIGITS 2]) # 包含整数部分的 3 和小数点这里n_terms我直接传了800。刚才说了大约714项够用多个几十项是为了让末尾更稳定。输出的时候因为getcontext().prec设置成了DIGITS 50得到的Decimal会多出一些尾部数字直接截断到我们需要的位数即可。2.4 为什么精度要设置成DIGITS 50很多第一次接触高精度计算的人会问10000位就10000位为什么要多算50位因为.这个公式虽然理论上每项贡献14位但每一步都涉及乘除法Decimal在运算中会对中间结果做舍入。如果精度卡在10000位整最后一次迭代的误差可能会污染最后几位。多出50位保护位等于给整个计算过程留了缓冲带。最终输出前再从第一位开始截取就不用担心末位抖动。有一个经验规律是保护位至少是迭代次数的两位数级别。我一般习惯取DIGITS / 100 20对于10000位来说就是120比较保守。但这里为了简洁用50也够。要注意的是保护位越多计算耗时越长因为高精度乘法耗时大约是位数的一次方到二次方之间没必要无限放大。2.5 输出后怎么验证输出一串数字很容易难的是确认这串数字是正确的。第一步先确认前几位是3.14159265。这能排除大部分错误。 第二步截取一段已知的正确圆周率片段做比对。比如我知道圆周率小数点后第1到20位是897932384626433832可以直接用字符串查找验证。 第三步更具说服力的办法是用两种完全不同的算法分别计算10000位再逐位比较。比如这一节的Chudnovsky算法算一遍下一节手写的基于整数和Machin公式的版本再算一遍两者比对。如果完全一致基本可以认定结果正确。3. 不用浮点数手写高精度整数运算算pi3.1 核心思路如果说上一部分你只是写了一个调用Decimal的脚本那这一部分才算真正理解高精度计算的底层原理。你要丢掉Decimal丢掉float只用Python原生的整数类型自己定义数组来模拟“任意精度实数”。一个很自然的想法是把pi表示成“整数部分 一个小数点后N位的整数”。比如我要算100位那我实际计算的是int(pi * 10^100)一个巨大的整数最后需要输出时我再在合适的位置插入小数点。这样一来所有高精度计算都转化成整数运算完全绕开浮点数精度限制。3.2 用整数实现高精度算术什么叫“用整数实现高精度算术”比如你要算两个“拥有10000位小数”的数之和最简单的方式是每个数用Python列表存每一位数字0-9。从最低位开始逐位相加记下进位。同理高精度乘法和除法都可以用类似手算的方法实现。Python的整数类型本身可以极大所以还有个取巧的做法不逐位存0-9而是把每9位数字拼成一个int作为“一坨”然后像数组元素一样处理从而大幅减少进位处理逻辑。对于10000位圆周率位数并不算夸张所以我直接用“整个数放大成一个大整数”的策略。比如我计算int(pi * 10^10000)这个数大约有10001位十进制数字Python的int处理起来毫无压力。3.3 基于Machin公式的整数实现Machin公式是pi/4 4*arctan(1/5) - arctan(1/239)而arctan(x)的泰勒展开是arctan(x) x - x^3/3 x^5/5 - x^7/7 ...当x 1/5或x 1/239时每一项都是有理数。我们把整个计算乘以一个足够大的缩放因子SCALE 10^(DIGITS 10)就能把小数运算变成整数运算。举个例子arctan(1/5)的第一项1/5乘以SCALE就是SCALE // 5。第二项1/(3*5^3)乘以SCALE就是SCALE // (3 * 5**3)。逐项累加直到新增项小到对DIGITS位不再有影响。实现代码如下DIGITS 10000 SCALE 10 ** (DIGITS 10) # 多出10位保护 def arctan_inverse(x_inv, scale): 计算 scale * arctan(1 / x_inv) 利用级数: arctan(1/x) 1/x - 1/(3x^3) 1/(5x^5) - ... x x_inv * x_inv # x^2 term scale // x_inv # 第一项: scale / x result term n 1 while term ! 0: term // x n 2 if n % 4 1: # n 1, 5, 9... 这一项系数为正需注意符号 pass if n % 4 3: result - term // n else: result term // n return result等等上面的符号判断有点绕。按标准写法更清晰的做法是def arctan_inverse(x_inv, scale): x x_inv * x_inv term scale // x_inv result term n 1 while term: term // x n 2 # 符号交替n3 为减n5 为加n7 为减... if n % 4 3: result - term // n else: result term // n return result然后pi_scaled 4 * (4 * arctan_inverse(5, SCALE) - arctan_inverse(239, SCALE))这里pi_scaled就是pi * 10^(DIGITS 10)的整数表示。最后要输出时只需要先除以10^10得到放大10^DIGITS倍的结果再转成字符串插入小数点。3.4 为什么这个版本更“手工”这个版本里没有Decimal.sqrt没有任意精度库全程只用了整数加、减、乘、除整除。每一步你都可以在纸上验算逻辑完全透明。对于想搞懂原理的人来说这是很好的学习材料。唯一的“性能问题”是SCALE 10^10010每一项的整数都很大Python的大整数乘法虽然做过优化但反复除法还是有点耗时。对于10000位这个量级实测在几秒到几十秒之间完全能接受。这也解释了一个常见误区“高精度计算一定要用高精度库”。其实库本质上是帮你封装了这些整数运算。当需求只是算10000位圆周率时直接用大整数模拟就够了。4. 常见问题与排查技巧实录4.1 为什么算出来的pi前面是3.14159后面全错了这是最常见的翻车现场。原因通常是两类SCALE不够大。如果你只放大10^(DIGITS)但中间运算需要再多几位来保存进位误差末尾就会错。解决办法就是留出保护位。迭代次数不够。比如Machin级数中term // x会越变越小当term变成0时说明当前项对最终结果已经没有贡献。如果你在while term ! 0的循环条件里退出了按理说是够的。但如果某个实现里用固定range(N)而N太小就会漏项。排查方法是临时打印一下循环结束时term的大小确认它已经完全为0或足够小。4.2 字符串插入小数点的位置错了10000位pi算出来是一个大整数表示的是pi * 10^(DIGITSextra)。打印前一定要把这个“额外放大倍率”还原。一个我自己经常用的做法是最终整数记为pi_big它等于floor(pi * 10^DIGITS * 10^extra)。输出时先pi_big // 10^extra得到floor(pi * 10^DIGITS)再转字符串第一位是整数部分3后面DIGITS位是小数部分。如果直接拿着最原始的整数去插小数点你会得到一堆莫名其妙的数字。4.3 算得太慢怎么办10000位其实不算大慢多半是因为实现太粗糙。两个最简单的优化手段一个是在算法层面用Chudnovsky而非Machin。Machin公式每项约贡献1.4位算10000位大约需要7000多次大整数除法虽然能跑完但明显更吃力。Chudnovsky每项约14位迭代次数少了一个数量级。另一个是在数据结构层面如果你真的在列表中逐位存数字那么每次加法都要遍历整个列表10000位还不算大但如果你哪天想算100万位这种写法会慢到你怀疑人生。改进办法就是“压位”每个列表元素存9位或18位整数而不是存0到9的一位数字这样列表长度缩短9到18倍循环次数相应减少。4.4 末位数字正确吗误差检验方法高精度计算最怕的就是“看起来对”但末尾错。这里给出三个实用的校验手段查表法从网上找已知的pi前几万位文本文件把自己算出的结果的前100位、中段100位、末尾100位拿出来比对。双算法法用Chudnovsky和Machin两个不同级数分别计算比对结果是否完全一致。因为两者逻辑完全不同同一位数同时出错的概率几乎不存在。增量法把保护位从10改成50重新跑一遍看前10000位是否保持不变。如果变了说明原结果的保护位不足需增加保护位。有一个很实用的自查点pi的小数点后前100位是公开的可以快速验证3.14159265358979323846264338327950288419716939937510 58209749445923078164062862089986280348253421170679拿自己算出来的结果从第1位到第100位对一遍如果全中说明大方向对了。5. 实操总结从10000位到更大规模的思考5.1 10000位之后还能玩什么算完10000位之后很多人会觉得“就这”其实这一步的意义在于你已经完整走了一遍高精度数值计算的流程公式选型、精度控制、整数化运算、结果校验。下一步可以尝试几件事把算法改成并行或分块挑战100万位。用C语言写一个大整数乘法把FFT快速傅里叶变换引入高精度乘法体验一下从O(n^2)到O(n log n)的飞跃。把输出结果做可视化、做艺术化展示。在这些方向上Chudnovsky算法的地位会越来越高因为极限计算场景下收敛速度决定一切。5.2 我个人的实操心得我最早自己写这个程序时犯过一个特别低级的错误忘了给中间结果保留保护位导致算出来的数最后20位全是乱的。当时我拿标准值一对比整个人都懵了排查了大半天才意识到原来不是公式问题而是精度不够。后来我养成了一个习惯就是所有高精度计算在开头先问自己三句话我的中间计算需要多少额外保护位我的最终结果会放大多少倍我的结果校验用什么标准把这三个问题想清楚基本不会翻车。最后再分享一个小技巧如果只是临时验证一下某段高精度结果可以在计算时故意把DIGITS设成100和10000各跑一次看前100位是否完全一致。如果一致说明算法稳定如果不一致说明存在精度或边界问题。这个方法几乎能定位90%以上的隐藏bug。