Por que a expansão de fração contínua do arco tangente combinada com a fórmula do meio-ângulo não funciona com séries do tipo Machin?

Desculpe pelo título longo. Não sei se isso é mais um problema de matemática ou de programação, mas acho que minha matemática está extremamente enferrujada e sou melhor em programação.

Então eu tenho essa expansão de fração contínua do arco tangente:

Eu peguei da Wikipedia

Tentei encontrar um algoritmo simples para calculá-lo:

E eu consegui, escrevi uma implementação de precisão infinita da expansão de frações contínuas sem usar nenhuma biblioteca, usando apenas aritmética básica de inteiros:

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

E converge mais rápido que a série arco-tangente de Newton que usei anteriormente .

Agora acho que se eu combinar isso com a fórmula do meio-ângulo do arco tangente, a convergência deve ser mais rápida.

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

E, de fato, converge ainda mais rápido:

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
        ]
    }
}

Esse bloco assert é bastante longo para um grande número de amostras, mas nunca gera exceções.

Então, acho que posso usar a expansão de fração contínua do arco tangente com séries do tipo Machin para calcular π. (Usei a última série na seção vinculada porque ela converge mais rápido)

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 contém os primeiros 1000001 dígitos de π, Pi + Milhão = Pillion.

E funciona, mas apenas parcialmente. A expansão básica de fração contínua funciona muito bem com fórmulas do tipo Machin, mas combinadas com fórmulas de meio ângulo, eu só consigo obter 9 casas decimais corretas não importa o que aconteça, e de fato, eu obtenho 9 dígitos corretos na primeira iteração, e então essa coisa toda não melhora nunca:

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')

Mas os dígitos de fato mudam:

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'

Por que a fração contínua com fórmula de meio-ângulo não está funcionando com fórmula tipo Machin? E é possível fazê-la funcionar, e se puder funcionar, então como? Eu quero uma prova de que é impossível, ou um exemplo funcional que prove que é possível.

Apenas uma verificação de sanidade, usando π/4 = arctan(1) consegui gerar half_arctan_cfdígitos de π, mas ele converge muito mais lentamente:

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')

E o mesmo problema se repete, ele atinge a precisão máxima de 15 dígitos na 10ª iteração:

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')

Acabei de reescrever minha implementação de fração contínua arco-tangente e fiz com que ela evitasse cálculos redundantes.

No meu código, em cada iteração, t1 aumenta em 2 * y_sq, então não há necessidade de multiplicar repetidamente y_sq pelo número ímpar; em vez disso, basta usar uma variável cumulativa e um passo de 2 * y_sq.

E a diferença entre cada par de números quadrados consecutivos são apenas os números ímpares, então posso usar uma variável cumulativa de uma variável cumulativa.

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)

Embora isso não melhore muito a velocidade, é definitivamente uma melhoria.