Saturday, January 5, 2013

Sieve of Eratosthenes in F#

"God may have not play dice with the universe, but something strange is going on with the prime numbers."
Paul Erdős

"Mathematicians  have tried in vain to this day to discover some order in the sequence of prime numbers, and we have reason to believe that it is a mystery into which the human mind will never penetrate."
Leonhard Euler

Before Caesar, before one anno domini, before the internet, there was the library of Alexander and  Eratosthenes.  Eratosthenes of Cyrene was the third librarian of the library of Alexander.  Among other things, he was the first person to calculate the circumference and tilt of the Earth!

If that was not enough, Eratosthenes summoned a simple prime sieve to find prime numbers.  The Sieve of Eratosthenes works in the following way:
  1. list integers from 2 to N
  2. initially mark 2 as prime (it is prime)
  3. cross off all of the integers which are increments of the marked number apart
  4. take the next non-crossed off number and mark it as prime
  5. repeat step 3 until all numbers are either crossed off or marked as prime

Visually here is what the sieve looks like:





skip to the end



We mark the prime numbers with green and those that have been crossed off with red.

What would the Sieve of Eratosthenes look like in F#?  Well we would want to have a list created from 2 to N, then a function which would take the head of the list and filter out the numbers from the list that are divisible by the head (I know this is not exactly the same as what is above) and loop back around to the top of the function with the new filtered list.


> let sieve n =
-   let rec aux list =      
-     match list with
-     | [] -> []
-     | hd :: tl -> hd :: aux (list |> List.filter (fun x -> x % hd <> 0) )
-   aux [2 .. n]
- ;;

val sieve : int -> int list


Let's test it (I'll blog about using FsUnit at a latter date, so we'll just use poking around testing for now).


> sieve 3;;
val it : int list = [2; 3]
> sieve 11;;
val it : int list = [2; 3; 5; 7; 11]
> sieve 100;;
val it : int list =
  [2; 3; 5; 7; 11; 13; 17; 19; 23; 29; 31; 37; 41; 43; 47; 53; 59; 61; 67; 71;
   73; 79; 83; 89; 97]
> sieve 2;;  
val it : int list = [2]
> sieve 1;;
val it : int list = []

Looks like it works.