Skip to content

Instantly share code, notes, and snippets.

@WillNess

WillNess/p0e81 2013-09-08.hs

Created November 27, 2013 09:10
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save WillNess/7672824 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save WillNess/7672824 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
p0e81 2013-09-08.hs
{-# 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]) -- prevent out-of-bounds (p*p)
sieve [] = []
mroot = 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
==========================================================================================
-}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment