如何使用类似机械的公式改善π的近似？ 这是π天，所以我认为将π的近似值作为编程挑战是美好的一天。 今年，我想使用类似Machin的公式来计算它，我将牛顿的加速系列用于MA ...

Question

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.。

我尝试了一些其他公式，来自

here的：

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) -> float:
    return Machin_Pi_worker1(((4, 5), (-1, 239)), lim)


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


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


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


def Machin_Pi_5(lim: int) -> float:
    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

，我想知道如何准确估计特定分数具有的正确数字的数量，因此我们可以相应地通过Precision参数。 解决此问题的标准方法是使用

核

内置软件包。

import decimal as dc

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

    square = x*x
    y = y_0 = 1 + square
    even_p = even = 2
    odd_p = odd = 3
    s = x / y

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

    return s


def Machin_Pi_worker(terms: list[int], lim: int) -> dc.Decimal:
    return 4 * sum(coef * Newton_arctan(1 / dc.Decimal(denom), lim) for coef, denom in terms)

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

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

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

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

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


if __name__ == "__main__":
    # setting context
    import sys
    if (len(sys.argv) >= 2) and sys.argv[1].isdecimal():
        dc.getcontext().prec = int(sys.argv[1])  # precision
    else:
        dc.getcontext().prec = 15  # default precision
    dc.getcontext().traps[dc.FloatOperation] = True  # prevents mixing floats and decimals in operations


    # running test
    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)
    print(test_result)

运行此代码时：

$ python 79507720.py
{
'Machin_Pi1': (Decimal('3.14159265358976'), 10),
'Machin_Pi2': (Decimal('3.14159265358978'), 8),
'Machin_Pi3': (Decimal('3.14159265358978'), 6),
'Machin_Pi4': (Decimal('3.14159265358979'), 5),
'Machin_Pi5': (Decimal('3.14159265358980'), 5)
}

$ python 79507720.py 32
{
'Machin_Pi1': (Decimal('3.1415926535897932384626433832793'), 22), 
'Machin_Pi2': (Decimal('3.1415926535897932384626433832792'), 17), 
'Machin_Pi3': (Decimal('3.1415926535897932384626433832795'), 13), 
'Machin_Pi4': (Decimal('3.1415926535897932384626433832794'), 10), 
'Machin_Pi5': (Decimal('3.1415926535897932384626433832795'), 9)
}

$ python 79507720.py 80
{
'Machin_Pi1': (Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862092'), 56), 
'Machin_Pi2': (Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862086'), 44), 
'Machin_Pi3': (Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862090'), 32), 
'Machin_Pi4': (Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862092'), 24), 
'Machin_Pi5': (Decimal('3.1415926535897932384626433832795028841971693993751058209749445923078164062862093'), 23)
}

Answer 1