Fermat’s Method
May 19, 2009
We have previously studied integer factorization by trial division and wheel factorization. In this exercise we will implement integer factorization by Fermat’s method. Pierre de Fermat was a French lawyer and amateur mathematician of the seventeenth century, a contemporary of René Descartes and Blaise Pascal, who did early work in calculus, number theory, analytic geometry and probability.
Fermat’s method works by noticing that if n, the odd number to be factored, can be written as the difference of two squares n = x2 – y2, then it can be factored as (x – y) × (x + y). Thus, to find the factors of n, we start with y = 0 and x as the smallest integer greater than or equal to the square root of n and increase y until x2 – y2 is equal to n, in which case x and y reveal the factors of n, or x2 – y2 is less than n, when we increase x by one and iterate. This process must terminate; if n is prime, it will stop with the factors 1 and n.
Though this method works, the repeated squarings are relatively expensive, so it is normal to work with u = 2x + 1 and v = 2y + 1 and to keep track of r = x2 – y2 – n; when r = 0, the algorithm terminates. The variable u keeps track of the amount r increases when x2 is replaced by (x + 1)2 and variable v keeps track of the amount r decreases when y is replaced by (y + 1)2; as x and y increase by 1, u and v increase by 2.
Your task is to implement integer factorization by Fermat’s method. When you are finished, you are welcome to read or run the suggested solution, or post your own solution or discuss the exercise in the comments below.
Programming Praxis 
[…] Praxis – Fermat’s Method By Remco Niemeijer Today’s Programming Praxis problem brings us yet another factorization problem, this time using […]
There seems to be a bug in your solution. If I call (display (factors 1)), or any power of 2 (because they all reduce to 1), I get a timeout on codepad.
My Haskell solution (see http://bonsaicode.wordpress.com/2009/05/19/programming-praxis-fermat%E2%80%99s-method/ for a version with comments):
isqrt :: Integer -> Integer isqrt = ceiling . sqrt . fromIntegral factor :: Integer -> [Integer] -> [Integer] factor _ [] = [] factor n (x:xs) | y*y /= x*x - n = factor n xs | x - y == 1 = [x + y] | otherwise = concatMap factors [x - y, x + y] where y = isqrt (x*x - n) factors :: Integer -> [Integer] factors n | even n = 2 : factors (div n 2) | otherwise = factor n [isqrt n..n] main :: IO () main = print $ factors 600851475143Remco pointed out that the suggested solution fails if the number to be factored is 1 or a power of 2. Here’s a corrected version:
(define (factors n)(cond ((= n 1) '())
((even? n) (cons 2 factors (/ n 2))))
(else (let ((s (+ (isqrt n) 1)))
(let loop ((u (+ s s 1)) (v 1) (r (- (* s s) n)))
(cond ((positive? r) (loop u (+ v 2) (- r v)))
((negative? r) (loop (+ u 2) v (+ r u)))
((= (- u v) 2) (list (/ (+ u v -2) 2)))
(else (append (factors (/ (+ u v -2) 2))
(factors (/ (- u v) 2)))))))))
There’s still a bug in your solution if the number to factor is a multiple of a perfect square:
(factors 25)
=> (25)
(factors 36)
=> (2 2 9)
Jamie:
This is a well-known limitation of Fermat’s method.