{-# 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
---------------------------------------------------------------------------------------------
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 = sieve [2..] where
sieve
(p:xs
) = p : sieve
[x
| x
<- xs
, rem x p
/= 0] -- (* Turner's Trial Division *)
primesE = sieve [2..] where
sieve (p:xs) = p : sieve (minus xs [p,p+p..]) -- (* Sieve of Eratosthenes *)
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 = 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
= 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 = 2 : sieve 3 9 primes_ 0 where
primes
_ = tail primesST
-- (* Segmented Turner's sieve *) sieve x q ~(_:t) k =
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
= 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 = 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 = 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 = 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 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 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 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 m = sieve [2..m] where
sieve (p:xs) = p : sieve (minus xs [p,p+p..m]) -- (* Sieve of Eratosthenes /lim *)
sieve [] = []
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 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 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
==========================================================================================
-}