Primes in Haskell

Question

I'm learning Haskell, and I've tried to generate an infinite list of primes, but I can't understand what my function is doing wrong.

The function:

prime = 2:3:filter (\x -> all (\y -> (mod x y) > 0) (init prime)) [5..]

I think it's the init prime, but the strange thing is that even if I set an upper bound to the range (5..10 for example), the function loops forever and never gets any result for prime !! 2

Can you please tell me what I'm doing wrong?

Answer 1

Well, for one let's look at what init does for a finite list:

init [1] == []
init [1,2] == [1]
init [1,2,3] == [1,2]

ok, so it gives us all but the last element of the list.

So what's init primes? Well, prime without the last element. Hopefully if we implemented prime correctly it shouldn't have a last element (because there are infinitely many primes!), but more importantly we don't quite need to care yet because we don't have the full list for now anyway - we only care about the first couple of elements after all, so for us it's pretty much the same as just prime itself.

Now, looking at all: What does this do? Well, it takes a list and a predicate and tells us if all the elements of the list satisfy the predicate:

all (<5) [1..4] == True
all even [1..4] == False

it even works with infinite lists!

all (<5) [1..] == False

so what's going on here? Well, here's the thing: It does work with infinite lists... but only if we can actually evaluate the list up to the first element of the list that violates the predicate! Let's see if this holds true here:

all (\y -> (mod 5 y) > 0) (init prime)

so to find out if 5 is a prime number, we'd have to check if there's a number in prime minus the last element of prime that divides it. Let's see if we can do that.

Now let's look at the definition of prime, we get

all (\y -> (mod 5 y) > 0) (2:3:filter (\x -> all (\y -> (mod x y) > 0) (init prime)) [5..])

So to determine whether 5 is a prime number, we only have to check if it's:

1. divisible by 2 - it's not, let's continue
2. divisible by 3 - still no
3. divisible by ...? Well, we're in the process of checking what the 3rd prime is so we don't know yet...

and there's the crux of the problem. With this logic, to determine the third prime number you need to know the third prime number! Of course logically, we actually don't want to check this at all, rather we only need to check if any of the smaller prime numbers are divisors of the current candidate.

So how do we go about doing that? Well, we'll have to change our logic unfortunately. One thing we can do is try to remember how many primes we already have, and only take as many as we need for our comparison:

prime = 2 : 3 : morePrimes 2 [5..]
  morePrimes n (x:xs)
    | all (\y -> mod x y > 0) (take n prime) = x : morePrimes (n+1) xs
    | otherwise                              = morePrimes n xs

so how does this work? Well, it basically does what we were just talking about: We remember how many primes we already have (starting at 2 because we know we have at least [2,3] in n. We then check if our next prime is divisible by any of the of n primes we already know by using take n, and if it is we know it's our next prime and we need to increment n - otherwise we just carry on.

There's also the more well known form inspired by (although not quite the same as) the Sieve of Eratosthenes:

prime = sieve [2..] where
  sieve (p:xs) = p : sieve (filter (\x -> mod x p > 0) xs)

so how does this work? Well, again with a similar idea: We know that the next prime number needs to be non-divisible by any previous prime number. So what do we do? Well, starting at 2 we know that the first element in the list is a prime number. We then throw away every number divisible by that prime number using filter. And afterwards, the next item in the list is going to be a prime number again (because we didn't throw it away), so we can repeat the process.

Neither of these are one liners like the one you were hoping for though.

Answer 2

If the code in the other answer is restructured under the identity

[take n primes | n <- [0..]]  ==  inits primes

eventually we get

import Data.List
              -- [ ([], 2), ([2], 3), ([2,3], 5), ... ]
primes = 2 : [ c | (ps, p) <- zip (inits primes) primes,  
                   c <- take 1 [c | c <- [p+1..], 
                                    and [mod c p > 0 | p <- ps]]]

Further improving it algorithmically, it becomes

primes = 2 : [ c | (ps, r:q:_) <- zip (inits primes)                  -- [] [3,4,...]
                                      (tails $ 3 : map (^2) primes),  -- [2] [4,9,...]
                   c <- [r..q-1], and [mod c p > 0 | p <- ps]]        -- [2,3] [9,25,...]