Python计算Pi：马青公式的高效实现技巧

一、马青公式的数学本质与加速原理

要高效率地计算圆周率，最核心的抉择在于选取收敛速度快且各项分母易于整数分解的展开式。初等微积分中最著名的反正切级数是格里高利-莱布尼茨级数（Gregory-Leibniz Series）：

arctan(x) = x - x³/3 + x⁵/5 - x⁷/7 + ... = ∑ (-1)ᵏ · x^(2k+1) / (2k+1)

当取 x = 1 时，可直接得出 π/4 = 1 - 1/3 + 1/5 - 1/7 + ...。然而，该级数的公比绝对值恒等于 1，仅靠分母奇数线性递增来缩小余项。数学推导表明，若想通过莱布尼茨公式计算出 1,000 位有效小数，需要累加约 10¹⁰⁰⁰ 项，在现代计算机硬件下完全不可行。

英国天文学家约翰·马青（John Machin）在 1706 年提出了著名的正切和角恒等式分解：

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

两边同乘 4，即可得到计算 π 的显式表达：

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

这一公式为何能够带来颠覆性的速度提升？原因在于它将自变量 x = 1 转化为了非常接近 0 的两个有理数 1/5 和 1/239：

• 分支一 arctan(1/5)：每推进一项，自变量乘方衰减 1/5² = 1/25 = 0.04。因为 log10(25) ≈ 1.398，意味着每计算一项，就能稳定获得约 1.4 位十进制有效数字。计算 1,000 位仅需约 715 项。
• 分支二 arctan(1/239)：每推进一项，自变量乘方衰减 1/239² = 1/57121 ≈ 1.75 × 10⁻⁵。因为 log10(57121) ≈ 4.757，每计算一项就能推进将近 4.8 位十进制精度，计算 1,000 位仅需约 210 项。

二、规避浮点陷阱：定点大整数缩放体系

在标准 Python 中，普通的浮点数 float 是基于 IEEE 754 双精度浮点格式实现的，其尾数仅有 53 位二进制位，对应十进制有效精度仅为 15 到 17 位。一旦超出此范围，尾数将发生不可逆的硬件截断。虽然标准库中的 decimal 模块支持任意精度，但对象开销和上下文管理较为厚重。

最轻量、执行速度最快的工程实践是定点大整数缩放（Fixed-Point Arithmetic）：

1. 设定目标位数与放大因子：若要计算 N 位小数，我们首先构造一个放大因子 Scale = 10^(N + G)，其中 G 为安全保护位（Guard Digits，通常设定为 8~10 位）。保护位的作用是吸收多次除法舍入累加产生的小数点末位误差。
2. 放大级数计算对象：将 arctan(1/x) 的求和整体放大 Scale 倍。原级数的首项分子直接变为 Scale // x，后续项的分子则每次整除 x²。
3. 剥离保护位并格式化：求得放大后的整数 π_scaled 后，通过纯整数除法 π_scaled // 10^G 剔除多余的保护位，再格式化输出为包含小数点的标准字符串。

三、核心算法实现与项递推优化

在编写级数求和循环时，初学者最常见的性能陷阱是每一轮循环都通过幂运算符 x ** (2 * k - 1) 从头计算分母幂次。这会导致大整数幂的计算量随着迭代深度线性膨胀，使时间复杂度恶化至 O(N³)。

高效的做法是项递推（Term Recurrence）：维护一个分子变量 num。第 1 项的分子是 Scale // x；此后每进入下一轮迭代，只需做一次单常数整除 num //= (x * x)。由于 x² 分别是固定的常数 25 和 57121，在 CPython 底层，大整数除以单字长整数（Single-word integer division）具有极高的机器执行效率。

完整 Python 实现代码

以下为优化后的完整生产级函数，包含自适应零项终止与精度截断处理：

# -*- coding: utf-8 -*-
"""
Python 高精度计算 Pi：马青公式优化实现
算法：pi/4 = 4 * arctan(1/5) - arctan(1/239)
"""

