Como implementar corretamente a fatoração de Fermat em Python?

Estou tentando implementar algoritmos eficientes de fatoração de primos em Python. Isso não é tarefa de casa ou trabalho, é completamente por curiosidade.

Aprendi que a fatoração prima é muito difícil :

Quero implementar algoritmos eficientes para isso como um desafio autoimposto. Eu me propus a implementar o método de fatoração de Fermat primeiro, pois parece simples o suficiente.

Código Python traduzido diretamente do pseudocódigo:

def Fermat_Factor(n):
    a = int(n ** 0.5 + 0.5)
    b2 = abs(a**2 - n)
    while int(b2**0.5) ** 2 != b2:
        a += 1
        b2 = a**2 - n

    return a - b2**0.5, a + b2**0.5

(Tenho que usar, abscaso contrário b2, será facilmente negativo e into cast falhará TypeErrorporque a raiz é complex)

Como você pode ver, ele retorna dois inteiros cujo produto é igual à entrada, mas ele retorna apenas duas saídas e não garante a primalidade dos fatores. Não tenho ideia de quão eficiente esse algoritmo é, mas a fatoração de semiprimos usando esse método é muito mais eficiente do que o método de divisão por tentativa usado na minha pergunta anterior: Por que a fatoração de produtos de primos próximos é muito mais lenta do que produtos de primos diferentes .

In [20]: %timeit FermatFactor(3607*3803)
2.1 μs ± 28.2 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

In [21]: FermatFactor(3607*3803)
Out[21]: [3607, 3803]

In [22]: %timeit FermatFactor(3593 * 3671)
1.69 μs ± 31 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

In [23]: FermatFactor(3593 * 3671)
Out[23]: [3593, 3671]

In [24]: %timeit FermatFactor(7187 * 7829)
4.94 μs ± 47.4 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

In [25]: FermatFactor(7187 * 7829)
Out[25]: [7187, 7829]

In [26]: %timeit FermatFactor(8087 * 8089)
1.38 μs ± 12.9 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

In [27]: FermatFactor(8087 * 8089)
Out[27]: [8087, 8089]

Então eu quero usar esse algoritmo para gerar todos os fatores primos de um inteiro qualquer dado (claro que eu sei que isso só funciona com inteiros ímpares, mas isso não é um problema, já que potências de dois podem ser fatoradas trivialmente usando hacking de bits). A maneira mais fácil que eu consigo pensar é chamar recursivamente Fermat_Factoraté que nseja primo. Eu não sei como verificar se um número é primo nesse algoritmo, mas eu notei algo:

In [3]: Fermat_Factor(3)
Out[3]: (1.0, 3.0)

In [4]: Fermat_Factor(5)
Out[4]: (1.0, 3.0)

In [5]: Fermat_Factor(7)
Out[5]: (1.0, 7.0)

In [6]: Fermat_Factor(11)
Out[6]: (1.0, 11.0)

In [7]: Fermat_Factor(13)
Out[7]: (1.0, 13.0)

In [8]: Fermat_Factor(17)
Out[8]: (3.0, 5.0)

In [9]: Fermat_Factor(19)
Out[9]: (1.0, 19.0)

In [10]: Fermat_Factor(23)
Out[10]: (1.0, 23.0)

In [11]: Fermat_Factor(29)
Out[11]: (3.0, 7.0)

In [12]: Fermat_Factor(31)
Out[12]: (1.0, 31.0)

In [13]: Fermat_Factor(37)
Out[13]: (5.0, 7.0)

In [14]: Fermat_Factor(41)
Out[14]: (1.0, 41.0)

O primeiro número na saída deste algoritmo para muitos primos é 1, mas não para todos, como tal, não pode ser usado para determinar quando a recursão deve parar. Aprendi isso da maneira mais difícil.

Então, eu apenas decidi usar a verificação de associação de um conjunto pré-gerado de primos. Naturalmente, isso ocorrerá RecursionError: maximum recursion depth exceededquando a entrada for um primo maior que o máximo do conjunto. Como não tenho memória infinita, isso deve ser considerado um detalhe de implementação.

Então implementei uma versão funcional (para algumas entradas), mas para algumas entradas válidas (produtos de primos dentro do limite), de alguma forma o algoritmo não fornece a saída correta:

import numpy as np
from itertools import cycle

TRIPLE = ((4, 2), (9, 6), (25, 10))
WHEEL = ( 4, 2, 4, 2, 4, 6, 2, 6 )
def prime_sieve(n):
    primes = np.ones(n + 1, dtype=bool)
    primes[:2] = False
    for square, double in TRIPLE:
        primes[square::double] = False
    
    wheel = cycle(WHEEL)
    k = 7
    while (square := k**2) <= n:
        if primes[k]:
            primes[square::2*k] = False
        
        k += next(wheel)
    
    return np.flatnonzero(primes)
    
PRIMES = list(map(int, prime_sieve(1048576)))
PRIME_SET = set(PRIMES)
TEST_LIMIT = PRIMES[-1] ** 2

def FermatFactor(n):
    if n > TEST_LIMIT:
        raise ValueError('Number too large')
    
    if n in PRIME_SET:
        return [n]
    
    a = int(n ** 0.5 + 0.5)
    if a ** 2 == n:
        return FermatFactor(a) + FermatFactor(a)
    
    b2 = abs(a**2 - n)
    while int(b2**0.5) ** 2 != b2:
        a += 1
        b2 = a**2 - n
    
    return FermatFactor(factor := int(a - b2**0.5)) + FermatFactor(n // factor)

Funciona para muitas entradas:

In [18]: FermatFactor(255)
Out[18]: [3, 5, 17]

In [19]: FermatFactor(511)
Out[19]: [7, 73]

In [20]: FermatFactor(441)
Out[20]: [3, 7, 3, 7]

In [21]: FermatFactor(3*5*823)
Out[21]: [3, 5, 823]

In [22]: FermatFactor(37*333667)
Out[22]: [37, 333667]

In [23]: FermatFactor(13 * 37 * 151 * 727 * 3607)
Out[23]: [13, 37, 727, 151, 3607]

Mas não todos:

In [25]: FermatFactor(5 * 53 * 163)
Out[25]: [163, 13, 2, 2, 5]

In [26]: FermatFactor(3*5*73*283)
Out[26]: [17, 3, 7, 3, 283]

In [27]: FermatFactor(3 * 11 * 29 * 71 *  137)
Out[27]: [3, 11, 71, 61, 7, 3, 3]

Por que é esse o caso? Como posso consertar?