But does the Fermat Factorization algorithm work?

I was testing the basic version of the Fermat Factorization method, so I'm not using any improvement whatsoever. The code I wrote in Python works as fast as expected for those numbers obtained by relatively close primes.

import math

N = int(input("Insert a large number: "))

def fermat_factorization(N: int) -> list[int]:
    if N % 2 == 0:
        return [2, N // 2]

    a = math.isqrt(N) + 1
    b_square = a**2 - N

    while not math.sqrt(b_square).is_integer():
        a += 1
        b_square = a**2 - N

    b = math.isqrt(b_square)
    return [a+b, a-b]

factors = fermat_factorization(N)
print(f"\nThe factor of {N} are {factors}")

However, upon checking for some large Mersenne primes (e.g. 2147483647), and playing with some more from here, the output missed the marked completely (2147483647 gets factorized into 268435456*8). Is there some error in the code or should I expect the algorithm to fail for large primes?

Solution

math.sqrt uses floating-point math, which isn’t precise enough for the size of integers you’re working with. You can discover this by doing some debugging:

while not math.sqrt(b_square).is_integer():
    a += 1
    b_square = a**2 - N

print(f"Done because {b_square} is a perfect square")

Done because 18014397435740177 is a perfect square

>>> import math
>>> x = math.sqrt(18014397435740177)
>>> x
134217724.0
>>> int(x)**2
18014397435740176
>>> 18014397435740177.0 == 18014397435740176.0
True

Using integer math (with the fix @MarkTolonen suggested for the square N issue) will work:

def fermat_factorization(N: int) -> list[int]:
    if N % 2 == 0:
        return [2, N // 2]

    a = math.isqrt(N)
    if a**2 != N:
        a += 1

    b_square = a**2 - N

    while math.isqrt(b_square)**2 != b_square:
        a += 1
        b_square = a**2 - N

    b = math.isqrt(b_square)
    return [a+b, a-b]

The factor of 2147483647 are [2147483647, 1]