def compute_pi_machin(digits: int) -> str:
    """
    使用马青公式与定点大整数运算计算指定位数的圆周率 Pi。

    :param digits: 期望输出的小数点后有效位数 (例如 1000)
    :return: 包含小数点的高精度圆周率字符串 (如 '3.14159...')
    """
    if digits <= 0:
        raise ValueError("计算位数 digits 必须为正整数")

    # 1. 设定安全保护位，抵御递推累加中的舍入误差
    guard_digits = 10
    scale = 10 ** (digits + guard_digits)

    # 2. 定点缩放反正切子函数：计算 scale * arctan(1/x)
    def arctan_scaled(x: int) -> int:
        x2 = x * x              # 常数递推底数 (x=5 时为 25，x=239 时为 57121)
        total = 0               # 累加总和
        num = scale // x        # 初项分子：scale / x
        k = 1                   # 项序计数器 (对应分母奇数 2k - 1)
        sign = 1                # 正负符号交替器

        while num > 0:
            term = num // (2 * k - 1)   # 当前项贡献值
            if term == 0:               # 自适应终止：当前贡献已落入精度保护阈值以下
                break
            total += sign * term
            num //= x2                  # 递推传递：无需幂乘，仅消耗单次常数整除
            sign = -sign                # 翻转符号
            k += 1

        return total

    # 3. 按马青公式加权合成：pi = 16 * arctan(1/5) - 4 * arctan(1/239)
    pi_scaled = 16 * arctan_scaled(5) - 4 * arctan_scaled(239)

    # 4. 剥离末尾 10 位保护位，得到纯大整数结果
    pi_int = pi_scaled // (10 ** guard_digits)

    # 5. 格式化输出为标准浮点数字符串
    pi_str = str(pi_int)
    return f"{pi_str[0]}.{pi_str[1:digits + 1]}"

代码关键细节说明

• 整除运算符 //：严格使用 Python 的地板除法操作符，确保参与计算的操作数全程保持为纯 int 类型，不发生隐式类型转换为 float。
• term == 0 退出条件：随着 num 的不断折半衰减，最终 num // (2 * k - 1) 会等于 0。此时说明后续所有级数项对当前设定的精度窗口已没有任何数值改变，循环自动安全退出，无需冗余估算最大迭代次数。

四、运行示例与精度结果验证

我们可以通过简单的调用命令来验证实现的正确性。以下展示计算 50 位精度并与数学界公认的 Pi 基准常数进行比对的操作过程：

# 验证执行脚本
if __name__ == "__main__":
    result_50 = compute_pi_machin(50)
    print("计算结果:", result_50)

    # 业界公认的标准圆周率前 50 位小数基准
    expected_50 = "3.14159265358979323846264338327950288419716939937510"
    
    assert result_50 == expected_50, "精度校验未通过！"
    print("校验状态: 完美匹配 (50位完全一致)")

终端控制台中的实际运行输出示例为：

计算结果: 3.14159265358979323846264338327950288419716939937510
校验状态: 完美匹配 (50位完全一致)

进一步对更大规模的精度进行实测，记录各级别精度在普通单核 CPU 下的耗时表现：

五、开发注意事项与常见问题排查

1. 为什么必须保留保护位（Guard Digits）？

在纯整数定点除法中，每次 num // (2k - 1) 都会向下舍入丢弃余数。当计算 10,000 位精度时，经历数千次累加，这些舍入误差如果不断在最末尾位累加，会导致最后 2~3 位数字产生进位偏差。通过在 Scale 初始化时额外放大 10^10（即 10 位安全保护余量），所有舍入误差都被牢牢锁在保护位区间内。最后一步通过 // (10**10) 整除剥离，便能确保目标位数 100% 准确无误。

2. 超过万位后的字符串转换限制（sys.set_int_max_str_digits）

import sys
# 允许转换高达 100,000 位的大整数为字符串
sys.set_int_max_str_digits(100000)

3. 马青公式与 Chudnovsky 算法的选型对比

读者在深入探索圆周率算法时经常会听说用于创造世界纪录的 Chudnovsky 算法。在实际选型中需明确适用场景：

• 马青公式（当前方案）：原理直观、仅需 20 余行纯 Python 代码即可完成，不依赖任何第三方二进制库或复杂的超几何展开，在 100 ~ 50,000 位 的计算需求下兼具极致的开发效率与极低的学习成本。
• Chudnovsky 算法：每项可提供约 14.18 位十进制精度，适合 数百万至数千亿位 的超大规模极值挑战，但实现上必须深度结合快速傅里叶变换（FFT）大数乘法和二进制拆分（Binary Splitting），工程复杂度极其高昂。

因此，在日常工程验证、教学演示、科学基准测试以及中等规模高精度数学建模中，马青公式配合定点整数优化无疑是综合性价比最高的实现选择。