Monday, September 12, 2011

Finding Primes in Haskell

I'm currently distracted from life, immersed in Concepts of Modern Mathematics.

Did you know that you can find out whether a number is prime by multiplying all the integers before it, adding one, and then dividing by the number of interest? Me neither! it's like this:

a number p is prime if and only if: (1*2*3*...*(p-1) + 1) mod p = 0

So, if the product of all positive integers less than p is one short of dividing evenly by p, then p is prime. Otherwise, not prime.

As the book points out, this is not a practical way to determine primeness, because if p is 17, that product is fourteen digits long. However, there's a sneaky way to make this more practical.

This formula came out of the concept of arithmetic modulus p, which is a strange sort of arithmetic that considers only the numbers 0 through p-1. Any higher or lower numbers circle around. Any number above p-1 reduces to its remainder when divided by p. It's a bit like when a two-year old is counting: "One, two, three, four, one, two, three..."

Since this equation came out of arithmetic modulus p, we can use this property to shrink the product as we do the multiplication. While we multiply 1 times 2 times 3 times 4, etc, we can stop after each multiplication and mod p. That way the product never gets crazy high -- at most (p-1)*(p-1) -- because at each step it drops to something less than p.

To simplify this equation, we take out the "mod p" because in arithmetic modulus p, "mod p" is free. Also, 0 and p are the same number. (yes, that's weird.)
We can knock a few steps off this procedure:

1*2*...*(p-1) + 1 = 0 = p

1*2*...*(p-1) = p - 1

1*2*...*(p-2) = (p-1)/(p-1) = 1

1*2*3*...*(p-2) = 1

In Haskell, we can write this as: for a list of numbers from 1 to p-2, accumulate a product, modding the product by 17 at each step; see whether the end result is 1.

let isPrime p = foldl (\acc x -> acc * x `mod` p) 1 [2..(p-2)] == 1

This says start with 1, accumulate the product-then-mod, and see whether we come out with 1.

Let's check all the numbers up to 100 and see which ones are prime.

>filter isPrime [2..100]


That's pretty much instantaneous. Finding every prime between 2 and 10000 -- there are 1229 -- takes about thirty seconds on Swirlybob in GHCi. It was twenty seconds on Carolyn, my work laptop.

A future goal: use monads to time the operation precisely. I wonder whether implementing it as a recursive function and explicitly checking for 0 in the accumulator will be faster. (This happens early in the process for all non-primes except perfect squares.)

No comments:

Post a Comment