Explain this chunk of haskell code that outputs a stream of primes

Contrary to what others have stated here, this function does not implement the true sieve of Eratosthenes. It does returns an initial sequence of the prime numbers though, and in a similar manner, so it's okay to think of it as the sieve of Eratosthenes.

I was about done explaining the code when mipadi posted his answer; I could delete it, but since I spent some time on it, and because our answers are not completely identical, I'll leave it here.

Firs of all, note that there are some syntax errors in the code you posted. The correct code is,

let sieve (p:xs) = p : sieve (filter (\x -> x `mod` p /= 0) xs) in sieve [2..]

let x in y defines x and allows its definition to be used in y. The result of this expression is the result of y. So in this case we define a function sieve and return the result of applying [2..] to sieve.
Now let us have a closer look at the let part of this expression:
```
sieve (p:xs) = p : sieve (filter (\x -> x `mod` p /= 0) xs)
```
1. This defines a function sieve which takes a list as its first argument.
2. (p:xs) is a pattern which matches p with the head of said list and xs with the tail (everything but the head).
3. In general, p : xs is a list whose first element is p. xs is a list containing the remaining elements. Thus, sieve returns the first element of the list it receives.
4. Not look at the remainder of the list:
```
sieve (filter (\x -> x `mod` p /= 0) xs)
```
  1. We can see that sieve is being called recursively. Thus, the filter expression will return a list.
  2. filter takes a filter function and a list. It returns only those elements in the list for which the filter function returns true.
  3. In this case xs is the list being filtered and
```
(\x -> x `mod` p /= 0)
```
    is the filter function.
  4. The filter function takes a single argument, x and returns true iff it is not a multiple of p.
Now that sieve is defined, we pass it [2 .. ], the list of all natural numbers starting at 2. Thus,
1. The number 2 will be returned. All other natural number which are a multiple of 2 will be discarded.
2. The second number is thus 3. It will be returned. All other multiples of 3 will be discarded.
3. Thus the next number will be 5. Etc.

It's actually pretty elegant.

First, we define a function sieve that takes a list of elements:

sieve (p:xs) =

In the body of sieve, we take the head of the list (because we're passing the infinite list [2..], and 2 is defined to be prime) and append it (lazily!) to the result of applying sieve to the rest of the list:

p : sieve (filter (\ x -> x 'mod' p /= 0) xs)

So let's look at the code that does the work on the rest of the list:

sieve (filter (\ x -> x 'mod' p /= 0) xs)

We're applying sieve to the filtered list. Let's break down what the filter part does:

filter (\ x -> x 'mod' p /= 0) xs

filter takes a function and a list on which we apply that function, and retains elements that meet the criteria given by the function. In this case, filter takes an anonymous function:

\ x -> x 'mod' p /= 0

This anonymous function takes one argument, x. It checks the modulus of x against p (the head of the list, every time sieve is called):

 x 'mod' p

If the modulus is not equal to 0:

 x 'mod' p /= 0

Then the element x is kept in the list. If it is equal to 0, it's filtered out. This makes sense: if x is divisible by p, than x is divisible by more than just 1 and itself, and thus it is not prime.

It defines a generator - a stream transformer called "sieve",

Sieve s = 
  while( True ):
      p := s.head
      s := s.tail
      yield p                             -- produce this
      s := Filter (nomultsof p) s         -- go next

primes := Sieve (Nums 2)

which uses a curried form of an anonymous function equivalent to

nomultsof p x = (mod x p) /= 0

Both Sieve and Filter are data-constructing operations with internal state and by-value argument passing semantics.

Here we can see that the most glaring problem of this code is not, repeat not that it uses trial division to filter out the multiples from the working sequence, whereas it could find them out directly, by counting up in increments of p. If we were to replace the former with the latter, the resulting code would still have abysmal run-time complexity.

No, its most glaring problem is that it puts a Filter on top of its working sequence too soon, when it should really do that only after the prime's square is seen in the input. As a result it creates a quadratic number of Filters compared to what's really needed. The chain of Filters it creates is too long, and most of them aren't even needed at all.

The corrected version, with the filter creation postponed until the proper moment, is

Sieve ps s = 
  while( True ):
      x := s.head
      s := s.tail
      yield x                             -- produce this
      p := ps.head
      q := p*p
      while( (s.head) < q ):
          yield (s.head)                  --    and these
          s := s.tail
      ps := ps.tail                       -- go next
      s  := Filter (nomultsof p) s

primes := Sieve primes (Nums 2)

or in Haskell,

primes = sieve primes [2..] 
sieve ps (x:xs) = x : h ++ sieve pt [x | x <- t, rem x p /= 0]
     where (p:pt) = ps
           (h,t)  = span (< p*p) xs

rem is used here instead of mod as it can be much faster in some interpreters, and the numbers are all positive here anyway.

Measuring the local orders of growth of an algorithm by taking its run times t1,t2 at problem-size points n1,n2, as logBase (n2/n1) (t2/t1), we get O(n^2) for the first one, and just above O(n^1.4) for the second (in n primes produced).

Just to clarify it, the missing parts could be defined in this (imaginary) language simply as

Nums x =            -- numbers from x
  while( True ):
      yield x
      x := x+1

Filter pred s =     -- filter a stream by a predicate
  while( True ):
      if pred (s.head) then yield (s.head)
      s := s.tail

Explain this chunk of haskell code that outputs a stream of primes

Related

Recent Posts