primes.hs

import Data.List
import Test.QuickCheck
import qualified Data.Map as Map

nats = [1..]
divides a b = b `mod` a == 0
sieve (n:t) m = case Map.lookup n m of
  Nothing -> n : sieve t (Map.insert (n*n) [n] m)
  Just ps -> sieve t (foldl reinsert (Map.delete n m) ps)
  where
    reinsert m p = Map.insertWith (++) (n+p) [p] m
primes = 2 : sieve (map (\n -> n*2+1) nats) Map.empty