Genuine Eratosthenes sieve [Was: Optimization fun]

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] -}

oleg@pobox.com wrote:
Still, in the interest of purity, here it is, in Haskell. As the original Eratosthenes sieve, this algorithm uses only successor and predecessor operations.
I don't think the Greeks had too much trouble with addition. If that is the case, then Rafael's definition is still valid after a slight modification, and still the clearest: \begin{code} -- Delete the elements of the first list from the second list, -- where both lists are assumed sorted and without repetition. deleteOrd :: Ord a => [a] -> [a] -> [a] deleteOrd xs@(x:xs') ys@(y:ys') | x > y = y : deleteOrd xs ys' | x < y = deleteOrd xs' ys | otherwise = deleteOrd xs' ys' deleteOrd _ ys = ys sieve (x:xs) = x : sieve (deleteOrd [x+x,x+x+x..] xs) sieve _ = [] primes = sieve [2..] \end{code} Regards, Yitz
participants (2)
-
oleg@pobox.com
-
Yitzchak Gale