[SOLVED] Why is the continued fraction expansion of arctangent combined with half-angle formula not working with Machin-like series?

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.