But does the Fermat Factorization algorithm work?

Question

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?

math.sqrt uses imprecise floating-point math; you should check if math.isqrt(b_square)**2 == b_square, for example. Not immediately sure if that’s the only issue, though. — Ry-
– Ry- ♦, Commented Jan 10, 2024 at 0:21
float.is_integer() is simply not as precise as you're hoping. Past a certain magnitude, all floats are detected as integers, simply because all of the significant digits that they're capable of representing are to the left of the decimal point. — jasonharper
– jasonharper, Commented Jan 10, 2024 at 0:23
another problem is a = math.isqrt(N) + 1 is not equivalent to math.ceil(math.sqrt(N)) which is what Fermat factorization needs to initialize a. Use a = math.isqrt(N) then if a**2 != N: a += 1. Perfect squares have a bug otherwise, e.g. 25 returns (25, 1) instead of (5, 5). — Mark Tolonen
– Mark Tolonen, Commented Jan 10, 2024 at 1:30
Thank you for your bits of advice. Despite having implemented what you said, and even using gmpy2 library instead of math, I find that the algorithm still factorizes 2147483647 incorrectly — lanzik03
– lanzik03, Commented Jan 10, 2024 at 9:52

Ry- · Accepted Answer · 2024-01-10 18:01:08Z

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]

Collectives™ on Stack Overflow

But does the Fermat Factorization algorithm work?

1 Answer 1

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Your Answer

Sign up or log in

Post as a guest

Related