[Python] Why is the continued fraction expansion of arctangent combined with half-angle...

Sorry for the long title. I don't know if this is more of a math problem or programming problem, but I think my math is extremely rusty and I am better at programming.

So I have this continued fraction expansion of arctangent:



I got it from Wikipedia

I tried to find a simple algorithm to calculate it:



And I did it, I have written an infinite precision implementation of the continued fraction expansion without using any libraries, using only basic integer arithmetic:

import json
import math
import random
from decimal import Decimal, getcontext
from typing import Callable, List, Tuple


Fraction = Tuple[int, int]


def arctan_cf(y: int, x: int, lim: int) -> Fraction:
y_sq = y**2
a1, a2 = y, 3 * x * y
b1, b2 = x, 3 * x**2 + y_sq
odd = 5
for i in range(2, 2 + lim):
t1, t2 = odd * x, i**2 * y_sq
a1, a2 = a2, t1 * a2 + t2 * a1
b1, b2 = b2, t1 * b2 + t2 * b1
odd += 2

return a2, b2


And it converges faster than Newton's arctangent series which I previously used.

Now I think if I combine it with the half-angle formula of arctangent it should converge faster.

def half_arctan_cf(y: int, x: int, lim: int) -> Fraction:
c = (x**2 + y**2) ** 0.5
a, b = c.as_integer_ratio()
a, b = arctan_cf(a - b * x, b * y, lim)
return 2 * a, b


And indeed, it does converge even faster:

def test_accuracy(lim: int) -> dict:
result = {}
for _ in range(lim):
x, y = random.sample(range(1024), 2)
while not x or not y:
x, y = random.sample(range(1024), 2)

atan2 = math.atan2(y, x)
entry = {"atan": atan2}
for fname, func in zip(
("arctan_cf", "half_arctan_cf"), (arctan_cf, half_arctan_cf)
):
i = 1
while True:
a, b = func(y, x, i)
if math.isclose(deci := a / b, atan2):
break

i += 1

entry[fname] = (i, deci)

result[f"{y} / {x}"] = entry

return result


print(json.dumps(test_accuracy(8), indent=4))

for v in test_accuracy(128).values():
assert v["half_arctan_cf"][0] <= v["arctan_cf"][0]


{
"206 / 136": {
"atan": 0.9872880750087898,
"arctan_cf": [
16,
0.9872880746658675
],
"half_arctan_cf": [
6,
0.9872880746018052
]
},
"537 / 308": {
"atan": 1.0500473287277563,
"arctan_cf": [
18,
1.0500473281360896
],
"half_arctan_cf": [
7,
1.0500473288158192
]
},
"331 / 356": {
"atan": 0.7490241118247137,
"arctan_cf": [
10,
0.7490241115996227
],
"half_arctan_cf": [
5,
0.749024111913438
]
},
"744 / 613": {
"atan": 0.8816364228048325,
"arctan_cf": [
13,
0.8816364230439662
],
"half_arctan_cf": [
6,
0.8816364227495634
]
},
"960 / 419": {
"atan": 1.1592605364805093,
"arctan_cf": [
24,
1.1592605359263286
],
"half_arctan_cf": [
7,
1.1592605371181872
]
},
"597 / 884": {
"atan": 0.5939827714677137,
"arctan_cf": [
7,
0.5939827719895824
],
"half_arctan_cf": [
4,
0.59398277135389
]
},
"212 / 498": {
"atan": 0.40246578425167584,
"arctan_cf": [
5,
0.4024657843859885
],
"half_arctan_cf": [
3,
0.40246578431841773
]
},
"837 / 212": {
"atan": 1.322727785860997,
"arctan_cf": [
41,
1.322727786922624
],
"half_arctan_cf": [
8,
1.3227277847674388
]
}
}


That assert block runs quite a bit long for large number of samples, but it never raises exceptions.

So I think I can use the continued fraction expansion of arctangent with Machin-like series to calculate π. (I used the last series in the linked section because it converges the fastest)

def sum_fractions(fractions: List[Fraction]) -> Fraction:
while (length := len(fractions)) > 1:
stack = []
for i in range(0, length - (odd := length & 1), 2):
num1, den1 = fractions
num2, den2 = fractions[i + 1]
stack.append((num1 * den2 + num2 * den1, den1 * den2))

