pythonprimesfactorization

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]