Question: 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[i] 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[i] != 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')   I just rewrote my arctangent continued fraction implementation and made it avoid doing redundant computations. In my code in each iteration t1 increases by 2 * y_sq, so there is no need to repeatedly multiply y_sq by the odd number, instead just use a cumulative variable and a step of 2 * y_sq. And the difference between each pair of consecutive square numbers is just the odd numbers, so I can use a cumulative variable of a cumulative variable. def arctan_cf_0(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 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 t1_step, t3_step = 2 * x, 2 * y_sq t1, t2 = 5 * x, 4 * y_sq t3 = t2 + y_sq for _ in range(lim): a1, a2 = a2, t1 * a2 + t2 * a1 b1, b2 = b2, t1 * b2 + t2 * b1 t1 += t1_step t2 += t3 t3 += t3_step return a2, b2  In [301]: arctan_cf_0(4, 3, 100) == arctan_cf(4, 3, 100) Out[301]: True In [302]: %timeit arctan_cf_0(4, 3, 100) 58.6 μs ± 503 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each) In [303]: %timeit arctan_cf(4, 3, 100) 54.3 μs ± 816 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)  While this doesn't improve the speed by much, this is definitively an improvement.

Answer: The loss of precision is here: c = (x**2 + y**2) ** 0.5 a, b = c.as_integer_ratio()  This code works with float type, which is the result of power calculation ** 0.5. As soon as you use float, you lose arbitrary precision. You might want to use the Decimal.sqrt method instead.