[Haskell-cafe] Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)

[Melissa O'Neill]

I'd like to bring the thread that began with a question about ``why  
the "classic sieve" example did so badly when compared to "a C  
implementation of the same thing"'' back to its starting point.  In  
the discussion that ensued, various people have spoken about various  
algorithms to find primes.  I've mentioned my own, and claimed it was  
a more faithful rendition of the Sieve of Eratosthenes in a  
functional setting (where we like to generalize it to produce a lazy  
and infinite list) than the classic one-liner.

But talk is cheap.  What about some actual numbers, and some code for  
some actual implementations...?

I've put together an archive of all the code we've talked about, and  
put it at

    http://www.cs.hmc.edu/~oneill/code/haskell-primes.zip

Here are some performance measurements of the primes algorithms we've  
been discussing, with the task of finding the 1000th prime (7919),  
the 10,000th prime (104729), the 100,000th prime (1299709) and the  
1,000,000th prime (15485863).  The times are with GHC on my PowerMac G5.

   (this table best viewed in a fixed-width font!)
    --------------------------------------------------------
                  Time (s) for Number of Primes
                  ------------------------------------------
    Algorithm     10^3    10^4      10^5    10^6     10^7
    --------------------------------------------------------
    O'Neill (#2)  0.01      0.09     1.45    22.41    393.28
    O'Neill (#1)  0.01      0.14     2.93    47.08    -
    "sieve" (#3)  0.01      0.25     7.28   213.19    -
    Naive         0.32      0.66    16.04   419.22    -
    Runciman      0.02      0.74    29.25   -         -
    Reinke        0.04      1.21    41.00   -         -
    Gale (#1)     0.12     17.99    -       -         -
    "sieve" (#1)  0.16     32.59    -       -         -
    "sieve" (#2)  0.01     32.76    -       -         -
    Gale (#2)     1.36    268.65    -       -         -
    --------------------------------------------------------

Notes:
- The dashes in the table mean "I gave up waiting" (i.e., > 500 seconds)
- "sieve" (#1) is the classic example we're all familiar with
- "sieve" (#2) is the classic example, but sieving a list without  
multiples of 2,3,5, or 7 -- notice how it makes no real difference
- "sieve" (#3) is the classic example, but generating a lazy-but- 
finite list (see below)
- O'Neill (#1) is the algorithm of mine discussed in http:// 
www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf
- O'Neill (#2) is a variant of that algorithm that generates a lazy- 
but-finite list of primes.
- Naive is the the non-sieve-based "divide by every prime up to the  
square root" algorithm for finding primes
- Runciman is Colin Runciman's algorithm, from his _Lazy Wheel Sieves  
and Spirals of Primes_ paper
- Reinke is the ``applyAt'' algorithm Claus Reinke posted here
- Gale (#1) is Yitz Gale's deleteOrd algorithm
- Gale (#2) is Yitz Gale's crossOff algorithm

Returning to where we began, if you want to calculate primes in  
Haskell, there are ways to do it -- ways based on the Sieve of  
Eratosthenes -- where you don't have to feel too bad about the  
performance of functional programs.

     Melissa.

P.S.  Here's the code for a "sieve" (#3) style prime finder.

primesToLimit :: Integer -> [Integer]
primesToLimit limit = 2 : lsieve [3,5..] where
     lsieve ps@(p : xs)
         | p < sqrtlimit = p : lsieve [x | x <- xs, x `mod` p > 0]
         | otherwise     = takeWhile (<= limit) ps
         where
             sqrtlimit = round (sqrt (fromInteger limit))