fork(22) download
copy 
  1. {-# OPTIONS_GHC -O2 -fno-cse #-}
  2. ---------------------------------------------------------------------------------------------
  3. ----- Sieve of Eratosthenes Comparison Table --- Treefold Merge with Wheel by Will Ness -----
  4.  
  5. --- original linear-fold merging idea due to Richard Bird, in M. O'Neill JFP article
  6. --- original tree-like folding idea due to Dave Bayer, on Haskell-cafe
  7. --- double primes feed idea to prevent memoization/space leaks due to Melissa O'Neill
  8. --- simplification of tree-folding formulation and wheel adaptation due to Will Ness
  9.  
  10. --- original Euler sieve one-liner due to Daniel Fischer on haskell-cafe
  11. ---------------------------------------------------------------------------------------------
  12.  
  13. main  = do spec <- (map read.take 2.words) `fmap` getLine
  14.            case spec of
  15.              [n]     -> print $ slice 10 primesTMWE        n    -- 10 up to n-th
  16.              [n,lim] -> print $ slice 10 (primesLimQE lim) n    -- 10 up to n-th, under limit
  17.  
  18. slice n xs ending = take n $ drop (ending-n) xs
  19.  
  20. ---------------- (on using INT: WolframAlpha( PrimePi[2147483647] ) => 105,097,565 :
  21. ----------------                   no chance to get there in 15 sec limit on ideone.com anyway)
  22.  
  23. primesT :: [Int]
  24. primesT = sieve [2..] where
  25.   sieve (p:xs) = p : sieve [x | x <- xs, rem x p /= 0]           -- (* Turner's Trial Division     *)
  26.  
  27. primesE :: [Int]
  28. primesE = sieve [2..] where
  29.   sieve (p:xs) = p : sieve (minus xs [p,p+p..])                   -- (* Sieve of Eratosthenes       *)
  30.  
  31. primesQE :: [Int]
  32. primesQE = 2 : sieve [3,5..] where                                -- (*           ... from squares  *)
  33.   sieve (p:xs)
  34.    | p > 46340 = p : xs                              -- prevent INT wraparound
  35.    | True      = p : sieve (minus xs [p*p,p*p+2*p..])
  36.  
  37. primesEU :: [Int]
  38. primesEU = 2 : sieve [3,5..] where                                -- (* Euler's Sieve               *)
  39.   sieve (p:xs)
  40.    | p > 46340 = p : xs                              -- Prime[4792,..] = 46337, 46349 ...
  41.    | True      = p : sieve (minus xs [p*x | x <- p:xs])
  42. -------------------------------------------------------------------------------------------------------
  43. ----- it is NOT about from *WHERE* to start removing the multiples, BUT **WHEN** to start doing it! ---
  44. -------------------------------------------------------------------------------------------------------
  45. primesPT :: [Int]
  46. primesPT = 2 : sieve [3,5..] 9 (tail primesPT) where
  47.   sieve (x:xs) q ps@ ~(p:t)                                        -- (* Postponed Turner's          *)
  48.    | x < q = x : sieve xs q ps
  49.    | True  =     sieve [x | x <- xs, rem x p /= 0] (head t^2) t
  50. -------------------------------------------------------------------------------------------------------
  51. primesST :: [Int]
  52. primesST = 2 : sieve 3 9 primes_ 0 where
  53.   primes_ = tail primesST                                          -- (* Segmented Turner's sieve    *)
  54.   sieve x q ~(_:t) k =
  55.     let fs = take k primes_
  56.     in  [x | x <- [x,x+2..q-2], and [rem x f /= 0 | f <- fs]]
  57.         ++ sieve (q+2) (head t^2) t (k+1)
  58. -------------------------------------------------------------------------------------------------------
  59. primesPE :: [Int]
  60. primesPE = 2 : sieve [3,5..] 9 (tail primesPE) where
  61.   sieve (x:xs) q ps@ ~(p:t)                                        -- (* Postponed Eratosthenes      *)
  62.    | x < q = x : sieve xs q ps
  63.    | True  =     sieve (xs `minus` [q,q+2*p..]) (head t^2) t
  64. -------------------------------------------------------------------------------------------------------
  65. primesLME :: [Int]
  66. primesLME = 2 : ([3,5..] `minus` fold [[p*p,p*p+2*p..] | p <- primes_])    -- separate feed
  67.  where                                                                      --   for low space
  68.    primes_ = 3 : ([5,7..] `minus` fold [[p*p,p*p+2*p..] | p <- primes_])
  69.    fold  ((x:xs):t)    = x : (xs `union` fold t)                   -- (* Linear Merge Eratosthenes   *)
  70. -------------------------------------------------------------------------------------------------------
  71. primesTME :: [Int]
  72. primesTME = 2 : ([3,5..] `minus` foldt [[p*p,p*p+2*p..] | p <- primes_])
  73.  where
  74.    primes_ = 3 : ([5,7..] `minus` foldt [[p*p,p*p+2*p..] | p <- primes_])
  75.    foldt ((x:xs):t)    = x : union xs (foldt (pairs t))
  76.    pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t               -- (* Treefold Merge Eratosthenes *)
  77. -------------------------------------------------------------------------------------------------------
  78. primesTMWE :: [Int]
  79. primesTMWE = 2:3:5:7: gaps 11 wheel (fold3t $ roll 11 wheel primes_)
  80.  where
  81.    primes_ = 11: gaps 13 (tail wheel) (fold3t $ roll 11 wheel primes_)     -- separate feed
  82.    fold3t ((x:xs): ~(ys:zs:t)) = x : union xs (union ys zs)
  83.                                       `union` fold3t (pairs t)              -- fold3t: 5% ~ 10% speedup
  84.    pairs ((x:xs):ys:t)         = (x : union xs ys) : pairs t
  85.    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:
  86.            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
  87.    gaps k ws@(w:t) cs@ ~(c:u) | k==c  = gaps (k+w) t u              -- (*  better fold, w/ Wheel!   *)
  88.                               | True  = k : gaps (k+w) t cs
  89.    roll k ws@(w:t) ps@ ~(p:u) | k==p  = scanl (\c d->c+p*d) (p*p) ws
  90.                                           : roll (k+w) t u
  91.                               | True  = roll (k+w) t ps
  92. -------------------------------------------------------------------------------------------------------
  93. minus :: [Int] -> [Int] -> [Int]
  94. minus xs@(x:xt) ys@(y:yt) = case compare x y of
  95.     LT -> x : minus xt ys
  96.     EQ ->     minus xt yt
  97.     GT ->     minus xs yt
  98. minus a         b         = a
  99.  
  100. union :: [Int] -> [Int] -> [Int]
  101. union xs@(x:xt) ys@(y:yt) = case compare x y of
  102.     LT -> x : union xt ys
  103.     EQ -> x : union xt yt
  104.     GT -> y : union xs yt
  105. union a         []        = a
  106. union []        b         = b
  107. -------------------------------------------------------------------------------------------------------
  108. ----- *it is about* removing ONLY what MUST be removed, not especially *WHEN* to start doing it -------
  109. -------------------------------------------------------------------------------------------------------
  110.  
  111. primesLimT  :: Int -> [Int]            ----- bounded versions, generating primes UP TO a LIMIT: ------
  112. primesLimT  m = sieve [2..m] where
  113.   sieve (p:xs) = p : sieve [x | x <- xs, rem x p /= 0]           -- (* Turner's / up a to limit    *)
  114.   sieve []     = []
  115.  
  116. primesLimGT :: Int -> [Int]
  117. primesLimGT m = sieve [2..m] where
  118.   sieve (p:xs)
  119.    | p*p > m   = p : xs
  120.    | True      = p : sieve [x | x <- xs, rem x p /= 0]           -- (*           ... guarded       *)
  121.   sieve []     = []
  122.  
  123. primesLimE  :: Int -> [Int]
  124. primesLimE  m = sieve [2..m] where
  125.   sieve (p:xs) = p : sieve (minus xs [p,p+p..m])                  -- (* Sieve of Eratosthenes /lim  *)
  126.   sieve []     = []
  127.  
  128. primesLimGE :: Int -> [Int]
  129. primesLimGE m = sieve [2..m] where                                -- (*           ... guarded       *)
  130.   sieve (p:xs)
  131.    | p*p > m   = p : xs                                -- to STOP EARLY is essential!
  132.    | True      = p : sieve (minus xs [p,p+p..m])
  133.   sieve []     = []
  134.  
  135. primesLimGQE :: Int -> [Int]
  136. primesLimGQE m = 2 : sieve [3,5..m] where                         -- (*            ... from squares *)
  137.   sieve (p:xs)
  138.    | p*p > m   = p : xs
  139.    | True      = p : sieve (minus xs [p*p,p*p+2*p..m]) -- to start from squares is a small improvement
  140.   sieve []     = []
  141.  
  142. primesLimQE :: Int -> [Int]              -- not Guarded:
  143. primesLimQE m = 2 : sieve [3,5..m] where                          -- (*            from squares     *)
  144.   sieve (p:xs)
  145.    | p > mroot = p : sieve (minus xs []) -- no early cut-off, just INT wraparound guarded off
  146.    | True      = p : sieve (minus xs [p*p,p*p+2*p..m]) -- (with ~ same effect though)
  147.   sieve []     = []
  148.   mroot = 46340                 -- was: floor $ sqrt $ fromIntegral m + 1 -- ??
  149.  
  150. primesLimGEU :: Int -> [Int]
  151. primesLimGEU m = 2 : sieve [3,5..m] where                         -- (* Euler's Sieve /lim /guard   *)
  152.   sieve (p:xs)
  153.    | p*p > m   = p : xs                                -- ... prevent INT wraparound too
  154.    | True      = p : sieve (minus xs [p*x | x <- p:xs])
  155.   sieve []     = []
  156.  
  157. {-========================================================================================
  158.           2,000           10,000           15,000             17,000             19,000   primes produced:
  159. ------------------------------------------------------------------------------------------
  160. T     0.08s-3.7MB      2.65s-4.7MB       7.01s-5.8MB        9.52s-5.8MB       12.83s-5.8MB     Turner's
  161.                  n^2.17            n^2.40            n^2.45            n^2.68
  162. LimT                    2.72-4.7         7.06-5.8                                                   /limit
  163. ------------------------------------------------------------------------------------------
  164. E     0.05s-4.8MB      2.10s-5.8MB       5.95s-6.8MB        8.82s-6.8MB       14.00s-7.8MB     Eratosthenes'
  165.                  n^2.30            n^2.60            n^3.14            n^4.15
  166. LimE                    1.17-4.7          3.15-4.7         25k:11.55-5.8                            /limit
  167.                                    n^2.44            n^2.54
  168. 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
  169.                  n^1.97            n^1.28            n^1.12            n^1.03
  170. ------------------------------------------------------------------------------------------
  171. EU    0.41s-48.8MB   4k:1.71s-167.6MB                                                          EUler's
  172.                  n^2.06
  173. LimGEU               4k:0.04- 8.8       15k:0.27-24.2      50k:2.06-166.5                         /limit/guard
  174.                                    n^1.44            n^1.69
  175. ==========================================================================================
  176.  
  177. ==========================================================================================
  178.          78,499           200,000          400,000          1,000,000          2,000,000  primes produced:
  179. ------------------------------------------------------------------------------------------
  180. LimGT   0.65-3.7       2.56s-3.7MB       6.93s- 3.7MB     650k:14.07-3.7                       T/guard /lim
  181.                  n^1.47            n^1.44            n^1.46
  182. PT    0.48s-5.8MB      1.85s-7.8MB       5.08s-11.9MB     800k:13.90-20.1                      Postponed T
  183.                  n^1.44            n^1.46            n^1.45
  184. ST    0.36s-5.8MB      1.38s-7.9MB       3.72s-12.0MB      14.04s-37.6MB                       Segmented T
  185.                  n^1.44            n^1.43            n^1.45
  186. ------------------------------------------------------------------------------------------
  187. LimGE   0.46-4.8       1.62s-4.8MB       4.35s- 4.8MB     900k:14.05-4.8                       E/guard /lim
  188.                  n^1.35            n^1.43            n^1.45
  189. LimQE   0.39-4.8       1.46s-4.8MB       3.88s- 4.8MB      14.81s- 4.8MB                       E/sQ    /lim
  190.                  n^1.41            n^1.41            n^1.46
  191. LimGQE  0.38-4.8       1.35s-4.8MB       3.59s- 4.8MB      13.79s- 4.8MB                       E/sQ /g /lim
  192.                  n^1.36            n^1.41            n^1.47
  193. PE    0.29s-6.8MB      1.12s-11.9MB      3.01s-22.1MB      11.60s-57.0MB                       Postponed E
  194.                  n^1.44            n^1.43            n^1.47
  195. ------------------------------------------------------------------------------------------
  196. LME   0.26s-4.8MB      0.95s- 4.8MB      2.52s- 4.8MB       9.33s- 4.8MB                       Linear Merge
  197.                  n^1.39            n^1.41            n^1.43
  198. TME   0.16s-4.8MB      0.49s- 4.8MB      1.10s- 4.8MB       3.35s- 4.8MB     7.70s- 4.8MB      Tree Merge
  199.                  n^1.20            n^1.17            n^1.22             n^1.20
  200. TMWE  0.07s-4.8MB      0.19s- 4.8MB      0.44s- 4.8MB       1.33s- 4.8MB     3.14s- 4.8MB         with Wheel
  201.                  n^1.07            n^1.21            n^1.21             n^1.24
  202. LimAOE 0.03-4.8MB      0.09s- 4.8MB      0.21s- 5.8MB       0.65s- 7.8MB     2.16s-11.9MB      Arr-odds,"RzqLP"
  203.                  n^1.17            n^1.22            n^1.23             n^1.73
  204.       ====================================================================================
  205.           2m                4m                 6m                7m              8m
  206. TME   7.70 - 4.8      3.4m:14.78-4.8
  207.                  n^1.23
  208. TMWE  3.14s- 4.8MB      7.44s- 4.8MB      12.23s- 4.8MB     14.81s-4.8MB
  209.                  n^1.24         2m->n^1.24         2m->n^1.24
  210.       2.74 - 4.8        6.23 - 4.8        10.14 - 4.8       12.06 -4.8      14.20S- 4.8MB      ONeill PQ-sieve
  211.                  n^1.19         2m->n^1.19         2m->n^1.18        2m->n^1.19                   entry "H8Y5V"
  212.       2.21 - 4.9        5.09 - 4.9         8.30 - 4.9       10.02 -4.9    9m:13.56s-4.9MB      improved ONeill
  213.                  n^1.20         2m->n^1.20         2m->n^1.21        2m->n^1.21                   entry "lKxjW"
  214. LimAOE 2.16-11.9        5.63s-24.2         9.84s-39.6       12.20-36.5      14.62s-40.6MB      Arr-odds,"RzqLP"
  215.                  n^1.38             n^1.38         2m->n^1.38        2m->n^1.38
  216. ==========================================================================================
  217. -}
  218.  
Success #stdin #stdout 1.42s 7328KB
comments (?)
stdin
copy 
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
copy 
[15485761,15485773,15485783,15485801,15485807,15485837,15485843,15485849,15485857,15485863]
https://ideone.com/P0E81
www.haskell.org/haskellwiki/Prime_numbers
gist.github.com/WillNess/7672824#file-p0e81-2013-09-08-hs
gist.github.com/WillNess/7672836#file-p0e81-2013-10-06-hs
language:
Haskell (ghc 8.4.4)
created:
10 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!