## 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.

Pages: 1 2

### 5 Responses to “Fermat’s Method”

1. […] Praxis – Fermat’s Method By Remco Niemeijer Today’s Programming Praxis problem brings us yet another factorization problem, this time using […]

2. Remco Niemeijer said

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.

```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 600851475143
```
3. programmingpraxis said

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)))))))))```

4. Jamie Hope said

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)

5. programmingpraxis said

Jamie:

This is a well-known limitation of Fermat’s method.