{-# OPTIONS_GHC -O2 -fno-cse #-} --------------------------------------------------------------------------------------------- ----- Sieve of Eratosthenes Comparison Table --- Treefold Merge with Wheel by Will Ness ----- --- original linear-fold merging idea due to Richard Bird, in M. O'Neill JFP article --- original tree-like folding idea due to Dave Bayer, on Haskell-cafe --- double primes feed idea to prevent memoization/space leaks due to Melissa O'Neill --- simplification of tree-folding formulation and wheel adaptation due to Will Ness --- original Euler sieve one-liner due to Daniel Fischer on haskell-cafe --------------------------------------------------------------------------------------------- main = do spec <- (map read.take 2.words) `fmap` getLine case spec of [n] -> print $ slice 10 primesTMWE n -- 10 up to n-th [n,lim] -> print $ slice 10 (primesLimQE lim) n -- 10 up to n-th, under limit slice n xs ending = take n $ drop (ending-n) xs ---------------- (on using INT: WolframAlpha( PrimePi[2147483647] ) => 105,097,565 : ---------------- no chance to get there in 15 sec limit on ideone.com anyway) primesT :: [Int] primesT = sieve [2..] where sieve (p:xs) = p : sieve [x | x <- xs, rem x p /= 0] -- (* Turner's Trial Division *) primesE :: [Int] primesE = sieve [2..] where sieve (p:xs) = p : sieve (minus xs [p,p+p..]) -- (* Sieve of Eratosthenes *) primesQE :: [Int] primesQE = 2 : sieve [3,5..] where -- (* ... from squares *) sieve (p:xs) | p > 46340 = p : xs -- prevent INT wraparound | True = p : sieve (minus xs [p*p,p*p+2*p..]) primesEU :: [Int] primesEU = 2 : sieve [3,5..] where -- (* Euler's Sieve *) sieve (p:xs) | p > 46340 = p : xs -- Prime[4792,..] = 46337, 46349 ... | True = p : sieve (minus xs [p*x | x <- p:xs]) ------------------------------------------------------------------------------------------------------- ----- it is NOT about from *WHERE* to start removing the multiples, BUT **WHEN** to start doing it! --- ------------------------------------------------------------------------------------------------------- primesPT :: [Int] primesPT = 2 : sieve [3,5..] 9 (tail primesPT) where sieve (x:xs) q ps@ ~(p:t) -- (* Postponed Turner's *) | x < q = x : sieve xs q ps | True = sieve [x | x <- xs, rem x p /= 0] (head t^2) t ------------------------------------------------------------------------------------------------------- primesST :: [Int] primesST = 2 : sieve 3 9 primes_ 0 where primes_ = tail primesST -- (* Segmented Turner's sieve *) sieve x q ~(_:t) k = let fs = take k primes_ in [x | x <- [x,x+2..q-2], and [rem x f /= 0 | f <- fs]] ++ sieve (q+2) (head t^2) t (k+1) ------------------------------------------------------------------------------------------------------- primesPE :: [Int] primesPE = 2 : sieve [3,5..] 9 (tail primesPE) where sieve (x:xs) q ps@ ~(p:t) -- (* Postponed Eratosthenes *) | x < q = x : sieve xs q ps | True = sieve (xs `minus` [q,q+2*p..]) (head t^2) t ------------------------------------------------------------------------------------------------------- primesLME :: [Int] primesLME = 2 : ([3,5..] `minus` fold [[p*p,p*p+2*p..] | p <- primes_]) -- separate feed where -- for low space primes_ = 3 : ([5,7..] `minus` fold [[p*p,p*p+2*p..] | p <- primes_]) fold ((x:xs):t) = x : (xs `union` fold t) -- (* Linear Merge Eratosthenes *) ------------------------------------------------------------------------------------------------------- primesTME :: [Int] primesTME = 2 : ([3,5..] `minus` foldt [[p*p,p*p+2*p..] | p <- primes_]) where primes_ = 3 : ([5,7..] `minus` foldt [[p*p,p*p+2*p..] | p <- primes_]) foldt ((x:xs):t) = x : union xs (foldt (pairs t)) pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t -- (* Treefold Merge Eratosthenes *) ------------------------------------------------------------------------------------------------------- primesTMWE :: [Int] primesTMWE = 2:3:5:7: gaps 11 wheel (fold3t $ roll 11 wheel primes_) where primes_ = 11: gaps 13 (tail wheel) (fold3t $ roll 11 wheel primes_) -- separate feed fold3t ((x:xs): ~(ys:zs:t)) = x : union xs (union ys zs) `union` fold3t (pairs t) -- fold3t: 5% ~ 10% speedup pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t wheel = 2:4:2:4:6:2:6:4:2:4:6:6:2:6:4:2:6:4:6:8:4:2:4:2: 4:8:6:4:6:2:4:6:2:6:6:4:2:4:6:2:6:4:2:4:2:10:2:10:wheel gaps k ws@(w:t) cs@ ~(c:u) | k==c = gaps (k+w) t u -- (* better fold, w/ Wheel! *) | True = k : gaps (k+w) t cs roll k ws@(w:t) ps@ ~(p:u) | k==p = scanl (\c d->c+p*d) (p*p) ws : roll (k+w) t u | True = roll (k+w) t ps ------------------------------------------------------------------------------------------------------- minus :: [Int] -> [Int] -> [Int] minus xs@(x:xt) ys@(y:yt) = case compare x y of LT -> x : minus xt ys EQ -> minus xt yt GT -> minus xs yt minus a b = a union :: [Int] -> [Int] -> [Int] union xs@(x:xt) ys@(y:yt) = case compare x y of LT -> x : union xt ys EQ -> x : union xt yt GT -> y : union xs yt union a [] = a union [] b = b ------------------------------------------------------------------------------------------------------- ----- *it is about* removing ONLY what MUST be removed, not especially *WHEN* to start doing it ------- ------------------------------------------------------------------------------------------------------- primesLimT :: Int -> [Int] ----- bounded versions, generating primes UP TO a LIMIT: ------ primesLimT m = sieve [2..m] where sieve (p:xs) = p : sieve [x | x <- xs, rem x p /= 0] -- (* Turner's / up a to limit *) sieve [] = [] primesLimGT :: Int -> [Int] primesLimGT m = sieve [2..m] where sieve (p:xs) | p*p > m = p : xs | True = p : sieve [x | x <- xs, rem x p /= 0] -- (* ... guarded *) sieve [] = [] primesLimE :: Int -> [Int] primesLimE m = sieve [2..m] where sieve (p:xs) = p : sieve (minus xs [p,p+p..m]) -- (* Sieve of Eratosthenes /lim *) sieve [] = [] primesLimGE :: Int -> [Int] primesLimGE m = sieve [2..m] where -- (* ... guarded *) sieve (p:xs) | p*p > m = p : xs -- to STOP EARLY is essential! | True = p : sieve (minus xs [p,p+p..m]) sieve [] = [] primesLimGQE :: Int -> [Int] primesLimGQE m = 2 : sieve [3,5..m] where -- (* ... from squares *) sieve (p:xs) | p*p > m = p : xs | True = p : sieve (minus xs [p*p,p*p+2*p..m]) -- to start from squares is a small improvement sieve [] = [] primesLimQE :: Int -> [Int] -- not Guarded: primesLimQE m = 2 : sieve [3,5..m] where -- (* from squares *) sieve (p:xs) | p > mroot = p : sieve (minus xs []) -- no early cut-off, just INT wraparound guarded off | True = p : sieve (minus xs [p*p,p*p+2*p..m]) -- (with ~ same effect though) sieve [] = [] mroot = 46340 -- was: floor $ sqrt $ fromIntegral m + 1 -- ?? primesLimGEU :: Int -> [Int] primesLimGEU m = 2 : sieve [3,5..m] where -- (* Euler's Sieve /lim /guard *) sieve (p:xs) | p*p > m = p : xs -- ... prevent INT wraparound too | True = p : sieve (minus xs [p*x | x <- p:xs]) sieve [] = [] {-======================================================================================== 2,000 10,000 15,000 17,000 19,000 primes produced: ------------------------------------------------------------------------------------------ T 0.08s-3.7MB 2.65s-4.7MB 7.01s-5.8MB 9.52s-5.8MB 12.83s-5.8MB Turner's n^2.17 n^2.40 n^2.45 n^2.68 LimT 2.72-4.7 7.06-5.8 /limit ------------------------------------------------------------------------------------------ E 0.05s-4.8MB 2.10s-5.8MB 5.95s-6.8MB 8.82s-6.8MB 14.00s-7.8MB Eratosthenes' n^2.30 n^2.60 n^3.14 n^4.15 LimE 1.17-4.7 3.15-4.7 25k:11.55-5.8 /limit n^2.44 n^2.54 QE 0.05s-3.7MB 1.19s-4.7MB 2.00s-4.7MB 50k:7.72-4.7 78,498:12.31-4.7 E. from sQuares n^1.97 n^1.28 n^1.12 n^1.03 ------------------------------------------------------------------------------------------ EU 0.41s-48.8MB 4k:1.71s-167.6MB EUler's n^2.06 LimGEU 4k:0.04- 8.8 15k:0.27-24.2 50k:2.06-166.5 /limit/guard n^1.44 n^1.69 ========================================================================================== ========================================================================================== 78,499 200,000 400,000 1,000,000 2,000,000 primes produced: ------------------------------------------------------------------------------------------ LimGT 0.65-3.7 2.56s-3.7MB 6.93s- 3.7MB 650k:14.07-3.7 T/guard /lim n^1.47 n^1.44 n^1.46 PT 0.48s-5.8MB 1.85s-7.8MB 5.08s-11.9MB 800k:13.90-20.1 Postponed T n^1.44 n^1.46 n^1.45 ST 0.36s-5.8MB 1.38s-7.9MB 3.72s-12.0MB 14.04s-37.6MB Segmented T n^1.44 n^1.43 n^1.45 ------------------------------------------------------------------------------------------ LimGE 0.46-4.8 1.62s-4.8MB 4.35s- 4.8MB 900k:14.05-4.8 E/guard /lim n^1.35 n^1.43 n^1.45 LimQE 0.39-4.8 1.46s-4.8MB 3.88s- 4.8MB 14.81s- 4.8MB E/sQ /lim n^1.41 n^1.41 n^1.46 LimGQE 0.38-4.8 1.35s-4.8MB 3.59s- 4.8MB 13.79s- 4.8MB E/sQ /g /lim n^1.36 n^1.41 n^1.47 PE 0.29s-6.8MB 1.12s-11.9MB 3.01s-22.1MB 11.60s-57.0MB Postponed E n^1.44 n^1.43 n^1.47 ------------------------------------------------------------------------------------------ LME 0.26s-4.8MB 0.95s- 4.8MB 2.52s- 4.8MB 9.33s- 4.8MB Linear Merge n^1.39 n^1.41 n^1.43 TME 0.16s-4.8MB 0.49s- 4.8MB 1.10s- 4.8MB 3.35s- 4.8MB 7.70s- 4.8MB Tree Merge n^1.20 n^1.17 n^1.22 n^1.20 TMWE 0.07s-4.8MB 0.19s- 4.8MB 0.44s- 4.8MB 1.33s- 4.8MB 3.14s- 4.8MB with Wheel n^1.07 n^1.21 n^1.21 n^1.24 LimAOE 0.03-4.8MB 0.09s- 4.8MB 0.21s- 5.8MB 0.65s- 7.8MB 2.16s-11.9MB Arr-odds,"RzqLP" n^1.17 n^1.22 n^1.23 n^1.73 ==================================================================================== 2m 4m 6m 7m 8m TME 7.70 - 4.8 3.4m:14.78-4.8 n^1.23 TMWE 3.14s- 4.8MB 7.44s- 4.8MB 12.23s- 4.8MB 14.81s-4.8MB n^1.24 2m->n^1.24 2m->n^1.24 2.74 - 4.8 6.23 - 4.8 10.14 - 4.8 12.06 -4.8 14.20S- 4.8MB ONeill PQ-sieve n^1.19 2m->n^1.19 2m->n^1.18 2m->n^1.19 entry "H8Y5V" 2.21 - 4.9 5.09 - 4.9 8.30 - 4.9 10.02 -4.9 9m:13.56s-4.9MB improved ONeill n^1.20 2m->n^1.20 2m->n^1.21 2m->n^1.21 entry "lKxjW" LimAOE 2.16-11.9 5.63s-24.2 9.84s-39.6 12.20-36.5 14.62s-40.6MB Arr-odds,"RzqLP" n^1.38 n^1.38 2m->n^1.38 2m->n^1.38 ========================================================================================== -}