在技术社区搜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/3+1/5-1/7+...)就更夸张了。它是一个交替调和级数，误差衰减速度是1/(2k+1)。每算一项，小数点后的有效位数增加约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 - ...，代入x=1/5和x=1/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 / Σ_{k=0}^∞ ( (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(ndigits=10000): # 留出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) / ((2k+1) * d^(2k+1))，我先把分子固定为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(ndigits=10000): 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^(2k+1)，每一步只需要一次大整数除法。整个过程中所有的数都是整数，不存在浮点舍入，误差只来自每一次整除的向下取整。由于我留了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.dps=10000; mpmath.pi，或者从OEIS、可信的开源仓库下载标准π文本文件，然后在随机位置分段切片做对比。
我习惯的做法是抽查三处：
- 第1000位附近，取第990到1010位；
- 第5000位附近，取第4990到5010位；
- 第9990位附近，取第9970到10000位。
如果这三段都能对上，那基本可以确认算到了第10000位。更进一步，把整个10000位字符串做一次SHA256哈希，与官方文本的哈希比对，一旦对上，连“中间某处错一位”的可能也被排除。在Python里做这个只是几行代码的事。
5. 性能实测与优化路径
5.1 三个精度档位的耗时对比
我在自己的笔记本（Intel i5，8GB内存，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位需要约7200+2200项，100000位就需要约72000+22000项，项数和精度成正比。
第二是每一项操作的大数规模。整数版里的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 - ...你以为后面把精度调到10050，x已经是一个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，这个经验值够用。但如果你图省力写n=ndigits，也能过；写n=ndigits//2，就会在很靠后的位置出现错误。重要的是理解估算逻辑：arctan(1/5)每项增加约1.4位，arctan(1/239)每项增加约4.8位，按精度需求反推项数，再留出10%左右的余量，就不会踩坑。
6.4 字符串输出时的长度与补零
整数版本里π乘以10^prec后，整数部分有prec+1位，字符串长度足够，一般不用补零。但如果你把prec设成ndigits+20，str长度是ndigits+21，切片时[1:ndigits+1]会正确拿到10000位。如果你在别的公式里遇到首项特别大或特别小的情况，建议还是加一句zfill兜底。这种细节在现场跑数据时最磨人，宁可多写一行防御代码。
我个人做这类项目，习惯先写一个简单的验证函数，把前缀、中间段、尾部段三处断言写进去，每次改完代码立刻全量自检。这样即使后面迭代了很多版本，也不会在某个深夜把一段错误的结果当成“正确的一万位”发布出去。
这个计算任务，表面上是玩数字，实际上把数值稳定性、大整数运算、算法复杂度分析全练了一遍。如果你也想动手试试，建议从马青公式的整数版本开始，一步步把代码写出来，再亲手踩一遍精度余量的坑。等你能稳定输出并验证10000位时，后面再接触Chudnovsky、二进制分割这些高阶技巧，会轻松很多。
