Fermat’s Method
May 19, 2009
The basic method described above finds two factors, but neither is necessarily prime, so we augment the method to call itself recursively until it proves that all the factors are prime:
(define (factors n)
(if (even? n) (cons 2 (factors (/ n 2)))
(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)))))))))
The factors are returned in the order they are found, which is not likely sorted; if you want them in order, you’ll have to sort them yourself:
> (factors 600851475143)
(1471 839 6857 71)
Since the inner loop has only addition and subtraction, it runs quickly, but since the number of cycles required is large, the algorithm as a whole runs rather slowly, essentially in time relative to the square root of n.
Factors uses isqrt from the Standard Prelude. You can run the code at http://codepad.org/Exr2kj9L.
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.