★ FEATURED ARTICLE 在技术社区搜pi跳出来多半是树莓派、PI控制器、pi agent这类内容真要搜“计算pi小数点后10000位”反而会掉进一堆年代久远的代码片段里有的用C语言全篇宏定义有的只贴出几千位就说“已算到一万位”。我自己动手完整做了一遍之后最大的感受是这个题目非常适合当作“任意精度计算”的入门实践它逼着你把算法收敛速度、中间截断误差、浮点数的精度天花板这些平时被框架掩盖掉的问题全部面对一遍。文章后面会给出可直接运行的两套Python实现一套用decimal模块逻辑直观一套用纯整数运算速度更好并讲清楚验证、性能、踩坑三个环节。无论是把它当面试题、项目引子还是性能基准测试这篇文章都能让你少走弯路。1. 10000位背后的真实难度一个看似简单的编程题1.1 目标不只是“算出来”而是“算对这10000位”很多人一上来就写while True: pi ...跑完把数字贴出来结果对前50位后面就乱了。这个问题的本质不是循环次数而是精度系统的搭建。我们把这个需求拆开看实际上是三个子目标得到至少10000个正确的十进制小数位而不是一个近似浮点数计算过程可以被验证别人能复现你的结果耗时可控不至于让一次计算变成等待两小时的煎熬。我把这三个子目标写进计划之后才意识到这不是一个“写个公式就完事”的题目。它横跨了数值分析选公式、估计截断误差、编程语言的数值模型float和Decimal的区别、大整数运算当数字变成10^10000级别时普通类型都失效三个层面。从规模上看10000位小数是double可表示精度的600倍左右。double在大多数语言里只有53位二进制有效数字换算成十进制大约是15到17位。也就是说用原生浮点类型你连第18位都保证不了。这个是后面所有坑的总根源。1.2 浮点数的“精度天花板”到底在哪IEEE 754规定C/C的double、Java的double、Python的float都使用64位存储其中1位符号、11位指数、52位尾数加上隐含位可以看作53位。53位二进制对应的十进制精度是log10(2^53)≈15.95所以通常说“double有16位有效数字”。用这样的类型去算π就算你键盘敲冒烟屏幕上永远只会显示3.141592653589793再往后都是噪音。要突破这个天花板只有两条路一是引入任意精度库比如GMP、MPFR、Java的BigDecimal、Python的decimal二是在整数空间里做运算把小数部分放大到10的N次幂全程用整数加减乘除最后再把小数点插回去。后面的整数版本用的就是第二条路。2. 算法选型哪些公式能撑起一万位圆周率公式在数学史上非常多但真正适合编程计算的就那么几类。我筛选时首先放弃的不是“错”的公式而是“收敛太慢导致物理意义上不可能”的公式。2.1 蒙特卡洛与莱布尼茨级数入门可以上万位不行蒙特卡洛法往正方形里随机撒点靠面积比估计π。随机采样的误差收敛速度是O(1/√N)也就是说要把误差压到10^(-10)需要10^20次采样。放到10000位需要的采样次数是10^20000级别宇宙毁灭都算不完。莱布尼茨级数π4(1-1/31/5-1/7...)就更夸张了。它是一个交替调和级数误差衰减速度是1/(2k1)。每算一项小数点后的有效位数增加约0.3位。想靠它算到10000位需要约3×10^10000项同样不可能。这两个例子说明一个关键判断标准当你要冲击极高精度时级数的项与精度的关系必须是对数级的或者至少是幂级数中收敛极快的否则就是死路。2.2 马青公式中等精度绕不开的经典马青公式Machin formula是1706年发现的π 16 arctan(1/5) - 4 arctan(1/239)把arctan展开成泰勒级数arctan(x) x - x^3/3 x^5/5 - ...代入x1/5和x1/239之后每一项的大小分别按1/25和1/57121的比例衰减。1/25的log10是-1.39794也就是说arctan(1/5)的级数每迭代一项大约能多1.4位小数而arctan(1/239)的每项衰减是4.756位。要达到10000位精度考虑上截断余量arctan(1/5)需要大约7200项arctan(1/239)需要大约2200项。这个计算量非常温和现代CPU毫秒级就能跑完。这也是为什么马青公式是“万位级精度”最实用的选择。2.3 再看一眼Chudnovsky精度更高但复杂度也更高Chudnovsky算法是1989年提出的公式长这样π 426880 √10005 / Σ_{k0}^∞ ( (6k)! (13591409 545140134k) ) / ( (3k)! (k!)^3 (-262537412640768000)^k )它的优点非常吓人每一项贡献约14.18位十进制有效数字。算10000位只需要约710项算1亿位也就600多万项。但代价是每一项都要做超大整数的阶乘、乘方和除法还涉及高精度的平方根计算。要用好它通常需要配合二进制分割binary splitting技术代码复杂度直接上一个台阶。对于10000位这个精度马青公式和Chudnovsky差距并不大。我的建议是如果你把这次任务当作算法学习马青公式足够如果你打算以后冲击百万位、千万位那直接学Chudnovsky更值。2.4 我最终选型先马青再用整数优化我的最终方案分成两步先用马青公式的Decimal版本把逻辑跑通验证前几百位正确再切换成整数运算版本把速度提上去。这样的好处是两个实现互为参照算出来的结果还可以交叉验证一旦有一个出问题立刻能发现。3. 从公式到代码两种可落地的实现方案我用的语言是Python 3。先声明一点Python内置的float完全不参与这次计算核心是decimal模块和大整数。3.1 Decimal版本最容易读懂的实现Python的decimal模块提供了任意精度的十进制浮点数核心是把精度上下文getcontext().prec设成目标位数。下面是完整实现from decimal import Decimal, getcontext def arctan_inv_decimal(x, n): 计算 arctan(x) 的泰勒级数x 必须是 Decimal total Decimal(0) term x xx x * x sign 1 for k in range(1, 2 * n, 2): total sign * term / k term * xx sign -sign return total def calc_pi_decimal(ndigits10000): # 留出20位余量避免中间舍入污染最后一位 getcontext().prec ndigits 20 # 迭代次数粗略估算arctan(1/5) 需要约 ndigits/1.397 项 n int(ndigits / 1.3) 300 a arctan_inv_decimal(Decimal(1) / Decimal(5), n) b arctan_inv_decimal(Decimal(1) / Decimal(239), n) pi 16 * a - 4 * b return str(pi)[:ndigits 2] if __name__ __main__: print(calc_pi_decimal(10000))注意两个关键点getcontext().prec必须在做除法之前设置。如果在默认精度28下先算Decimal(1) / Decimal(239)得到的是一个只有28位有效数字的数后续无论怎么加精度误差已经埋进去了。迭代次数n不用算得特别精确取大一点不亏最多多跑几千次循环但取小了最后若干位就是错的。3.2 整数运算版本更快、更可控Decimal版本容易理解但每次循环都做Decimal除法本质上是模拟十进制浮点运算开销不低。更贴近底层、也更快的方式是把整个结果放大10^prec倍用纯整数来算泰勒级数。思路是这样的arctan(1/d)的第k项是(-1)^(k) / ((2k1) * d^(2k1))我先把分子固定为10^prec用一个整数term表示当前项放大后的值def arctan_int(den, prec): 计算 arctan(1/den) * 10^prec 的整数近似值 total 0 term 10 ** prec // den k 1 sign 1 den2 den * den while term: total sign * (term // k) term // den2 k 2 sign -sign return total def calc_pi_int(ndigits10000): prec ndigits 20 # 余量留大一点更稳 a arctan_int(5, prec) b arctan_int(239, prec) pi_int 16 * a - 4 * b return pi_int这里term // den2的作用是让当前项从1/5^(2k-1)过渡到1/5^(2k1)每一步只需要一次大整数除法。整个过程中所有的数都是整数不存在浮点舍入误差只来自每一次整除的向下取整。由于我留了20位余量向下取整带来的损失会被控制在最后十几位以内不会污染前10000位。3.3 输出格式与运行效果整数版本算出来的pi_int是一个大约有10020位数字的大整数第一位是3后面跟着10019位小数部分。输出时只需要把它转成字符串然后在第一位后面插入小数点def pi_to_string(pi_int, ndigits): s str(pi_int) # 防止某些极端情况下整数位数不够先补零 if len(s) ndigits 1: s s.zfill(ndigits 1) return s[0] . s[1:ndigits 1] pi_int calc_pi_int(10000) print(pi_to_string(pi_int, 10000))在我的笔记本上跑一遍前几行输出是3.14159265358979323846264338327950288419716939937510 58209749445923078164062862089986280348253421170679 ...第一眼看到这个结果我就知道整个流程跑通了。但“看到了π”和“确认这一万位全对”是两码事我单独把验证环节拎出来说。4. 验证结果算出来的10000位怎么保证没错很多人算出结果就结束了但如果你真要把这个结果用于基准测试、算法对比或者教学演示一定要做验证。这里分享几种我用下来觉得靠谱的方式。4.1 前缀对比前100位一眼定胜负π的前100位是公开常数随手可查3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679我在代码里固定存了一段前缀字符串算完后直接startswith检查。这一步能过滤掉90%的明显错误公式抄错、泰勒展开符号错、小数点位置错基本都逃不过这双火眼金睛。4.2 交叉验证用两个独立实现互算我前面特意保留了两套实现Decimal版和整数版它们各有各的舍入来源。用同一个马青公式分别算10000位再把字符串做一次全量对比如果完全一致那基本可以判定正确。交叉验证里有个容易被忽略的细节两套实现要尽量独立不要复制同一份代码。我的Decimal版和整数版从数据结构、循环方式到误差来源都不一样交叉验证才有意义。如果你只是改改变量名验证就是自欺欺人。4.3 分段切片核对与哈希校验前缀对比只能证明开头对交叉验证能证明两套代码一致但还不能证明“两套代码一起错了”这种极端情况。为了彻底打消疑虑可以引入第三方结果。方法很简单找一个与你的代码完全无关的高精度计算工具比如gmpy2.const_pi()、mpmath的mp.dps10000; mpmath.pi或者从OEIS、可信的开源仓库下载标准π文本文件然后在随机位置分段切片做对比。我习惯的做法是抽查三处第1000位附近取第990到1010位第5000位附近取第4990到5010位第9990位附近取第9970到10000位。如果这三段都能对上那基本可以确认算到了第10000位。更进一步把整个10000位字符串做一次SHA256哈希与官方文本的哈希比对一旦对上连“中间某处错一位”的可能也被排除。在Python里做这个只是几行代码的事。5. 性能实测与优化路径5.1 三个精度档位的耗时对比我在自己的笔记本Intel i58GB内存Python 3.11上分别跑了1000位、10000位、100000位耗时量级大致如下目标位数Decimal版耗时整数版耗时1,000约0.05秒约0.02秒10,000约0.9秒约0.2秒100,000约60秒约7秒这个数据不是精确基准不同机器差异很大但量级关系是稳定的Decimal版在10万位时有明显吃力感整数版快了近一个数量级却也开始逼近秒级。5.2 性能瓶颈到底在哪里马青公式的计算量由两部分组成第一是级数项数前面算过10000位需要约72002200项100000位就需要约7200022000项项数和精度成正比。第二是每一项操作的大数规模。整数版里的term有10^prec量级也就是10万位时需要处理一个十万位的整数每做一次整除开销跟大数的字节长度成正比。项数乘上每次操作的大数长度总复杂度大致是O(n^2)。这就是为什么1000位时感觉不到时间100000位时明显卡顿。Python的大整数虽然有C语言底层优化但O(n^2)的曲线摆在那里位数每翻10倍时间要翻约100倍。5.3 更进一步从O(n^2)往O(n log n)走如果只是算到10000位优化空间已经不大。但如果想继续冲更高精度可以考虑三条路用gmpy2库替换Python原生int和decimal。gmpy2.mpz底层是GMP大整数乘除法比Python原生快数倍到数十倍同样的马青公式代码改成gmpy2之后10万位能压进1秒以内。换Chudnovsky公式加二进制分割。二进制分割能把阶乘和级数求和变成分治形式总复杂度降到接近O(n log n)这是目前百万位以上的主流做法。如果只是临时验证直接用gmpy2.const_pi(prec)它会调用MPFR把π算到任意精度一行代码速度还快得离谱。不过对于10000位这个目标我的结论是没必要为了性能引入复杂方案马青整数运算已经是性价比最高的组合。6. 踩坑记录几个容易让结果悄悄出错的地方最后这部分是这次实操里最值钱的经验。下面每个坑我都实际踩过或者看着它们让结果悄悄出错。6.1 Decimal精度设置在“操作之前”而不是“表达式之前”Python的decimal精度上下文是全局状态不是表达式属性。最容易犯的错误是getcontext().prec 28 x Decimal(1) / Decimal(239) # 这里的除法已经在28位精度下算完了 getcontext().prec 10050 # 改晚了 pi 16 * x - ...你以为后面把精度调到10050x已经是一个28位精度的数后续计算结果的有效位数最多28位。正确做法是先把getcontext().prec设为目标精度再执行任何除法、开方等会产生舍入的操作。6.2 整数版本的整除截断误差与余量设计整数版本每一步term // den2都会丢掉一点余数项数越多向下取整的累计误差越大。我最初用prec ndigits去跑结果第9990位开始就和参考值对不上。后来把余量从10加到20尾端才稳定下来。所以整数版本的余量不能省。prec ndigits 20是我实测够用的值但如果你的机器环境不同建议算完后抽查尾部100位。6.3 迭代次数不足比超跑更可怕Decimal版本里我把迭代次数设成int(ndigits / 1.3) 300这个经验值够用。但如果你图省力写nndigits也能过写nndigits//2就会在很靠后的位置出现错误。重要的是理解估算逻辑arctan(1/5)每项增加约1.4位arctan(1/239)每项增加约4.8位按精度需求反推项数再留出10%左右的余量就不会踩坑。6.4 字符串输出时的长度与补零整数版本里π乘以10^prec后整数部分有prec1位字符串长度足够一般不用补零。但如果你把prec设成ndigits20str长度是ndigits21切片时[1:ndigits1]会正确拿到10000位。如果你在别的公式里遇到首项特别大或特别小的情况建议还是加一句zfill兜底。这种细节在现场跑数据时最磨人宁可多写一行防御代码。我个人做这类项目习惯先写一个简单的验证函数把前缀、中间段、尾部段三处断言写进去每次改完代码立刻全量自检。这样即使后面迭代了很多版本也不会在某个深夜把一段错误的结果当成“正确的一万位”发布出去。这个计算任务表面上是玩数字实际上把数值稳定性、大整数运算、算法复杂度分析全练了一遍。如果你也想动手试试建议从马青公式的整数版本开始一步步把代码写出来再亲手踩一遍精度余量的坑。等你能稳定输出并验证10000位时后面再接触Chudnovsky、二进制分割这些高阶技巧会轻松很多。
阅读完成 · 觉得有帮助？
