python - How to improve the approximation of π using Machin-like formulas? - Stack Overflow

admin管理员组

文章数量:1390729

I want to compute π using Machin-like formulae, I use Newton's accelerated series to make the approximation converge faster:

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

This time it also took 11 iterations to reach maximum precision, and all digits are correct, though this time it has only 15 decimal places, interestingly this value is the same as math.pi.

I tried a bunch of other formulae from 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)}

The later series converged faster, but they reached float underflow before they reached the maximum possible precision in double precision floating point format.

Now I think to minimize the impact of float underflow I want to make it so the numerator and denominator are computed separately as integers, so that we don't lose precision before final division.

It has been many years since I last used pen and paper and my math is extremely rusty, but I did the following computation:

And I re-implemented the whole thing:

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

It works, but I don't know if I have reduced repeated calculations to the minimum, and I don't know what library I can use to speed up the execution, though I am not asking for software recommendations, so you can just implement in pure Python like me, but the answer is required to run faster than mine code.

Though I am really interested in getting more digits, I use gmpy2.mpfr, and I want to know how to accurately estimate the number of correct digits a particular fraction has, so we can pass precision argument accordingly.

I want to compute π using Machin-like formulae, I use Newton's accelerated series to make the approximation converge faster:

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

This time it also took 11 iterations to reach maximum precision, and all digits are correct, though this time it has only 15 decimal places, interestingly this value is the same as math.pi.

I tried a bunch of other formulae from 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)}

The later series converged faster, but they reached float underflow before they reached the maximum possible precision in double precision floating point format.

Now I think to minimize the impact of float underflow I want to make it so the numerator and denominator are computed separately as integers, so that we don't lose precision before final division.

It has been many years since I last used pen and paper and my math is extremely rusty, but I did the following computation:

And I re-implemented the whole thing:

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

It works, but I don't know if I have reduced repeated calculations to the minimum, and I don't know what library I can use to speed up the execution, though I am not asking for software recommendations, so you can just implement in pure Python like me, but the answer is required to run faster than mine code.

Though I am really interested in getting more digits, I use gmpy2.mpfr, and I want to know how to accurately estimate the number of correct digits a particular fraction has, so we can pass precision argument accordingly.

• python
• math
• pi
• taylor-series

Share Improve this question Follow edited Mar 17 at 0:19 desertnaut 60.5k3232 gold badges155155 silver badges182182 bronze badges asked Mar 13 at 21:57 Ξένη ΓήινοςΞένη Γήινος 3,64011 gold badge1818 silver badges6060 bronze badges 5

• A fairly comprehensive collection of Machin-type formulas maintained by Jörg Arndt can be found here. – njuffa Commented Mar 13 at 22:05
• 3 I'm having a hard time understanding what problem you're asking for help with. Please focus you post to be a much boring "I'm trying to do X, this is what I already did, this is the problem I'm seeing that have, and this is the debugging I already did to try to fix it". Most of the text in your post is just padding that reads well in a blog post, but does not belong in a question post. – Mike 'Pomax' Kamermans Commented Mar 14 at 3:25
• Not a definitive answer, but instead of bothering with fractions I'd use the Decimal built-in package (that does the work for you), which allows for arbitrarily large precision and behaves as expected from numbers when computing operations ; at the cost of a slight computational overhead – globglogabgalab Commented Mar 14 at 7:55
• 1 Why the downvote now? I have absolutely shown more than enough amount of research, I have deleted unnecessary text, and the downvote is cast right after I got an answer so I can't easily delete this question... – Ξένη Γήινος Commented Mar 14 at 8:47
• The vote may not be a direct response to your edit but simply because the question gained new visibility due to the edit. – mkrieger1 Commented Mar 14 at 8:50

Add a comment  | 

1 Answer 1

Sorted by: Reset to default 1

The standard way to address these precision issues is to use the decimal built-in package, that allows for arbitrary fixed-point arithmetic precision, at the cost of a (small) computational overhead.

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)

When running this code:

$ 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 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)
}

$ python 79507720.py 999
{
'Machin_Pi1': (Decimal('3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642027'), 705), 
'Machin_Pi2': (Decimal('3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642018'), 550), 
'Machin_Pi3': (Decimal('3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642020'), 397), 
'Machin_Pi4': (Decimal('3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642015'), 295), 
'Machin_Pi5': (Decimal('3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807446237996274956735188575272489122793818301194912983367336244065664308602139494639522473719070217986094370277053921717629317675238467481846766940513200056812714526356082778577134275778960917363717872146844090122495343014654958537105079227968925892354201995611212902196086403441815981362977477130996051870721134999999837297804995105973173281609631859502445945534690830264252230825334468503526193118817101000313783875288658753320838142061717766914730359825349042875546873115956286388235378759375195778185778053217122680661300192787661119590921642018'), 284)
}

本文标签： pythonHow to improve the approximation of π using Machinlike formulasStack Overflow