fork(18) download
{-# 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
==========================================================================================
-}
Success #stdin #stdout 1.42s 7328KB
stdin
1000000

1000000 15485863

1000000 15485863      --  PrimePi[3935] = 546
 990000 15318907      --  PrimePi[3913] = 541
 950000 14657917 
 900000 13834103 
 800000 12195257 
 700000 10570841 
 650000 9763393 
 600000 8960453 
 400000 5800079 
 250000 3497861 
 200000 2750159 
 78499  1000003 
 50000  611953 
 25000  287117 
 15000  163841 
 10000  104729 
 4000   37813 
stdout
[15485761,15485773,15485783,15485801,15485807,15485837,15485843,15485849,15485857,15485863]