使用阿基米德几何法估算圆周率时循环报错求排查

核心问题分析 #

出现inf/nan及invalid value encountered in double_scalars错误，本质是大n值下的浮点数计算溢出或精度丢失，结合阿基米德法的计算逻辑，具体原因包括：

• 直接计算大n对应的三角函数（如sin(π/n)、tan(π/n)）时，当n过大（超过双精度浮点数的有效范围，约2^53），π/n的数值会小到无法被精确表示，导致三角函数计算结果异常，甚至触发溢出。
• 若循环未设置n的上限，会导致n无限增大，最终超出浮点数能表示的最大整数范围，直接变成inf，后续基于n的计算（如2n*sin(π/n)）自然也会变成inf，误差(S2-S1)/2则因inf-inf的无意义运算变成nan。
• 若采用逐次加1的方式增大n，收敛速度极慢，会在达到精度要求前就触发数值异常；而阿基米德法原本的边数翻倍递推逻辑被忽略，进一步加剧了数值稳定性问题。

具体解决方案 #

1. 改用阿基米德递推公式替代直接三角函数计算
   放弃直接用n计算周长，改用边数翻倍的递推关系（从n=6开始，这是阿基米德最初的选择），数值稳定性远高于直接计算大n的三角函数：

   import math
   
   # 初始值：n=6时的内接、外接多边形边长（单位圆）
   a = 1.0  # 内接6边形边长为1
   b = 2 / math.sqrt(3)  # 外接6边形边长
   error = float('inf')
   n = 6
   target_precision = 1e-6
   
   while error > target_precision:
       # 递推计算2n边形的内接、外接边长
       a_new = math.sqrt(2 - 2 * math.sqrt(1 - (a/2)**2))
       b_new = 2 * a / (a + math.sqrt(4 - a**2))
       # 计算周长对应的π近似值
       pi_inner = n * a
       pi_outer = n * b
       error = (pi_outer - pi_inner) / 2
       # 更新参数
       a, b = a_new, b_new
       n *= 2
       # 防止极端情况，设置n上限
       if n > 1e10:
           break
   
   print(f"满足精度的n: {n}")
   print(f"π的近似范围: [{pi_inner}, {pi_outer}]")
   print(f"误差: {error}")
   

2. 添加循环终止的双重条件
   除了精度判断，必须设置n的最大值（如1e10），避免因收敛过慢导致无限循环和数值溢出。

3. 避免逐次加1的n增长方式
   阿基米德法通过边数翻倍（n→2n）的方式收敛，每一次迭代精度都会显著提升，远快于逐次加1，能在n达到数值危险值前就满足精度要求。

4. 排查中间计算过程
   若仍有问题，在循环中打印每一步的n、pi_inner、pi_outer、error，定位首次出现异常的步骤，确认是哪一步的计算触发了数值错误。

内容的提问来源于stack exchange，提问作者gal peled