Prime Number generator with recursion and list comprehension

Question

I am a newbie to Haskell programming and have trouble understanding how the below list comprehension expands.

primes = sieve [2..] 
sieve (p:xs) = p : sieve [x | x <-xs, x `mod` p /= 0]

Can someone correct me how the sieve expansion works:

• As we are pattern matching in sieve, p would associate to 2 and xs from [3..].
• Next, in the list comprehension x<-3, but then how can we call sieve with just 3 when there is no short-circuit.

Other thing I do not understand is how recursion works here.

I think it would be clear if one could expand the above one step at a time for first few numbers, say until 5.

Answer 1

Here is an operational description of what sieve does.

To compute sieve (x:xs):

1. Emit the leading element x.
2. From the tail xs, let ys be the list xs with all of the multiples of x removed.
3. To generate the next element, recursively call sieve on ys defined in step 2.

Here is how the first couple of terms are computed:

sieve [2..]
  = sieve (2:[3..])              -- x = 2, xs = [3..]
  = 2 : sieve ys
      where ys = [3..] with all of the multiples of 2 removed
               = [3,5,7,9,...]
  = 2 : sieve [3,5,7,9,...]

and then:

sieve [3,5,7,9,...]              -- x = 3, xs = [5,7,9,11,...]
  = 3 : sieve ys
      where ys = [5,7,9,11,13,15,17,...] with all of the multiples of 3 removed
               = [5,7,  11,13,   17,...]
  = 3 : sieve [5,7,11,13,17,...]

and then:

sieve [5,7,11,13,17,...]         -- x = 5, xs = [7,11,13,17..]
  = 5 : sieve ys
      where ys = [7, 11,13, 17,19,...]  with all of the multiples of 5 removed
               = [7, 11,13, 17,19,...]  (the first one will be 25, then 35,...)
  = 5 : sieve [7,11,13,17,19,...]

etc.

Answer 2

Using an auxiliary function

transform (p:xs) = [x | x <- xs, mod x p /= 0]
                 = filter (\x-> mod x p /= 0) xs  -- remove all multiples of p
                 = xs >>= noMult p                -- feed xs through a tester
-- where
noMult p x = [x | rem x p > 0]      -- keep x if not multiple of p

we can rewrite the sieve function as

       .=================================================+
       |                                                 |
       |  sieve input =                                  |
       |                  .===========================+  |
       |                  |                           |  |
 <~~~~~~~~~~ head input : | sieve (transform input )  |  |
       |                  |                           |  |
       |                  \___________________________|  |
       |                                                 |
       \_________________________________________________|

In an imperative pseudocode, we could write it as

sieve input =
    while (True) :
        yield (head input)           -- wait until it's yanked, and then
        input := transform input     -- advance and loop

This pattern of repeated applications is known as iteration:

iterate f x = loop x
  where
    loop x = x : loop (f x)   -- [x, a:=f x, b:=f a, c:=f b, ...]

So that

sieve xs = map head ( iterate transform xs )

Naturally, head element of each transformed sequence on each step will be a prime, since we've removed all the multiples of the preceding primes on previous steps.

Haskell is lazy, so transformations won't be done in full on each step, far from it -- only as much of the work will be done as needed. That means producing only the first element, and "making a notice" to perform the transformation further when asked:

<--  2 ---<    [2..]
<--  3 ---<    [3..] >>= noMult 2 
<--  5 ---<   ([4..] >>= noMult 2) >>= noMult 3 
<--  7 ---<  (([6..] >>= noMult 2) >>= noMult 3) >>= noMult 5
<-- 11 ---< ((([8..] >>= noMult 2) >>= noMult 3) >>= noMult 5) >>= noMult 7
            ...............

---

Incidentally, this should give us an idea: 3 needn't really be tested by 2; 4..8 needn't really be tested by 3, let alone 5 or 7; 9..24 shouldn't really be tested by 5; etc. What we want is the following:

<-- 2
<-- 3         --<  2(4),    [3..]
<-- 5,7       --<  3(9),    [4..]  >>= noMult 2 
<-- 11,...,23 --<  5(25),  ([9..]  >>= noMult 2) >>= noMult 3 
<-- 29,...,47 --<  7(49), (([25..] >>= noMult 2) >>= noMult 3) >>= noMult 5
                   ......................

i.e. we want the creation of each >>= noMult p filter postponed until p*p is reached in the input.