用Python实现五种圆周率算法：从蒙特卡洛到Chudnovsky的精度与性能对比-北京尧图网络

前一阵子有个朋友在群里发需求“谁能用代码写个圆周率计算”结果底下清一色回复“import math然后print(math.pi)”。说真的这种回复也没错Python 内置的math.pi已经精确到双精度浮点能表达的极限了。但“圆周率”这三个字放在程序员面前从来就不只是“拿到结果”那么简单。pi 这个只有两个字母的标题背后其实藏着一整条数值算法的主线从随机模拟到无穷级数从经典公式到现代高性能算法每一步都在教你同一件事——同一道数学题工程实现上的取舍能差出多少个数量级。这篇文章我想用 Python 把五种计算 pi 的经典算法全部跑一遍蒙特卡洛、莱布尼茨级数、马青公式、Chudnovsky 算法还有 BBP 公式。每写一个我都会给出代码、实测精度、运行耗时和它背后的数学原理最后再横向对比一次。适合的人群很明确想理解数值分析基本概念的 Python 开发者、正在准备算法相关面试的人以及想拿“计算 pi”当练手项目但不知道从哪个算法切入的朋友。1. 为什么值得亲手算一次圆周率从“调函数”到“理解误差”你当然可以永远只调math.pi但一旦开始自己实现你碰到的第一个问题就是float 到底是什么精度Python 里的 float 对应 C 语言的 double总共 53 位二进制尾数折算下来大约是 15 到 16 位十进制有效数字。也就是说math.pi本身就只有小数点后十几位是可信的再往后打印出来的那一串数字其实是二进制到十进制的转换产物并没有真实的数学意义。为什么要搞清楚这件事因为“算 pi”这个需求在不同场景下的答案完全不同如果你只想在 UI 上画个圆math.pi绰绰有余如果你想演示“算法如何收敛到正确答案”你需要的是能看出误差变化的算法如果你想挑战高精度就需要改用十进制浮点甚至任意精度整数运算这时候不是“调库”能解决的。1.1 五种算法的定位差异我先按自己的理解把五种算法排了个序方便你建立整体印象算法核心思路收敛速度实现难度典型用途蒙特卡洛随机落点统计1/sqrt(N)极慢低教学演示、并行计算实验莱布尼茨级数反正切展开取 x11/N很慢极低理解级数收敛概念马青公式多个反正切组合每项约 1.4 位十进制低双精度范围内的经典选择Chudnovsky模群相关的快速级数每项约 14 位十进制中打破世界纪录的主力BBP 公式计算任意十六进制位每项 1/16^k中无需前置位直接提取某一位这张表其实反映了数值计算的核心矛盾速度不是唯一的评价指标。蒙特卡洛虽然精度差但它教会你“误差来自随机性”莱布尼茨虽然慢得离谱但它是最容易理解的级数Chudnovsky 快得惊人但代码复杂度也直线上升。下面我按“从易到难”的顺序逐个拆开讲。1.2 动手前的环境准备本文所有示例代码都用纯 Python 标准库完成不需要安装任何第三方包。唯一建议你准备好的是 Jupyter Notebook 或任意一个 .py 脚本环境因为后面有几个实验需要反复修改迭代次数和输出位数。我用的 Python 版本是 3.11理论上 3.8 以上都能正常运行。2. 蒙特卡洛方法随机落点带来的误差课蒙特卡洛大概是所有人最熟悉的 pi 算法在一个边长为 2 的正方形里画一个半径为 1 的内切圆然后随机往正方形里撒点统计落在圆内的比例。因为圆的面积是 π正方形的面积是 4所以比例 p 应该等于 π/4反推一下就是 π 4p。2.1 基础实现与关键代码import random import math import time def estimate_pi_monte_carlo(num_points: int, seed: int 42) - float: random.seed(seed) inside 0 for _ in range(num_points): x random.uniform(-1.0, 1.0) y random.uniform(-1.0, 1.0) if x * x y * y 1.0: inside 1 return 4.0 * inside / num_points start time.perf_counter() pi_est estimate_pi_monte_carlo(100_000) elapsed time.perf_counter() - start print(fpi_est {pi_est:.6f}, error {abs(pi_est - math.pi):.6f}, time {elapsed:.4f}s)我固定了随机种子跑 10 万次的结果大致在 3.1446 附近误差在 0.003 这个数量级耗时大约 0.03 秒。但你换一个 seed结果可能变成 3.1398也可能变成 3.1452。这就是蒙特卡洛方法最核心的特征它不是确定性算法每次运行结果都会波动。2.2 为什么误差只能按 1/sqrt(N) 的速度下降很多刚开始学的人会有一个直觉点数多 100 倍误差就应该缩小 100 倍。但实际上蒙特卡洛误差的收敛速度是 O(1/sqrt(N))也就是点数增加到原来的 100 倍误差才缩小到原来的 1/10。这个结论的来源是二项分布的方差。每个点是否落在圆内是一个伯努利随机变量落在圆内的概率是 π/4所以单次实验的方差是 p(1-p)。对 N 次独立实验取平均均值估计的标准差是 sqrt(p(1-p)/N)。把这个标准差乘以 4大约就是 π 估计值的标准误。代入 p ≈ 0.7854算出来的结果是标准误 ≈ 4 * sqrt(0.7854 * 0.2146 / N) ≈ 1.642 / sqrt(N)N100,000 时标准误大约是 0.0052。也就是说一次运行结果和真值差 0.003完全在正常波动范围内。蒙特卡洛方法想要小数点后 6 位精度理论上需要约 2.7 万亿个点——在普通笔记本上基本不现实。2.3 蒙特卡洛的实用场景别拿它拼精度我实际用下来觉得蒙特卡洛最大的价值不是算 pi而是让人第一次直观地感受到“误差传播”是怎么回事。你可以做一个简单的小实验固定 seed分别用 10、100、1000、10000、100000 个点跑一次把结果画成曲线。曲线下降的速度会让你对“随机模拟的成本”产生非常真实的体感。另外如果你在做并行计算入门实验蒙特卡洛也是一个很好的例子因为它天然可以拆成多个独立任务每个任务各算一部分落点最后汇总。我在多进程版本里跑过 8 个进程每个 100 万点总耗时反而比单进程快很多但精度依然只提高了不到一个数量级——这就是随机方法的宿命。3. 莱布尼茨级数收敛慢不是缺点是教材莱布尼茨公式长这样π/4 1/1 - 1/3 1/5 - 1/7 1/9 - ...所以 π 4 * (1 - 1/3 1/5 - 1/7 ...)。它的推导过程很优雅反正切函数 arctan(x) 的泰勒展开在 x1 处恰好就是 1 - 1/3 1/5 - 1/7 ...而 arctan(1) 正好是 π/4。3.1 代码实现与一个容易被忽略的优化def pi_leibniz(num_terms: int) - float: total 0.0 # 倒序累加先加小项后加大项减少大数吃小数的误差 for k in range(num_terms - 1, -1, -1): term 1.0 / (2 * k 1) total term if k % 2 0 else -term return 4.0 * total start time.perf_counter() pi_leibniz_est pi_leibniz(1_000_000) elapsed time.perf_counter() - start print(fpi_est {pi_leibniz_est:.9f}, error {abs(pi_leibniz_est - math.pi):.9f}, time {elapsed:.4f}s)这里有个细节我故意让循环从大到小累加。原因很简单浮点数相加时如果一个大数和很多小数加在一起小数的尾数可能直接丢失。倒序累加可以尽量保住小数的贡献。虽然在这个示例里影响不一定肉眼可见但这是数值计算里的一个通用习惯后面 Chudnovsky 算法里的 decimal 大数累加也一样适用。实测跑 100 万个项我得到的结果大约是 3.14159065误差约 1.6e-6耗时 0.4 秒左右。如果你只叠加 10 万项误差只到小数点后第 5 位。3.2 用泰勒余项解释“为何如此之慢”交错级数的误差有个经典结论如果截断在第 n 项截断误差的绝对值小于等于被舍弃的第一项。对莱布尼茨级数来说第 n 项大约是 1/(2n1)。要判断这个过程能算到多少位只需要解一个不等式4 / (2n1) 10^(-d)也就是说每多要一位十进制精度n 至少变成原来的 10 倍。想要小数点后 1000 位理论上要循环 10^1000 次这是个天文数字。所以莱布尼茨级数最大的价值就是让你对“线性收敛的极限”产生敬畏。它教会我们一个道理一个公式在数学上是对的不代表它在工程上是可用的。3.3 别把“欧拉公式”和“莱布尼茨级数”混为一谈我在评论区见过不少人说“明明欧拉公式算 pi 很快你怎么说莱布尼茨慢”。这里要澄清一下欧拉自己确实改进过反正切类公式但通常说的“欧拉公式”是指 π/4 arctan(1/2) arctan(1/3)这本质上已经属于“多个反正切相加”的思路和马青公式是同一个家族。它的收敛速度确实比莱布尼茨快很多但那是公式设计的功劳不是“欧拉”这个名字带来的魔力。学习时先把最基本的莱布尼茨吃透再去玩组合公式脉络会清晰很多。4. 马青公式双精度时代最优雅的折中马青公式长这样π/4 4 * arctan(1/5) - arctan(1/239)也就是说π 16 * arctan(1/5) - 4 * arctan(1/239)为什么这个公式快因为 arctan(x) 的泰勒展开是 x - x^3/3 x^5/5 - ...而这里两个自变量的绝对值分别是 1/5 和 1/239远小于 1。级数的每一项会按 x^2 的比例衰减1/5 的平方是 1/251/239 的平方约等于 1/57121所以每一项都能带来大约两位十进制精度。相比之下莱布尼茨级数的 x1衰减比例是 1/1本质上不衰减当然慢得离谱。4.1 用 float 实现逼近 math.pi 的极限def arctan_taylor(x: float, terms: int 10) - float: total 0.0 x_power x for k in range(terms): term x_power / (2 * k 1) total term if k % 2 0 else -term x_power * x * x # x^(2k3) return total def pi_machin(terms: int 10) - float: return 16.0 * arctan_taylor(1.0 / 5.0, terms) - 4.0 * arctan_taylor(1.0 / 239.0, terms) pi_machin_est pi_machin(10) print(fpi_est {pi_machin_est:.16f}) print(fmath.pi {math.pi:.16f}) print(ferror {abs(pi_machin_est - math.pi):.2e})我实际运行的结果是当terms10时pi_machin_est与math.pi在小数点后 15 位都完全一致误差小于 1e-15。再往下比较已经没有意义了因为第 16 位开始 float 本身就开始出现舍入误差。4.2 为什么马青是“双精度终点”正如前面说的float 只有约 15 到 16 位有效数字马青公式在 terms10 时已经把这个精度榨干了。所以如果只是想在双精度范围内快速拿一个准值马青公式是五个算法里性价比最高的代码短、原理清楚、速度极快terms10 时耗时不到 0.0001 秒。但如果你想要更多位数比如 50 位、100 位就必须换工具了。下一节的 Chudnovsky 算法才是真正用于高精度计算的选手。5. Chudnovsky 算法与 BBP 公式高精度与随机访问的两条路5.1 Chudnovsky 算法一节更比六节强Chudnovsky 公式是当前计算机计算 π 的世界纪录最常用的公式。它的标准形式让人一看就头大1/π 12 * Σ_{k0}^{∞} (-1)^k * (6k)! * (13591409 545140134k) / ((3k)! * (k!)^3 * 640320^(3k 3/2))但这个公式每一次迭代能产出大约 14 位十进制精度所以只需要迭代几次就能轻松超过 float 的精度上限。实现时需要注意一点640320^(3k 3/2)里有 3/2 次方如果直接用普通 float 参与运算高位小数会丢失。所以我改成如下等价形式把所有大整数保留在分子分母里只在最后一步做一次开方π 426880 * sqrt(10005) / S S Σ_{k0}^{∞} (-1)^k * (6k)! * (13591409 545140134k) / ((3k)! * (k!)^3 * 640320^(3k))这样分母里的640320^(3k)是整数可以精确计算最后的sqrt(10005)用 Decimal 处理即可。from decimal import Decimal, getcontext import math def pi_chudnovsky(iterations: int, precision: int 60) - Decimal: getcontext().prec precision 10 total Decimal(0) for k in range(iterations): numerator math.factorial(6 * k) * (13591409 545140134 * k) denominator math.factorial(3 * k) * (math.factorial(k) ** 3) * (640320 ** (3 * k)) term Decimal(numerator) / Decimal(denominator) total term if k % 2 0 else -term pi Decimal(426880) * Decimal(10005).sqrt() / total return pi pi_high pi_chudnovsky(5, precision70) print(pi_high)我用 5 次迭代、精度设置为 70 位时得到的结果是3.1415926535897932384626433832795028841971693993751058209749445923078164这串数字前 70 位和公开的 π 值完全一致。注意5 次迭代产出了差不多 70 位换算一下相当于每次迭代稳定提供 14 位这个效率是莱布尼茨那种线性收敛完全无法想象的。5.2 BBP 公式不关心前面的位直接抓某一位BBP 公式是 1995 年由 Bailey、Borwein 和 Plouffe 提出的它长这样π Σ_{k0}^{∞} 1/16^k * (4/(8k1) - 2/(8k4) - 1/(8k5) - 1/(8k6))这个公式最大的特别之处在于它可以直接计算 π 的十六进制小数的某一位不需要知道前面的位。原理是在求和时如果只关心第 n 位小数那么对于 k ≤ n 的那些项可以先做模运算把整数部分扔掉只保留小数部分对于 k n 的项则正常累积一个很小的量。这样就能在 O(n log n) 的时间复杂度内算出第 n 位。下面是我用来展示 BBP 公式前 20 位十六进制小数的简化代码适合理解原理from decimal import Decimal, getcontext getcontext().prec 60 pi_val Decimal(0) for k in range(30): base Decimal(16) ** k term (Decimal(4) / (8 * k 1) - Decimal(2) / (8 * k 4) - Decimal(1) / (8 * k 5) - Decimal(1) / (8 * k 6)) / base pi_val term frac pi_val - int(pi_val) hex_digits 0123456789ABCDEF hex_frac for _ in range(24): frac * 16 digit int(frac) hex_frac hex_digits[digit] frac - digit print(fpi integer part: 3) print(fpi hex fraction: {hex_frac})输出结果应该是pi integer part: 3 pi hex fraction: 243F6A8885A308D313198A2E对照公开的十六进制 π 值3.243F6A8885A308D313198A2E...前 24 位完全匹配。这种“跳着读位”的能力是 BBP 公式独有的其他所有级数都不具备。虽然它算的是十六进制位而不是十进制位工程上并不常用但在密码学、随机数检验等领域它是一把很好的技术锤子。5.3 两个现代算法的选型建议如果只是个人做着玩想突破 float 精度我建议直接用 Chudnovsky因为它的 decimal 实现思路直观迭代次数少容易验证。如果你对“随机访问某一位”这个能力感兴趣或者在做一个需要生成 pi 的不同区段做校验的项目BBP 才是合适的选择。日常业务代码里这两者都用不上但作为算法素养的训练它们比背一堆 API 有价值得多。6. 实测经验与踩坑记录误差、性能与停止条件这部分写点我实际跑实验时遇到的坑以及我最后留下的工程习惯。6.1 验证结果不能“看着像就行”我最开始跑 Chudnovsky 时发现输出从第 50 位开始和网上的 π 常量串对不上第一反应以为公式写错了。后来查了半天才发现问题出在 Decimal 的精度设置上——getcontext().prec不只是控制输出的小数位数它会影响每一次乘法和除法的中间精度。如果你设置 precision50但分母里那些大整数转换成 Decimal 的过程中已经损失了精度那结果自然从某一位开始漂移。我的做法是验证用的精度一定要比目标位数多留 10 位。比如想要 60 位结果getcontext().prec就设成 70。同理如果一项计算里同时有math.factorial这种整数和 Decimal 除法要确保整数到 Decimal 的转换发生在精度足够的环境下。6.2 为什么“倒序累加”是数值计算的基本功前面莱布尼茨那里我提过倒序累加在 Chudnovsky 的实现里同样适用。虽然 Decimal 的尾数比 float 大得多但只要精度有限大数吃小数的现象就依然会存在。特别是当你累计到后面几十项时每一项的量级可能比前几项小几十个数量级顺序累加会白白丢掉很多有效数字。我的习惯是无论用什么精度只要是在做高精度求和一律先算完所有项再排序累加或者干脆倒序循环。6.3 性能测量与“什么时候该停手”我在表格里汇总一下我实际跑出来的平均耗时不同机器会略有差异算法迭代规模大约耗时达到精度蒙特卡洛100,000 点0.03s约 0.003 绝对误差莱布尼茨1,000,000 项0.40s约 1.6e-6马青公式10 项0.0001s约 1e-15Chudnovsky5 项0.01s约 70 位十进制BBP 公式30 项0.02s前 24 位十六进制性能上最值得注意的是别为了追求更多位数无限提高迭代数量。Chudnovsky 我只跑了 5 项就得到了 70 位10 项能到 140 位但其中 100 位之后已经是在验证“输出位数是否正确”而不是在实际解决问题。如果你不是在做破纪录项目一个朴素的原则是目标位数乘以 0.8再留出验证余量就足够停手了。6.4 把 float、Decimal 和分数三种数值类型分清楚很多初学者会在一个函数里混用 float 和 Decimal这是我见过最多的问题。写这种实验代码之前先想清楚每一层的数值类型float适合蒙特卡洛、莱布尼茨、马青因为它们的目标本来就在双精度范围Decimal适合 Chudnovsky因为要输出几十上百位小数分数 Fraction适合做严格推导和验证但计算极慢不适合跑大迭代。我自己的流程是先用 Fraction 做小规模验证确认公式本身的整数运算没有逻辑错误再切换到 Decimal 跑最终版本。这样能避免公式写错之后在高精度里浪费排查时间。最后再分享一个小技巧如果你也想写一篇类似的实验博客或者技术分享我有一个经过多次验证的展示顺序先丢出蒙特卡洛的误差曲线让大家直观看到“随机方法的天花板”再用莱布尼茨展示“简单级数的慢性子”接着用马青公式让大家意识到“双精度其实很够用”最后用 Chudnovsky 和 BBP 收尾把视野拉到高精度和十六进制世界。这样安排的节奏会让读者感受到每一种算法都不是孤立的它们共同描绘了“算 pi 这件事的完整技术地图”。这套代码跑完我自己最大的体会是工程上选算法本质上是在“精度、速度、代码复杂度”三个维度里找平衡点。下一次当有人再和你说“直接用 math.pi 不就行了”你就可以把这篇文章甩过去然后问他那你觉得 3.141592653589793 后面的第 60 位是什么

