Tuesday 22 April 2014

Haskell Prime Number Generator (part 3 of 3)

This is the third of three posts looking at prime number generators in Haskell. The first post looked at an elegant but inefficient one-liner. Part 2 looked at longer but much more efficient solution. The code for those two posts was taken from users gauchopuro and vineet in the Project Euler forum.

In this, the final post on the topic, I tweak vineet's version firstly to make it a bit more readable, and then to improve its efficiency.


A more readable, but still efficient Prime Number Generator

The previously discussed solution uses three separate functions, but using the where keyword and partial application (which are both explained below), this can be reduced to the single function:

primes = 2 : filter (isPrime primes) [3..] where
  isPrime (x:xs) a
    | x*x > a        = True
    | a `mod` x == 0 = False
    | otherwise      = isPrime xs a

The where keyword allows you to use a term that was previously undefined, and which you then define. For example

volumeOfCylinder radius length = length * crossSectionalArea where
  crossSectionalArea = pi * radius * radius

where helps the readability in the above example and it also keeps the variable crossSectionalArea local to volumeOfCylinder. So, for example, adding:

flowRate height width volumetricFlowRate = volumetricFlowRate / crossSectionalArea where
  crossSectionalArea = height * width

to the same Haskell source file is fine and there's no name clash with crossSectionalArea from the original volumeOfCylinder function.


Partial Function Application in Haskell

Partial application is when a function is called without supplying all the arguments needed to evaluate it. For example, if I define the following function:

multiply x y = x * y

then I can define the function double through partial application of multiply like this:

double = multiply 2

In the definition of double I partially apply the multiply function supplying only one argument. The type of double is therefore a function that takes a Num (the second argument to 'multiply 2') and returns a Num. I could have specified the multiply's type to be a function that takes two Nums and returns a Num but if I let Haskell infer its type it will consider it to be a function that takes a Num and returns a function that takes a Num and returns a Num. This is expressed in Haskell notation as:

multiply :: Num a => a -> a -> a

and it makes more sense when considering the possibility of partial application.


Using Partial Application in the Prime Number Generator

The prime number generator at the top of this post contains:

filter (isPrime primes) [3..]

Remember that filter expects a predicate and a list. The predicate is applied to each item in the list and a new list is constructed from those items that return True. isPrime is a predicate, but it needs two arguments and so it can't be used by filter. However, (isPrime primes) is a partially applied function. It's already been supplied with a list of primes and so it only needs a number to test. It is therefore a suitable predicate and can be used by filter.

Along with the details in the previous two posts you should now have enough knowledge to understand this prime number generator.


A Final Efficiency Improvement

The prime number generator at the top of this post is quite clean and also pretty efficient. The final part of this post makes a small tweaks to give a slight efficiency improvement. It does make the code twice as long though, so you may not think it's not worth it.

The improvement is based on the fact that all prime numbers greater than 3 can be written in the form 6k +/-1. This follows from 'sieving' out all the multiples of 2 and 3. If we use this sieved list in the original code in place of [3..] then we'll only be testing prime candidates that can be written as 6k +/-1. This can be achieved by zipping the list of numbers that can be written as 6k-1 together with the numbers that can be written as 6k+1. The final solution therefore looks like this:

primes = 2 : 3 : filter (isPrime primes) times6PlusMinus1 where
  times6PlusMinus1 = concat $ map (\(x,y) -> [x,y]) $ zip times6Minus1 times6Plus1 where
    times6Minus1 = (map (\x -> 6*x-1) [1..])
    times6Plus1 = (map (\x -> 6*x+1) [1..])
  isPrime (x:xs) a
    | x*x > a        = True
    | a `mod` x == 0 = False
    | otherwise      = isPrime xs a

Again, you may not feel that the slight efficiency gain is worth the extra code. You could also think of other ways to zip the lists together, and you could go further and sieve out multiples of more primes. If I come up with or find any major improvements on this I'll do another post but for now I'll finish with this Haskell prime number generator.