[Haskell-cafe] Genuine Eratosthenes sieve [Was: Optimization fun]

[oleg at pobox.com]

It has been already remarked that any algorithm of finding prime
numbers that uses division or `mod` operations cannot be called
(Eratosthenes) sieve. The insight of Eratosthenes is finding primes
without resorting to division or multiplication. In his time, doing
either of those operations was quite expensive (especially given the
way numbers were written in Ancient Greece). Thus Eratosthenes devised
one of the early examples of trading space for time.

These days it's the other way around: it is far easier for a modern
CPU to divide numbers than to access memory. So, the original
algorithm is of historical interest it seems. Still, in the interest
of purity, here it is, in Haskell. As the original Eratosthenes sieve,
this algorithm uses only successor and predecessor operations.

-- repl_every_n n l replaces every (n+1)-th element in a list (_:l) with 0
{- The following leaks memory!
repl_every_n :: Int -> [Int] -> [Int]
repl_every_n n l = case splitAt n l of (l1,_:l2) -> l1 ++ 0:(repl_every_n n l2)
-}

-- But this seems to run in constant space...
repl_every_n :: Int -> [Int] -> [Int]
repl_every_n n l = repl_every_n' n l
 where repl_every_n' 0 (_:t) = 0: repl_every_n n t
       repl_every_n' i (h:t) = h: repl_every_n' (pred i) t


primes = 2:(loop [3,5..])
 where loop l = case dropWhile (==0) l of
		  (h:t) -> h:(loop (repl_every_n (pred h) t))

main = putStrLn $ "Last 10 primes less than 100000: " ++ 
       show (take 10 $ reverse  $ takeWhile (< 100000) primes)

{-
Last 10 primes less than 1000: [997,991,983,977,971,967,953,947,941,937]
-}
{-
Last 10 primes less than 10000: [9973,9967,9949,9941,9931,9929,9923,9907,9901,9887]
-}
{-
Last 10 primes less than 100000: [99991,99989,99971,99961,99929,99923,99907,99901,99881,99877]
-}