Monday, August 4, 2008

F# Sieve of Eratosthenes

My neighbor obsesses over rare candies. He'll travel to small towns with 5 and Dimes to buy the the wax lips and Mallo Cups he remembers growing up with. White people famously obsess over accumulating pirated music. I think it represents modern man's search for identity in a post-industrial society where we now have the luxury to contemplate more than just feeding our families and surviving the winter. It's intrinsically a play activity normally abandoned in young adulthood, but the affluence of our society has afforded us to become a culture of overgrown man-childs still tinkering with childhood past-times. So what do you obsessively accumulate? Oh, you're above it all, I understand. You're content with adult, intellectual pursuits... like spending Summer plowing through an endless list of fiction and memoirs (an accumulation of vague insights and good party conversation?). Perhaps you have a spiritual escape from it, and have melded a Christian upbringing with Kabbalah/ Scientology/ Pastafarianism) into a compassionate and pious life. Just be careful to not fall into Spirtual Materialism, Chögyam Trungpa's phrase describing the accumulation of spiritual detritus that simply serves to build up the ego.

If this sounds negative then you've gotten the wrong idea... I LOVE accumulating crap. Coins, comics, and crummy implementations of toy algorithms. Collecting useless junk is my birthright as a child of the modern age. My favorite: The Sieve of Eratosthenes. There's so much joy in retreading a useless relic from the past. Without further ado, here goes!

Maybe you remember Eratosthenes, the 3rd Century Greek who gave us the first accurate estimate of the Earth's circumference in the 3rd century bc. Perhaps not. He also gave us a method for finding prime numbers which, up until the 1970s, was the best way to find them. The algorithm from Wikipedia is a good starting point. They have you start with a finite list of numbers from 2 (2..100, say). The first element of the list is always prime, so 2 is prime. Now remove all multiples of 2 from the list, leaving only odd numbers. Again, the first element of the remaining list is always prime, so 3 is prime. Now remove all multiples of 3 from the list. First element (5) is prime, remove multiples of 5 from list, etc. etc. Implementing this algorithm will leave you with a list where the is always a prime number and the list is filtered each time the head is taken.

This method certainly works, and is the basis for Chris Smith's F# Sieve of Eratosthenes. Personally, I like mine better because all the data is immutable, there is no state, and it uses only built-in F# types (but kudos to Chris for getting the ball rolling):
#light

let primes max =
let rec next potentialPrimes =
match potentialPrimes with
| [] -> []
| n :: tail -> n :: next (List.filter (fun x -> x % n <> 0) tail)

next [ 2 .. max]

printfn "%A" (primes 1000000)
Reading through this quickly: "primes" is a function that returns all the prime numbers under the specified maximum. It is implemented as an inner, recursive function called "next", which is initailized with a list of all natural numbers from 2 to the maximum. "next" takes the list it is given as a parameter and a) returns an empty list when given an empty list, or b) concatenates the head of the list (remember, the head is always prime) with the result of a recursive call to next passing in the tail of the list with multiples of the head value filtered out. Writing this was every bit as satisfying as wax lips on a summer day. Mmm.

So why not end it here?

The problem with all this is that the Wikipedia algorithm starts with a finite list of numbers. What if you wanted to calculate all the prime numbers until you ran out of resource? What would you initialize the "next" function with? Long.MAX_VALUE? That isn't high enough if you have dedicated, purpose built hardware (which they had in the 1970s in order to run this exact algorithm). Infinity? The first step of this algorithm is to create an in-memory list of all numbers from 2 to max... if you set max to infinity then you'd run out of resources building the initial list before the first prime was even pulled. Ouch. Hopefully that bug is caught before it makes it to production.

A much cooler approach is to use infinite streams. And F# provides the perfect starting point: Seq. An infinite stream is a data structure that looks and behaves like a List type, but produces elements on demand using a function. The Seq.init_sequence function creates a seq instance that produces elements from the offset requests. So all natural numbers can be expressed lazily and infinitely using an identity function like so:

#light

// creates a Sequence of all numbers
let allNumbers = Seq.init_infinite (fun i -> i)
So what we need is a sieve based on lazy evaluation, bigint and Seq types. I humbly submit the following:
#light

