Python实战：用割圆法、蒙特卡洛等5种算法手算圆周率（附完整代码与避坑指南）

Python实战：5种圆周率算法从原理到调优全解析

记得第一次接触圆周率计算时，我被各种算法搞得晕头转向——为什么简单的π能有这么多算法？每种算法输出的精度为何差异这么大？直到把割圆法的多边形画到纸上，用蒙特卡洛法手动撒豆子模拟后，才真正理解算法背后的数学之美。本文将带你用Python实现五种经典算法，并深入分析它们的适用场景和优化技巧。

1. 算法原理与核心思想对比

1.1 割圆法：几何逼近的鼻祖

刘徽在公元263年提出的割圆术，核心思想非常简单：用正多边形逼近圆形。当边数足够多时，多边形周长就越接近圆周长。其数学表达为：

π ≈ (边长 × 边数) / 直径

实现时需要注意两个关键点：

• 边长迭代公式的推导准确性
• 初始边数设置（通常从正六边形开始）

收敛速度：每次分割边数翻倍，精度提升约4倍（线性收敛）

1.2 无穷级数法：微积分的礼物

莱布尼茨级数是最简单的无穷级数表示：

π/4 = 1 - 1/3 + 1/5 - 1/7 + ...

但实际应用中更常用马青公式：

π = 16arctan(1/5) - 4arctan(1/239)

注意事项：

• 需要设置合理的停止阈值
• 交替级数的收敛速度较慢（300项才精确到3.14）

1.3 蒙特卡洛法：概率的力量

这个方法的奇妙之处在于用随机性解决确定性问题。在单位正方形内随机撒点，统计落在1/4圆内的比例：

π ≈ 4 × (圆内点数/总点数)

关键参数：随机数种子影响结果可复现性，点数越多精度越高但计算量越大

2. 代码实现与性能优化

2.1 割圆法实现技巧

原始代码可以优化为：

import math def archimedes_pi(iterations): sides = 6 length = 1.0 for _ in range(iterations): length = math.sqrt((1 - math.sqrt(1 - (length/2)**2))**2 + (length/2)**2) sides *= 2 return length * sides / 2 # 测试10次分割（正6×2^10边形） print(archimedes_pi(10)) # 输出3.141592345570118

常见问题排查：

1. 结果不收敛？检查边长迭代公式是否正确
2. 精度不足？增加分割次数（但注意浮点数精度限制）

2.2 蒙特卡洛法的并行优化

原始串行实现效率低下，可用multiprocessing加速：

from multiprocessing import Pool import random def monte_carlo_chunk(size): count = 0 for _ in range(size): x, y = random.random(), random.random() if x**2 + y**2 <= 1: count += 1 return count def parallel_monte_carlo(total_points, workers=4): chunk_size = total_points // workers with Pool(workers) as p: results = p.map(monte_carlo_chunk, [chunk_size]*workers) return 4 * sum(results) / total_points print(parallel_monte_carlo(10_000_000)) # 千万级点数只需几秒

性能对比：

2.3 拉马努金公式的数值优化

原始实现直接计算阶乘会导致数值溢出，应使用对数变换：

import math def ramanujan_pi(n): total = 0 ln_constant = math.log(2*math.sqrt(2)/9801) for k in range(n): numerator = math.lgamma(4*k + 1) + math.log(1103 + 26390*k) denominator = 4*math.lgamma(k + 1) + 4*k*math.log(396) total += math.exp(ln_constant + numerator - denominator) return 1 / total # 仅需5次迭代即可达到15位精度 print(ramanujan_pi(5)) # 输出3.141592653589793

3. 算法选择指南

3.1 精度需求场景

不同算法的收敛速度差异巨大：

教学建议：初学者建议从蒙特卡洛法入手，直观理解概率与面积关系

3.2 计算资源考量

• 受限环境（如嵌入式设备）：选择割圆法
• 多核服务器：并行蒙特卡洛法
• 超高精度需求：拉马努金公式

实际案例：在树莓派上测试，拉马努金公式计算100位π仅需2秒，而蒙特卡洛法达到相同精度需要1小时

4. 调试技巧与常见陷阱

4.1 浮点数精度问题

所有算法都会遇到浮点精度限制。解决方法：

• 使用Python的decimal模块
• 对关键中间结果取对数运算
• 避免大数相减（如1 - 0.999999）

示例改进：

from decimal import Decimal, getcontext def precise_archimedes(iterations): getcontext().prec = 50 # 设置50位精度 sides = 6 length = Decimal(1) for _ in range(iterations): term = 1 - (1 - (length/2)**2).sqrt() length = (term**2 + (length/2)**2).sqrt() sides *= 2 return length * sides / 2

4.2 随机性控制

蒙特卡洛法中：

• 设置随机种子保证结果可复现
• 注意random模块的线程安全问题
• 大样本时考虑使用numpy.random提高性能

import numpy as np def np_monte_carlo(points): np.random.seed(42) # 固定种子 x = np.random.random(points) y = np.random.random(points) return 4 * np.sum(x**2 + y**2 <= 1) / points

4.3 收敛性判断

无穷级数法需要特别注意停止条件。改进方案：

def leibniz_pi(tolerance): pi_estimate = 0 denominator = 1 sign = 1 while True: term = sign / denominator pi_estimate += term if abs(term) < tolerance: break sign *= -1 denominator += 2 return 4 * pi_estimate

5. 扩展应用与可视化

5.1 割圆法动态演示

使用matplotlib展示逼近过程：

import matplotlib.pyplot as plt import numpy as np def plot_archimedes(iterations): fig, ax = plt.subplots(figsize=(8,8)) ax.set_aspect('equal') sides = 6 * 2**iterations angles = np.linspace(0, 2*np.pi, sides+1) x = np.cos(angles) y = np.sin(angles) ax.plot(x, y, label=f'{sides}边形') ax.set_title(f'割圆法逼近π ≈ {np.pi:.10f}') plt.legend() plt.show() plot_archimedes(5) # 展示5次分割后的效果

5.2 蒙特卡洛法散点图

直观展示随机点分布：

def plot_monte_carlo(points): inside_x, inside_y = [], [] outside_x, outside_y = [], [] for _ in range(points): x, y = random.random(), random.random() if x**2 + y**2 <= 1: inside_x.append(x) inside_y.append(y) else: outside_x.append(x) outside_y.append(y) plt.scatter(inside_x, inside_y, color='blue', s=1) plt.scatter(outside_x, outside_y, color='red', s=1) plt.title(f'π ≈ {4*len(inside_x)/points:.5f} (n={points})') plt.show() plot_monte_carlo(5000)

6. 进阶话题：算法变种与创新

6.1 布丰针问题

另一种基于概率的π计算方法：在画有平行线的平面上随机投针，通过相交概率计算π：

def buffon_needle(trials, line_distance=2, needle_length=1): crosses = 0 for _ in range(trials): y = random.uniform(0, line_distance/2) angle = random.uniform(0, math.pi/2) if y <= (needle_length/2)*math.sin(angle): crosses += 1 return (2*needle_length*trials) / (line_distance*crosses)

6.2 Chudnovsky算法

目前最快的π计算算法之一，每项可产生约14位正确小数：

from decimal import Decimal, getcontext def chudnovsky(precision): getcontext().prec = precision + 2 C = 426880 * Decimal(10005).sqrt() M = 1 L = 13591409 X = 1 K = 6 S = L for i in range(1, precision//14 + 2): M = (K**3 - 16*K) * M // (i+1)**3 L += 545140134 X *= -262537412640768000 S += Decimal(M * L) / X K += 12 return C / S