fork(11) download
  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]) -- prevent out-of-bounds (p*p)
  147. sieve [] = []
  148. mroot = 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
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]