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.

About these ads

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.

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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 611 other followers

%d bloggers like this: