优化类梅钦公式的π近似计算：提升性能、最小化重复计算及精度评估方法

我想用类梅钦公式来计算π，为了让近似收敛得更快，我采用了牛顿加速级数：

def Newton_arctan(x: int | float, lim: int) -> float:
    if not (lim and isinstance(lim, int)):
        raise ValueError(f"Argument lim must be a positive integer, received {lim}")

    square = x**2
    y = y_0 = 1 + square
    even_p = even = 2
    odd_p = odd = 3
    s = x / y
    for _ in range(lim - 1):
        s += even_p / odd_p * (x := x * square) / (y := y * y_0)
        even += 2
        odd += 2
        even_p *= even
        odd_p *= odd

    return s


def Machin_Pi_worker(terms: list, lim: int) -> float:
    return 4 * sum(coef * Newton_arctan(1 / denom, lim) for coef, denom in terms)


def Machin_Pi1(lim: int) -> float:
    return Machin_Pi_worker(((4, 5), (-1, 239)), lim)

我做了如下测试：

In [178]: old = Machin_Pi1(i := 1)
     ...: while True:
     ...:     if (new := Machin_Pi1(i := i + 1)) == old:
     ...:         break
     ...:
     ...:     old = new

In [179]: i -= 1; print(i, Machin_Pi1(i))
11 3.141592653589793

这次同样用了11次迭代就达到了最大精度，所有数字都是正确的——不过这次只有15位小数，有意思的是这个值和math.pi完全相同。

我还尝试了其他一些类梅钦公式：

def Machin_Pi2(lim: int) -> float:
    return Machin_Pi_worker(((6, 8), (2, 57), (1, 239)), lim)


def Machin_Pi3(lim: int) -> float:
    return Machin_Pi_worker(((12, 18), (8, 57), (-5, 239)), lim)


def Machin_Pi4(lim: int) -> float:
    return Machin_Pi_worker(((12, 49), (32, 57), (-5, 239), (12, 110443)), lim)


def Machin_Pi5(lim: int) -> float:
    return Machin_Pi_worker(((44, 57), (7, 239), (-12, 682), (24, 12943)), lim)


test_result = {}
for i in range(1, 6):
    func = globals()[fname := f"Machin_Pi{i}"]
    old = func(j := 1)
    while True:
        if (new := func(j := j + 1)) == old:
            break

        old = new

    test_result[fname] = (new, j - 1)

测试结果如下：

{'Machin_Pi1': (3.141592653589793, 11),
 'Machin_Pi2': (3.1415926535897936, 9),
 'Machin_Pi3': (3.1415926535897927, 7),
 'Machin_Pi4': (3.1415926535897927, 5),
 'Machin_Pi5': (3.1415926535897922, 5)}

后面的公式收敛得更快，但它们在达到双精度浮点数的最大可能精度之前就出现了浮点数下溢的问题。

于是我想，为了最小化浮点数下溢的影响，应该把分子和分母分开作为整数计算，这样在最终除法之前就不会丢失精度。

我已经很多年没动笔算数学题了，数学功底也生疏了不少，但还是推导了一下：

之后我重新实现了整个计算逻辑：

from typing import List, Tuple

Fraction = Tuple[int, int]


def Newton_arctan_xr(i: int | float, lim: int) -> float:
    if not (lim and isinstance(lim, int)):
        raise ValueError(f"Argument lim must be a positive integer, received {lim}")

    cur_hi = dividend = i_sqr = i * i
    i_sqr_p = i_sqr + 1
    divisor = i * i_sqr_p
    even = 2
    odd = 3
    for _ in range(lim - 1):
        cur_hi *= even * i_sqr
        divisor *= (prod := odd * i_sqr * i_sqr_p)
        dividend = dividend * prod + cur_hi
        even += 2
        odd += 2

    return dividend, divisor


def add_fractions(frac1: Fraction, frac2: Fraction) -> Fraction:
    a, b = frac1
    c, d = frac2
    return (a * d + b * c, b * d)


def sum_fractions(fractions: List[Fraction]) -> Fraction:
    result = fractions[0]
    for frac in fractions[1:]:
        result = add_fractions(result, frac)

    return result


def gcd(x: int, y: int) -> int:
    while y != 0:
        (x, y) = (y, x % y)

    return x


def Machin_Pi_worker1(terms: List[Tuple[int, int]], lim: int) -> Fraction:
    fractions = []
    for coef, inv in terms:
        dividend, divisor = Newton_arctan_xr(inv, lim)
        fractions.append((coef * dividend, divisor))

    dividend, divisor = sum_fractions(fractions)
    dividend *= 4
    extra = gcd(dividend, divisor)
    return dividend // extra, divisor // extra


def Machin_Pi_1(lim: int) -> Fraction:
    return Machin_Pi_worker1(((4, 5), (-1, 239)), lim)


def Machin_Pi_2(lim: int) -> Fraction:
    return Machin_Pi_worker1(((6, 8), (2, 57), (1, 239)), lim)


def Machin_Pi_3(lim: int) -> Fraction:
    return Machin_Pi_worker1(((12, 18), (8, 57), (-5, 239)), lim)


def Machin_Pi_4(lim: int) -> Fraction:
    return Machin_Pi_worker1(((12, 49), (32, 57), (-5, 239), (12, 110443)), lim)


def Machin_Pi_5(lim: int) -> Fraction:
    return Machin_Pi_worker1(((44, 57), (7, 239), (-12, 682), (24, 12943)), lim)

测试这个实现的结果如下：

In [230]: Machin_Pi_5(5)
Out[230]:
(1279457632672435538478197124236187110232840682131383545616,
 407264013432945209516129385309101616710788249969482421875)

In [231]: 1279457632672435538478197124236187110232840682131383545616/407264013432945209516129385309101616710788249969482421875
Out[231]: 3.141592653589793

这个实现是有效的，但我不知道有没有把重复计算降到最少，也不确定用什么库能加速执行——不过我并不是在求软件推荐，所以你可以像我一样用纯Python实现，但要求你的代码运行速度比我写的快。

另外我很想得到更多位数的π值，我用了gmpy2.mpfr，但不知道该如何准确估计某个分数对应的正确位数，这样我就能相应地传递精度参数了。

备注：内容来源于stack exchange，提问作者Ξένη Γήινος