Programming Praxis – Fermat’s Method

Today’s Programming Praxis problem brings us yet another factorization problem, this time using Fermat’s method. Originally I tried using the algorithm from the assignment, but that resulted in pretty much a direct translation of his Scheme solution, since I couldn’t find a way to express it that was both more elegant and not slower. Since that would be boring, we’re going to use a slightly different approach to Fermat’s method that is slightly shorter and faster.

First we need a way to take the square root of integers:

isqrt :: Integer -> Integer
isqrt = ceiling . sqrt . fromIntegral

What this algorithm does is check whether x2 – n is square for any xs between sqrt(n) and n. If so, we return the factors of n (which are further factorized if needed).

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)

Like the Scheme solution, we divide by 2 until we have an odd number to work with, since Fermat’s method only works on odd numbers.

factors :: Integer -> [Integer]
factors n | even n    = 2 : factors (div n 2)
          | otherwise = factor n [isqrt n..n]

And, as usual, we test our code:

main :: IO ()
main = print $ factors 600851475143

Nine lines of code. Not bad.

Tags: , , , , , ,

3 Responses to “Programming Praxis – Fermat’s Method”

  1. programmingpraxis Says:

    I like mine better. Sorry. My inner loop uses only addition and subtraction. Your inner loop uses squaring and integer square root. That must be less efficient. Not to mention that your integer square root converts from integer to floating point and back; mine uses integers only, and is based on Newton’s approximation.

    By the way, your comment is backwards; Fermat’s method only works on odd numbers (you said even). But your code properly handles odd/even numbers.


  2. Remco Niemeijer Says:

    Thanks for catching the odd/even typo; fixed.

    As for speed: this version was faster than the direct translation of your version for me (ca. 0.9 seconds vs. 1.3 seconds). I would’ve expected your version to be faster as well, but benchmarking proved otherwise.

    As for conceptual purity, we both make sacrifices based on computer architecture: I convert between ints and floats because computers have no infinite-precision floats and you stick to addition and subtraction because it’s faster than multiplication. It’s a toss-up as far as I’m concerned.

    Also, my solution doesn’t get stuck in an infinite loop when trying to factor 1 or a power of 2 (see also my comment on your blog) ;)

  3. programmingpraxis Says:

    Yeah. I’ll fix my bug.

    I’m surprised at the speed difference, also. I’ll do some checking with my Scheme system.

Leave a Reply

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

You are commenting using your 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


Get every new post delivered to your Inbox.

Join 34 other followers

%d bloggers like this: