May 19, 2009 9:00 AM
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.
Posted by programmingpraxis
Categories: Exercises
Tags:
Mobile Site | Full Site
Get a free blog at WordPress.com Theme: WordPress Mobile Edition by Alex King.
[…] Praxis – Fermat’s Method By Remco Niemeijer Today’s Programming Praxis problem brings us yet another factorization problem, this time using […]
By Programming Praxis – Fermat’s Method « Bonsai Code on May 19, 2009 at 2:42 PM
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 600851475143By Remco Niemeijer on May 19, 2009 at 2:43 PM
Remco 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)))))))))
By programmingpraxis on May 19, 2009 at 3:16 PM
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)
By Jamie Hope on May 25, 2009 at 11:32 PM
Jamie:
This is a well-known limitation of Fermat’s method.
By programmingpraxis on May 26, 2009 at 12:59 PM