Wie verbessert man die Näherung von π unter Verwendung von Maschin-ähnlichen Formeln? ⇐ Python

[Anonymous]

Ich möchte π unter Verwendung von machinähnlichen Formeln berechnen. Ich verwende Newtons beschleunigte Serie, um die Näherung schneller zu konvergieren:

Code: Select all

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)
< /code>
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
< /code>
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

.
Ich habe von hier eine Reihe anderer Formeln ausprobiert:

Code: Select all

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)
< /code>
{'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)}
< /code>
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:
[img]https://i.sstatic.net/fSDC846t.jpg[/img]

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)
< /code>
In [230]: Machin_Pi_5(5)
Out[230]:
(1279457632672435538478197124236187110232840682131383545616,
407264013432945209516129385309101616710788249969482421875)

In [231]: 1279457632672435538478197124236187110232840682131383545616/407264013432945209516129385309101616710788249969482421875
Out[231]: 3.141592653589793
< /code>
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

, und ich möchte wissen, wie man die Anzahl der korrekten Ziffern, die ein bestimmter Bruch hat, genau schätzen, damit wir das Präzisionsargument entsprechend übergeben können.