if odd:
stack.append(fractions[-1])

fractions = stack

return fractions[0]


MACHIN_SERIES = ((44, 57), (7, 239), (-12, 682), (24, 12943))


def approximate_loop(lim: int, func: Callable) -> List[Fraction]:
fractions = []
for coef, denom in MACHIN_SERIES:
dividend, divisor = func(1, denom, lim)
fractions.append((coef * dividend, divisor))

return fractions


def approximate_1(lim: int) -> List[Fraction]:
return approximate_loop(lim, arctan_cf)


def approximate_2(lim: int) -> List[Fraction]:
return approximate_loop(lim, half_arctan_cf)


approx_funcs = (approximate_1, approximate_2)


def calculate_pi(lim: int, approx: bool = 0) -> Fraction:
dividend, divisor = sum_fractions(approx_funcs[approx](lim))
dividend *= 4
return dividend // (common := math.gcd(dividend, divisor)), divisor // common

getcontext().rounding = 'ROUND_DOWN'
def to_decimal(dividend: int, divisor: int, places: int) -> str:
getcontext().prec = places + len(str(dividend // divisor))
return str(Decimal(dividend) / Decimal(divisor))


def get_accuracy(lim: int, approx: bool = 0) -> Tuple[int, str]:
length = 12
fraction = calculate_pi(lim, approx)
while True:
decimal = to_decimal(*fraction, length)
for i, e in enumerate(decimal):
if Pillion != e:
return (max(0, i - 2), decimal[:i])

length += 10


with open("D:/Pillion.txt", "r") as f:
Pillion = f.read()


Pillion.txt contains the first 1000001 digits of π, Pi + Million = Pillion.

And it works, but only partially. The basic continued fraction expansion works very well with Machin-like formula, but combined with half-angle formula, I can only get 9 correct decimal places no matter what, and in fact, I get 9 correct digits on the very first iteration, and then this whole thing doesn't improve ever:

In [2]: get_accuracy(16)
Out[2]:
(73,
'3.1415926535897932384626433832795028841971693993751058209749445923078164062')

In [3]: get_accuracy(32)
Out[3]:
(138,
'3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231')

In [4]: get_accuracy(16, 1)
Out[4]: (9, '3.141592653')

In [5]: get_accuracy(32, 1)
Out[5]: (9, '3.141592653')

In [6]: get_accuracy(1, 1)
Out[6]: (9, '3.141592653')


But the digits do in fact change:

In [7]: to_decimal(*calculate_pi(1, 1), 32)
Out[7]: '3.14159265360948500093515231500093'

In [8]: to_decimal(*calculate_pi(2, 1), 32)
Out[8]: '3.14159265360945286794831052938917'

In [9]: to_decimal(*calculate_pi(3, 1), 32)
Out[9]: '3.14159265360945286857612896472974'

In [10]: to_decimal(*calculate_pi(4, 1), 32)
Out[10]: '3.14159265360945286857611676794770'

In [11]: to_decimal(*calculate_pi(5, 1), 32)
Out[11]: '3.14159265360945286857611676818392'


Why is the continued fraction with half-angle formula not working with Machin-like formula? And is it possible to make it work, and if it can work, then how? I want either a proof that it is impossible, or a working example that proves it is possible.


Just a sanity check, using π/4 = arctan(1) I was able to make half_arctan_cf spit out digits of π but it converges much slower:

def approximate_3(lim: int) -> List[Fraction]:
return [half_arctan_cf(1, 1, lim)]


approx_funcs = (approximate_1, approximate_2, approximate_3)


In [28]: get_accuracy(16, 2)
Out[28]: (15, '3.141592653589793')

In [29]: get_accuracy(16, 0)
Out[29]:
(73,
'3.1415926535897932384626433832795028841971693993751058209749445923078164062')


And the same problem recurs, it reaches maximum precision of 15 digits at the 10th iteration:

In [37]: get_accuracy(9, 2)
Out[37]: (14, '3.14159265358979')

In [38]: get_accuracy(10, 2)
Out[38]: (15, '3.141592653589793')

In [39]: get_accuracy(11, 2)
Out[39]: (15, '3.141592653589793')

In [40]: get_accuracy(32, 2)
Out[40]: (15, '3.141592653589793')

Continue reading...