1

I have written this code (in python) for factorize a integer in prime numbers (Fermat's theorem).

#!/usr/bin/python2

import random,math


n=590632926550117049 

a=math.ceil(math.sqrt(n))
b2=a*a-n

while math.sqrt(b2)!=math.floor(math.sqrt(b2)): 
    a=a+1
    b2=a*a-n

b=math.sqrt(b2)

p=a+b
q=a-b

print("p =",p)
print("q =",q)

The number n=590632926550117049 is the product of 57848543*10209987943 but my program returns: 1156469901*510720535. Why ?

EDIT: i.e. with 187 or 15 or another number works fine.

Dwight J. Browne
  • 378
  • 2
  • 9
Ellipticat
  • 198
  • 1
  • 4
  • 11

1 Answers1

2

math.sqrt() uses standard IEEE 64-bit values. It can only calculate accurately for arguments less than ~2**53. Your value for n is larger than that.

If you want exact integer square roots for large numbers, I would recommend gmpy2.

Disclaimer: I maintain gmpy2.

Edit: Here is an updated version of your program.

import gmpy2

n = 590632926550117049

a = gmpy2.isqrt(n) + 1
b2 = a * a - n

while not gmpy2.is_square(b2):
    a = a + 1
    b2 = a * a - n

b = gmpy2.isqrt(b2)

p = a + b
q = a - b

print("p =", p)
print("q =", q)
casevh
  • 11,093
  • 1
  • 24
  • 35
  • I replaced all "math.sqrt" in my code with "gmpy2.sqrt" but the result isn't correct. p and q are wrong. – Ellipticat Dec 29 '14 at 21:44
  • Please use "gmpy2.isqrt" to get an exact integer result. "gmpy2.sqrt" calculates a floating-point result with a default precision of 53 bits. You could increase the precision but using "gmpy2.isqrt" is a better choice. – casevh Dec 29 '14 at 22:02
  • If I use gmpy2.isqrt i have an error when i calculate b2=a*a-n, because a*a is least than n (now i have b2 that is a negative number !). – Ellipticat Dec 29 '14 at 22:19
  • gmpy2.isqrt is equivalent to floor(sqrt). Since you need ceil(sqrt), you will want to add 1. (You should first check that n isn't a perfect square.) I've added modified version of your program to my answer. – casevh Dec 29 '14 at 22:31