let sieve max =

// creates a Sequence based on a "multiple of" function
let multiplesOf n = Seq.init_infinite (fun i ->
bigint.FromInt32(i) * n)

// checks if a given infinite bigint
// sequence contains a given bigint
let seqContains x seq =
let rec innerSeqContains count =
match x with
| _ when Seq.nth count seq = x -> true
| _ when Seq.nth count seq > x -> false
| _ -> innerSeqContains (count+1)
innerSeqContains 0

// checks if a list of infinite bigint
// sequences contains a given bigint
let rec seqsContain(x, seqs:seq<bigint> list) =
List.exists (
fun s -> seqContains x s) seqs

// generates a list of prime numbers, accumulating all
// non-prime numbers as a list of infinite sequences
let rec innerSieve (count : bigint) nonPrimes =
match count with
| _ when count = max -> []
| _ when count = 1I -> 1I :: innerSieve (count+1I) nonPrimes
| _ when seqsContain (count, nonPrimes) ->
innerSieve (count+1I) nonPrimes
| _ ->
let nextSeq = multiplesOf count
let newNonPrimeList = nextSeq :: nonPrimes
count :: innerSieve (count+1I) newNonPrimeList
innerSieve 2I []

printfn "%A" (sieve 500I)
Some details:

...multiplesOf is a helper function to create a Seq representing all multiples of a given bigint.
...seqContains determines if a given seq contains said bigint. The builtin contains can't be used without an infinite loop.
...seqsContain determines if any element of a seq of seq of bigint contains said bigint
...innerSieve simply produces the integers (this could in itself be expressed as a seq)

So what's up the initializing this thing to 500? Hey, I don't have specialized hardware, how do you expect me to test this thing? Anything over a couple thousand primes and my machine starts to crawl.

The sieve is a cool trick, but since the 70s it has been superseded by replaced as a prime finder by probabilistic techniques.

4 comments:

Anonymous said...

Firstly, thanks for posting these code snippets. They've been of much help to an F# newbie such as myself. I was however just testing your first (primes) function and it appears that your implementation of the Sieve of Eratosthenes is only half complete. The algorithm is sieving using all prime numbers up to the max when it in fact only needs to go up to sqrt(max), which clearly greatly reduces the time complexity.

Here's the modified code in case you or anyone else is interested.

/// Generates prime numbers using Sieve of Eratosthenes.
let primes max =
let sqrt_max = int64 (sqrt (float max))
let rec next potentialPrimes =
match potentialPrimes with
| [] -> []
| n :: tail when n > sqrt_max -> n :: tail
| n :: tail -> n :: next (tail |> List.filter (fun x -> x % n <> 0L))
next [2L .. max]

I expect your machine shouldn't be crawling until millions rather than thousands of primes now!

try2bGood said...

An interesting study!

But...

First, this isn't really the Sieve of Eratosthenes even though the algorithm is often erroneously called that, including Wikipedia at http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes. Melissa E. O'Neill at http://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf properly identifies the algorithm for what it is: an "unfaithful" incompletely constructed sieve by trial division (unfaithful in that it doesn't eliminate the squares as pointed out by Alex). A proper Sieve of Eratosthenes cancels composite numbers by a pattern shift operation using additions and not modulus or multiplication operations and uses other tricks in order to minimize the number of operations.

Secondly, your original implementation of the trial division algorithm and Alex's improved version that respects stopping at the square root of the maximum test number are not really in the spirit of functional declarative programming in that they don't work for an infinite series, although neither would the "pure" sequence based solution as it would have time constraints as you've seen.

Here is a very much more concise lazy evaluation "pure" seq based version of the true trial division as per O'Neil:

let rec tdsieve possibles =
seq {
let h = Seq.hd possibles
yield h
yield!
(tdsieve
((possibles |> Seq.skip 1) |> Seq.filter (fun i ->
if i < (h * h) then
true
else
(i % h) > 0))) }
let tdprimes =
Seq.append
(seq {2..2})
(tdsieve
(Seq.init_infinite
(fun i -> (i + i + 3))))

