[Python] How to improve the approximation of π using Machin-like formulas?

It is π-day so I think it is a good day to compute approximations of π as a programming challenge.

Of course we don't need astronomically many digits of π, but as you may well know this impracticality doesn't stop people from computing trillions of digits of π and keep setting world records of π approximation lengths. And as a matter of fact, I know computing π is a popular programming challenge, and I think trying to compute π may help one become a better programmer.

This year I want to compute it using Machin-like formulae, I chose these in homage of William Shanks, he computed 707 digits of decimal expansion of π, in the age well before electronic computers, it was tremendously difficult, and he calculated 527 correct digits. Why did he do it? Well, why do I do it?

So, naive implementation:

def arctan_series(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 * x
sign = 1
return x + sum(
(x := x * square) / i * (sign := sign * -1) for i in range(3, 2 * lim + 1, 2)
)


def Machin_Pi(lim: int) -> float:
return 16 * arctan_series(0.2, lim) - 4 * arctan_series(1 / 239, lim)


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

In [176]: i
Out[176]: 12

In [177]: Machin_Pi(11)
Out[177]: 3.1415926535897936


As you can see, it took 11 iterations to reach the maximum precision before floating point inaccuracy kicks in, and the last digit should be 2. And it is well known the Taylor series of arctangent converges extremely slowly.

So I used Newton's accelerated series instead:

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:



I am not a native English speaker so the following may seem strange to you.

In order to keep the fractions as fractions we need to unify the denominator, to unify the denominator we have to find the least common multiple of all the denominators (though for simplicity I didn't do that as you can see), then multiply the numerator by least common multiple divided by the denominator.

In this way we can simply add up the new numerators as they have the same denominator. And the reciprocals of powers of 5 should really go into the denominator so that we keep the numerator and denominator as separate integers.

You can see the result on the left.

Now we can compute 5^26 from 5^24 and then 5^28 from 5^26 et cetera, we can compute all these while computing 5^35 so as to minimize the redundant computations.

And in every addition term we are computing the product of every other odd number except one odd number, this is terribly inefficient, there are a lot of repeated calculations, we should try to eliminate these repeated calculations. And during the computation of the denominator we already did a good chunk of computation need for the numerators.

I wanted to find a chain of multiplications that does all the required multiplications and minimizes repeated multiplications but I can't:

odds = list(range(3, 14, 2))
odd_pairs = [(e, f) for i, e in enumerate(odds[:-1]) for f in odds[i + 1 :]]


def recur_find_runs_worker(tree: dict, excluded: set = set()) -> None:
for pair in odd_pairs:
if not (duo := set(pair)) & excluded:
tree[pair] = {}
recur_find_runs_worker(tree[pair], excluded | duo)


def flatten(tree: dict) -> list:
result = []
for k, v in tree.items():
if not v:
result.append([k])
else:
result.extend([[k] + i for i in flatten(v)])

return result


def recur_find_runs(excluded: set = set()) -> list:
tree = {}
recur_find_runs_worker(tree, excluded)
return flatten(tree)


every_other = {e: odds[:i] + odds[i + 1 :] for i, e in enumerate(odds)}


I got stuck and I am unsure what to do next.

What this question is really about, is rewriting Newton_arctan so that throughout the execution no float is used, all calculations are in int domain, the numerator and denominator are computed separately, so that there isn't precision loss. Because it doesn't do division.

In short Newton_arctan should return a pair of integers, and it should do so with least redundant computations. This is what I am asking.

You must use Machin-like formulae with Taylor-like series to do the approximation, like I do in this question.

And finally, you can use library code in your answer, but you can't use arctangent and any other trigonometry function for obvious reasons, that would defeat the purpose of this exercise. And of course you can't do things like math.pi. (And if you use library code, it would be better if you can get more than 15 (16?) correct digits, though this isn't a requirement.)

Continue reading...