用Python实现五种圆周率算法：从蒙特卡洛到Chudnovsky的精度与性能对比

相关推荐

VPX62电源标准中文解读：拓扑、键控与合规判读

NGT维修手册实战：60V双电池断电、拆解与故障排查指南

小程序支付取消后再支付的高可用方案设计

更多资讯

THK精密运动部件选型与失效分析全指南

智能体任务环境建模与架构设计实战指南

智能体系统架构与核心运行机制解析

ASP.NET在线考试系统源码解析：从环境搭建到二次开发避坑指南

t3code：基于T3 Stack的全栈类型安全脚手架实践与避坑指南

t3code 类型优先全栈实践：TypeScript 类型安全开发指南

今日更新

微信个人号API开发案例：图片/语音/视频消息怎么下载解析存到本地？

Python股票自动交易系统源码解析：从行情获取到风控的完整链路

视酷酷信IM源码部署指南：Java 17+MySQL 8.0.33实战配置

本周精选

防幻觉智能体架构：电源硬件设计中的大模型编排与工具化实践

OpenGL面试核心考点：渲染管线、矩阵变换与性能排查全解析

全国计算机二级C语言真题PDF高效刷题指南

本月热门

我发现了一个新思路：用 Remotion + Claude Code 像写代码一样自动化生成短视频

Windows下 Codex 中 Chrome 和 Computer Use 插件不可用问题排查及解决参考方式：TaoToken 统一 Key 配置与验证

2026 大模型集体涨价：用 Python 做企业 Token 成本测算与选型避坑（附配置）

读完文章，想聊聊您的网站？