Note that the "if" clause in the filter lambda function compensates for starting the composite compensations at the square of the found prime although it does this in reverse of the imperative version that uses the square root. The difference is that these functions will be applied in the future when the lazy expression is actually evaluated. This explains why the algorithm is still so slow in that all kinds of nested future computations are set up that aren't used until the iteration reaches those high number and the late evaluation kicks in.

This program can calculate your first 500 primes in about 35 seconds with only a one second saving due from the extra square test due to these still being quite low prime values, as predicted by O'Neill.

Alex's program can calculate about 300,000 primes in about 30 seconds.

Note as well that this is the usual functional declarative lazy evaluation version with no requirement of a maximum setting and would theoretically work up to infinity. However, practically it is still quite slow taking about 30 seconds to compute the first 500 values unlike the more imperative solutions. However, if this were the true Sieve of Eratosthenese, it would be "a lot" to "quite a lot" faster, depending on the algorithm, even if written using declarative function style.

To see a F# version of the real Sieve of Eratosthenese, refer to Dustin Campell's blog at http://diditwith.net/2009/01/20/YAPESProblemSevenPart2.aspx, which calculates 300,000 primes in about 32 seconds and even that could be improved with further algorithm tweaks. So the pure lazy evaluation method comes in at about the same speed as the more imperative solution with the right algorithm.

It's all in the algorithm!

Hamlet D'Arcy said...

Thanks for the comment, very informative. I did see O'Neill's Sieve paper thru LtU, but you explained it in much better terms (more accessible that is).

I'm amazed how much I still have to learn from simply toying with the Sieve!

try2bGood said...

Since you are a collector of Sieve programs, perhaps you will be interested in the following direct translation to F# of the Haskell code that is so often cited (ie. Wikipedia) as being such a good persuader of why functional languages are so elegant but in fact is one of the worst uses of code there could possibly be as to its purpose of filtering prime numbers. In the often used "all on two lines format (it will wrap but you'll get the idea and can extract it):

let rec sieve p = seq {let hd = Seq.hd p in yield hd; yield! (p |> Seq.skip 1 |> Seq.filter (fun x -> x % hd > 0) |> sieve)}

let primes = Seq.init_infinite (fun i -> i + 2) |> sieve

In words the algorithm is: Find the sequence of all primes by passing the infinite sequence of possible prime numbers which include all integers greater than 1 to a sieve filter. The sieve filter yields the first number in the sequence as a known prime and filters the rest of the sequence by eliminating all that can be evenly divided by the just found prime. The resulting sequence is passed recursively to the sieve filter.

This is an elegant expression of what prime numbers are by definition, but is a terrible implementation of code for two reasons:

1) Implementation: the code can not easily be optimized to a loop by tail recursion and thus puts heavy loads on the stack such that the program fails with default stack sizes at something between 400 and 500 primes found. Even to find the 600th prime in the sequence took about 36 seconds on my fairly fast machine due to all the recursive nesting in the lazy evaluation of sequences, which is about what your more complex program takes. This is pretty much useless for any practical purpose, as well as not being a Sieve of Eratosthenese as explained by O'Neill.

2) Algorithm: This is the very worst possible way to find primes even using trial division in that none of the known properties of primes are respected. First, two is the only even prime number so there is no need to scan other even numbers higher than it. Second, it is obvious that one does not have to check for primes higher than the root of the test value as other composite numbers would have already been eliminated by smaller primes.

The results of translating O'Neill's improved Haskell algorithm into F#, including the above optimizations, is as follows (divide and indent at the end of the '=' and 'in' to make 4 lines for #light syntax):

let primes =
let rec mprimes dummy =
let isprime x = mprimes dummy |> Seq.take_while (fun p -> p * p <= x) |> Seq.for_all (fun p -> x % p > 0) in
Seq.append {2..2} (Seq.init_infinite (fun i -> i + i + 3) |> Seq.filter isprime) in
mprimes None


Note that the callees are sub functions of the main prime function because F# does not allow forward declarations as required by this recursion. Also, to avoid a #nowarn "40" compiler metacommand, I've nested a "mprimes" function with a dummy variable so that compiler can determine the "mprimes" is a function rather than an object for recursion.

This new algorithm and program can find 10000 primes in about 34 seconds on the same machine, and it still isn't a true Sieve of Eratosthenes!

It's all in the algorithm!!!