fork download
copy 
  1. {-# OPTIONS_GHC -O2 -fno-cse #-}
  2. {-# LANGUAGE BangPatterns #-}
  3.  
  4. ---------------------------------------------------------------------------------------------
  5. --- Sieve of Eratosthenes Comparison Table --- Treefold Merge with Wheel --- by Will Ness ---    *)
  6. ---------------------------------------------------------------------------------------------
  7.  
  8. import Array
  9.  
  10. main0 = 
  11.  let pairs (x:y:t) = [x,y]:pairs t
  12.  in  
  13.     mapM_ print $ take 50 $ pairs $ zip [1..] $ zip primesTMWE primesFPQW    -- -fno-cse 
  14.  
  15.  
  16. main  = do spec <- (map read.take 2.words) `fmap` getLine  
  17.            case spec of  
  18.              [n]     -> print $ slice 10 primesTMWE         n    -- 10 up to n-th
  19.              [n,lim] -> print $ slice 10 (primesLimQE lim) n    -- 10 up to n-th, under limit
  20.  
  21. slice n xs ending = take n $ drop (ending-n) xs
  22.  
  23. ----- (on using INT: WolframAlpha( PrimePi[2147483647] ) => 105,097,565 :       
  24. -----                 no chance to get there in 15 sec limit on ideone.com anyway)
  25.  
  26. primesLimGQAE :: Int -> [Int]              
  27. primesLimGQAE m = sieve 3 $
  28.                     accumArray (\b c->False) True (3,m) 
  29.                                 [(i,()) | i<-[4,6..m]]
  30.  where                                                            -- (*  ... on array, from squares *)  
  31.   sieve !p !a 
  32.    | p*p > m   = 2 : [i | (i,True) <- assocs a]
  33.    | a!p       = sieve (p+2) (accum (\b c->False) a [(i,()) | i <- [p*p,p*p+2*p..m]])
  34.    | otherwise = sieve (p+2) a
  35.  
  36. {-   well, this is disappointing. I hoped for an array update ((//) a b)
  37.      to take up time ~ O(|b|), to give it TRUE complexity of a Sieve of E.,
  38.      but it is apparently ~ O(|a|) and so gives it even WORSE 
  39.      complexity than that of simple listMinus versions (where length
  40.      of first-arg list diminishes w/ iterations, and here it always
  41.      stays the same, for the same array) ... so to update ONE PLACE in array
  42.      has a cost of LENGTH_OF(ARRAY) ?????????? what is this????
  43.      ... well of course the documentation says "makes an altered copy of ..."
  44.      but since it is passed along to further iterations, one would expect
  45.      the implementation to recognize the fact that the previous one
  46.      EXISTS NO MORE and just UPDATE the dang ***UNIQUE*** thing!
  47.      and this is  not at all unreasonable: making a copy of an array 
  48.      should be an O(1) operation after all!!!!!!!!!!! not O(|a|) in any case.
  49.  
  50.      ... will bangPatterns help ? ... no, it just runs slower...
  51.  
  52.  
  53. here's the working Array version, which needs to be segmented
  54. into the lmited-sized chunks of array betwen the primes squares
  55. because as it is, its memory usage grows rapidly:
  56.  
  57.  
  58. primes () = 2: primes'   -- marking off composites on Array segments
  59.  where                   -- between consecutive primes squares, on odds
  60.    primes' = 3: sieve primes' 3 []
  61.    sieve (p:ps) x fs = 
  62.      let 
  63.          q = (p*p-x)`div`2                  
  64.          a = accumArray (\b c-> False) True 
  65.                  (1,q-1)
  66.                  [(i,()) | (s,y) <- fs, i <- [y+s,y+s+s..q]]
  67.      in  
  68.         [i*2 + x | (i,e) <- assocs a, e]
  69.         ++ sieve ps (p*p)
  70.                  ((p,0):[(s,rem (y-q) s) | (s,y) <- fs])
  71. -}
  72.  
  73.  
  74.  
  75. primesT :: [Int]            
  76. primesT = sieve [2..] where
  77.   sieve (p:xs) = p : sieve [x | x <- xs, rem x p /= 0]           -- (* Turner's Trial Division     *)
  78.  
  79. primesE :: [Int] 
  80. primesE = sieve [2..] where
  81.   sieve (p:xs) = p : sieve (minus xs [p,p+p..])                   -- (* Sieve of Eratosthenes       *)  
  82.  
  83. primesQE :: [Int]              
  84. primesQE = 2 : sieve [3,5..] where                                -- (*           ... from squares  *)  
  85.   sieve (p:xs) 
  86.    | p > 46340 = p : xs                              -- prevent INT wraparound 
  87.    | True      = p : sieve (minus xs [p*p,p*p+2*p..])
  88.  
  89. primesEU :: [Int] 
  90. primesEU = 2 : sieve [3,5..] where                                -- (* Euler's Sieve               *) 
  91.   sieve (p:xs) 
  92.    | p > 46340 = p : xs                              -- Prime[4792,..] = 46337, 46349 ...
  93.    | True      = p : sieve (minus xs [p*x | x <- p:xs]) 
  94. -------------------------------------------------------------------------------------------------------
  95. ----- it is NOT about *WHERE* to start removing the multiples, BUT **WHEN** to start doing it! ... ----
  96. -------------------------------------------------------------------------------------------------------
  97. primesPT :: [Int]                             
  98. primesPT = 2 : sieve [3,5..] (tail primesPT) 9
  99.  where
  100.    sieve (x:xs) ps@ ~(p:ps') q                                     -- (* Postponed Turner's          *)
  101.     | x < q = x : sieve xs ps q
  102.     | True  =     sieve [x | x <- xs, rem x p /= 0] ps' (head ps'^2)
  103. -------------------------------------------------------------------------------------------------------
  104. primesST :: [Int]
  105. primesST = 2 : primes' 
  106.   where                                                            -- (* Segmented Turner's sieve    *)
  107.     primes' = 3 : sieve primes' 3 0
  108.     sieve (p:ps) x k = let pfx = take k primes' in
  109.        [n | n <- [x+2,x+4..p*p-2], and [rem n f /= 0 | f <- pfx]]
  110.        ++ sieve ps (p*p) (k+1)
  111. -------------------------------------------------------------------------------------------------------
  112. primesPE :: [Int]                             
  113. primesPE = 2 : sieve [3,5..] (tail primesPE) 9
  114.  where
  115.    sieve (x:xs) ps@ ~(p:ps') q                                     -- (* Postponed Eratosthenes      *)
  116.     | x < q = x : sieve xs ps q
  117.     | True  =     sieve (xs `minus` [q,q+2*p..]) ps' (head ps'^2)
  118. -------------------------------------------------------------------------------------------------------
  119. primesLME :: [Int]                             
  120. primesLME = 2 : ([3,5..] `minus` fold [[p*p,p*p+2*p..] | p <- primes'])  -- separate feed 2save space
  121.  where
  122.    primes' = 3 : ([5,7..] `minus` fold [[p*p,p*p+2*p..] | p <- primes'])
  123.    fold  ((x:xs):t)    = x : (xs `union` fold t)                   -- (* Linear Merge Eratosthenes   *)
  124. -------------------------------------------------------------------------------------------------------
  125. primesTME :: [Int]                             
  126. primesTME = 2 : ([3,5..] `minus` foldt [[p*p,p*p+2*p..] | p <- primes']) -- iterate (+ (2*p)) (p*p)
  127.  where
  128.    primes' = 3 : ([5,7..] `minus` foldt [[p*p,p*p+2*p..] | p <- primes']) -- enumFromThen (p*p) (p*p+2*p)
  129.    foldt ((x:xs):t)    = x : union xs (foldt (pairs t))
  130.    pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t               -- (* Treefold Merge Eratosthenes *)
  131. -------------------------------------------------------------------------------------------------------
  132. primesTMWE :: [Int]
  133. primesTMWE = 2:3:5:7: gaps 11 wheel (fold3t $ roll 11 wheel primes')
  134.  where                                                             
  135.    primes' = 11: gaps 13 (tail wheel) (fold3t $ roll 11 (2:tail wheel) primes') -- for -fno-cse
  136.    fold3t ((x:xs): ~(ys:zs:t)) = x : union xs (union ys zs)           
  137.                                       `union` fold3t (pairs t)     --  fold3t: 5% speedup vs foldt
  138.    pairs ((x:xs):ys:t)         = (x : union xs ys) : pairs t
  139.    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:  
  140.            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 
  141.    gaps k ws@(w:t) cs@ ~(c:u) | k==c  = gaps (k+w) t u             -- (*  better fold, w/ Wheel!    *)
  142.                               | True  = k : gaps (k+w) t cs  
  143.    roll k ws@(w:t) ps@ ~(p:u) | k==p  = scanl (\c d->c+p*d) (p*p) ws 
  144.                                           : roll (k+w) t u              
  145.                               | True  = roll (k+w) t ps                 
  146.  
  147.  
  148. primesAME :: [Int]                                                 -- (* Array Merge Eratosthenes *)
  149. primesAME = 2 : ([3,5..] `minus` ajoin [[p*p,p*p+2*p..] | p <- primes'])
  150.  where
  151.    primes' = 3 : ([5,7..] `minus` ajoin [[p*p,p*p+2*p..] | p <- primes'])
  152.    ajoin ((x:xs):t) = [x] -- x : go [xs] t
  153.    m = 10000
  154.    -- go sts ((x:xs):t) = 
  155.  
  156.  
  157.  
  158.  
  159.  
  160. {-
  161.  
  162. -- reified tree
  163.  
  164. data T = N {c::Int, p::Int, ws::[Int]} | L {p::Int, ws::[Int]}
  165.  
  166. primesTMWRE :: [Int]
  167. primesTMWRE = 2:3:5:7: gaps 11 wheel (fold3t mults)
  168.  where                                                             --  double feed, to fix space leak
  169.    primes' = 11: gaps 13 (tail wheel) (fold3t mults) 
  170.    mults = roll 11 wheel primes'
  171.    fold3t (L p ws: ~(b:c:t)) = x : union xs (union ys zs)           
  172.                                       `union` fold3t (pairs t) 
  173.    pairs ((x:xs):ys:t)         = (x : union xs ys) : pairs t 
  174.  
  175.    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:  
  176.            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 
  177.    gaps k ws@(w:t) cs@ ~(c:u) | k==c  = gaps (k+w) t u             -- (*  fold3, 2x feed, w/ Wheel!  *)
  178.                               | True  = k : gaps (k+w) t cs  
  179.    roll k ws@(w:t) ps@ ~(p:u) | k==p  = L p ws -- scanl (\c d->c+p*d) (p*p) ws 
  180.                                           : roll (k+w) t u              
  181.                               | True  = roll (k+w) t ps                 
  182. -}
  183. -------------------------------------------------------------------------------------------------------
  184. minus :: [Int] -> [Int] -> [Int]
  185. minus xs@(x:xt) ys@(y:yt) = case compare x y of
  186.     LT -> x : minus xt ys
  187.     EQ ->     minus xt yt
  188.     GT ->     minus xs yt
  189. minus a         b         = a
  190.  
  191. union :: [Int] -> [Int] -> [Int]
  192. union xs@(x:xt) ys@(y:yt) = case compare x y of
  193.     LT -> x : union xt ys
  194.     EQ -> x : union xt yt
  195.     GT -> y : union xs yt
  196. union a         []        = a
  197. union []        b         = b
  198. -------------------------------------------------------------------------------------------------------
  199. ----- *it is about* removing ONLY what MUST be removed, not especially *WHEN* to start doing it -------
  200. -------------------------------------------------------------------------------------------------------
  201.  
  202. -------- bounded versions, generating primes UP TO a LIMIT: -------------------------------------------
  203.  
  204. primesLimT  :: Int -> [Int]            
  205. primesLimT  m = sieve [2..m] where
  206.   sieve (p:xs) = p : sieve [x | x <- xs, rem x p /= 0]           -- (* Turner's Trial Division/lim *)
  207.   sieve []     = []
  208.  
  209. primesLimGT :: Int -> [Int]            
  210. primesLimGT m = sieve [2..m] where
  211.   sieve (p:xs) 
  212.    | p*p > m   = p : xs
  213.    | True      = p : sieve [x | x <- xs, rem x p /= 0]           -- (*           ... guarded       *) 
  214.   sieve []     = []
  215.  
  216. primesLimE  :: Int -> [Int] 
  217. primesLimE  m = sieve [2..m] where
  218.   sieve (p:xs) = p : sieve (minus xs [p,p+p..m])                  -- (* Sieve of Eratosthenes/lim   *)  
  219.   sieve []     = []
  220.  
  221. primesLimGE :: Int -> [Int] 
  222. primesLimGE m = sieve [2..m] where                                -- (*           ... guarded       *) 
  223.   sieve (p:xs)
  224.    | p*p > m   = p : xs                                -- to STOP EARLY is essential!
  225.    | True      = p : sieve (minus xs [p,p+p..m]) 
  226.   sieve []     = []
  227.  
  228. primesLimGQE :: Int -> [Int]              
  229. primesLimGQE m = 2 : sieve [3,5..m] where                         -- (*            ... from squares *)  
  230.   sieve (p:xs) 
  231.    | p*p > m   = p : xs                  -- early cut-off on first p such that p*p>m
  232.    | True      = p : sieve (minus xs [p*p,p*p+2*p..m]) -- to start from squares is just small improvement
  233.   sieve []     = []
  234.  
  235. primesLimQE :: Int -> [Int]              -- no Guarded:
  236. primesLimQE m = 2 : sieve [3,5..m] where                          -- (*            ... from squares *)  
  237.   sieve (p:xs) 
  238.    | p > mroot = p : sieve (minus xs []) -- no early cut-off, just INT wraparound guarded off
  239.    | True      = p : sieve (minus xs [p*p,p*p+2*p..m]) -- prevent out-of-bounds (p*p)
  240.   sieve []     = []  
  241.   mroot = floor $ sqrt $ fromIntegral m + 1
  242.  
  243. primesLimGEU :: Int -> [Int] 
  244. primesLimGEU m = 2 : sieve [3,5..m] where                         -- (* Euler's Sieve/lim, guarded  *) 
  245.   sieve (p:xs) 
  246.    | p*p > m   = p : xs                                -- ... prevent INT wraparound too
  247.    | True      = p : sieve (minus xs [p*x | x <- p:xs]) 
  248.   sieve []     = []
  249.  
  250. {-========================================================================================
  251.           2,000           10,000           15,000             17,000             19,000
  252. ------------------------------------------------------------------------------------------
  253. T     0.08s-3.7MB      2.65s-4.7MB       7.01s-5.8MB        9.52s-5.8MB       12.83s-5.8MB     Turner's
  254.                  n^2.17            n^2.40            n^2.45            n^2.68
  255. LimT                    2.72-4.7         7.06-5.8                                                   /limit 
  256. ------------------------------------------------------------------------------------------
  257. E     0.05s-4.8MB      2.10s-5.8MB       5.95s-6.8MB        8.82s-6.8MB       14.00s-7.8MB     Eratosthenes'
  258.                  n^2.30            n^2.60            n^3.14            n^4.15 
  259. LimE                    1.17-4.7          3.15-4.7         25k:11.55-5.8                            /limit
  260.                                    n^2.44            n^2.54
  261. 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
  262.                  n^1.97            n^1.28            n^1.12            n^1.03
  263. ------------------------------------------------------------------------------------------
  264. EU    0.41s-48.8MB   4k:1.71s-167.6MB                                                          EUler's
  265.                  n^2.06 
  266. LimGEU               4k:0.04- 8.8       15k:0.27-24.2      50k:2.06-166.5                          /limit/guard
  267.                                    n^1.44            n^1.69
  268. ==========================================================================================
  269.  
  270. ==========================================================================================
  271.          78,499           200,000          400,000          1,000,000          2,000,000
  272. ------------------------------------------------------------------------------------------
  273. LimGQAE 3.59-20.2    150k:10.34-36.5                                                           E/arrays
  274.                  n^1.63
  275. LimGT   0.65-3.7       2.56s-3.7MB       6.93s- 3.7MB     650k:14.07-3.7                       T/guard /lim
  276.                  n^1.47            n^1.44            n^1.46
  277. PT    0.48s-5.8MB      1.85s-7.8MB       5.08s-11.9MB     800k:13.90-20.1                      Postponed T
  278.                  n^1.44            n^1.46            n^1.45 
  279. ST    0.36s-5.8MB      1.38s-7.9MB       3.72s-12.0MB      14.04s-37.6MB                       Segmented T
  280.                  n^1.44            n^1.43            n^1.45
  281. ------------------------------------------------------------------------------------------
  282. LimGE   0.46-4.8       1.62s-4.8MB       4.35s- 4.8MB     900k:14.05-4.8                       E/guard /lim
  283.                  n^1.35            n^1.43            n^1.45
  284. LimQE   0.39-4.8       1.46s-4.8MB       3.88s- 4.8MB      14.81s- 4.8MB                       E/sQ    /lim
  285.                  n^1.41            n^1.41            n^1.46
  286. LimGQE  0.38-4.8       1.35s-4.8MB       3.59s- 4.8MB      13.79s- 4.8MB                       E/sQ /g /lim
  287.                  n^1.36            n^1.41            n^1.47
  288. PE    0.29s-6.8MB      1.12s-11.9MB      3.01s-22.1MB      11.60s-57.0MB                       Postponed E
  289.                  n^1.44            n^1.43            n^1.47 
  290. ------------------------------------------------------------------------------------------
  291. LME   0.26s-4.8MB      0.95s- 4.8MB      2.52s- 4.8MB       9.33s- 4.8MB                       Linear Merge
  292.                  n^1.39            n^1.41            n^1.43                                          2x feed
  293. ------------------------------------------------------------------------------------------
  294. TME   0.16s-4.8MB      0.49s- 4.8MB      1.10s- 4.8MB       3.35s- 4.8MB     7.70s- 4.8MB      Tree Merge
  295.                  n^1.20            n^1.17            n^1.22             n^1.20 
  296. TMWE  0.07s-4.8MB      0.19s- 4.8MB      0.44s- 4.8MB       1.33s- 4.8MB     3.14s- 4.8MB      TME w/Wheel
  297.                  n^1.07            n^1.21            n^1.21             n^1.24 
  298. PQMWE 0.04s-3.8MB      0.14s- 3.8MB      0.32s- 3.8MB       0.98s- 4.8MB     2.25s- 4.8MB   Join w/O'Neill's PQ
  299.                  n^1.34            n^1.19            n^1.22             n^1.20 
  300.       ====================================================================================
  301.           2m                4m                 6m                7m              8m
  302. TME   7.70 - 4.8      3.4m:14.78-4.8
  303.                  n^1.23 
  304. TMWE  3.14s- 4.8MB      7.44s- 4.8MB      12.23s- 4.8MB     14.81s-4.8MB
  305.                  n^1.24            n^1.23             n^1.24 
  306. H8Y5V 2.74 - 4.8        6.23 - 4.8        10.14 - 4.8       12.06 - 4.8      14.20 - 4.8       ONeill original
  307.                  n^1.19            n^1.20             n^1.12            n^1.22           ^1.19
  308. PQMWE 2.25s- 4.8MB      5.19s- 4.8MB       8.48s- 4.8MB     10.18s-4.8MB     11.99s-4.8MB      PQW w/ONeill-PQ
  309.                  n^1.21            n^1.21             n^1.19            n^1.23           ^1.21
  310. lKxjW 2.21 - 4.9        5.09 - 4.9         8.30 - 4.9       10.02 -4.9       11.79s-4.9MB      imprvd ONeill sv
  311.                  n^1.20            n^1.21             n^1.22            n^1.22           ^1.21
  312. ==========================================================================================-}
  313.  
  314. -- *) --- double-feed idea by Melissa O'Neill --- initial tree-folding idea by Dave Bayer 
  315.  
  316.  
  317.  
  318. -- now THIS -- PQ-join instead of Tree-join:
  319.                    --  (... ran and produced results OK on the FIRST TRY!!)
  320.                    -- using MINUS instead of GAPS is BAD: 6.41 vs 5.21s for 4,000,000th prime
  321.                    -- try separating out MULTS... YES! *BETTER!*
  322.                    -- may try to fuse `gaps` with `joinPQ` -- see if it makes it faster .........
  323.  
  324. primesPQMWE :: [Int]                                           
  325. primesPQMWE = 2:3:5:7: gaps 11 wheel (joinPQ mults)
  326.  where                                                             --  sep mults:
  327.    primes' = 11:13: gaps 17 (drop 2 wheel) (joinPQ mults)          -- 4m:  NO CHANGE
  328.    mults   = roll 11 wheel primes'                                 -- 6m: 8.48s-4.8MB vs 8.53s-4.8 
  329.    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:        -- 7m: 10.18s-4.8 vs 10.29s-4.8
  330.            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 
  331.    gaps k ws@(w:t) cs@ ~(c:u) | k==c  = gaps (k+w) t u             -- (*  PQ, 2x feed, w/ Wheel!   *)
  332.                               | True  = k : gaps (k+w) t cs  
  333.    roll k ws@(w:t) ps@ ~(p:u) | k==p  = (p,ws) -- scanl (\c d->c+p*d) (p*p) ws  -- with scanl, sep mults
  334.                                           : roll (k+w) t u                     -- is MUCH WORSE:
  335.                               | True  = roll (k+w) t ps                       --  slower and MEM shoots UP!
  336.  
  337. joinPQ ((p,ws) : sts@((p2,_):_) ) = go (singletonPQ (p*p) (p,ws)) (p2*p2) sts
  338.   where
  339.     go pq q sts =                                -- pq is good up to q, the 1st one in sts
  340.        let  k1 = minKeyPQ pq                     -- next composite to produce
  341.        in  if k1 < q 
  342.              then k1 : go (next k1 pq) q sts
  343.              else 
  344.                   let (a : t@((p2,_):_) ) = sts
  345.                   in  go (insertPQ q a pq) (p2*p2) t
  346.     next k1 pq =                                 -- k1 is the top -- get all equals off 
  347.        let  (c,(p,d:ws)) = minKeyValuePQ pq      -- and reinsert with next values
  348.        in  if c==k1 
  349.              then next k1 $ deleteMinAndInsertPQ (c+p*d) (p,ws) pq
  350.              else pq
  351.  
  352.  
  353.  
  354. -- next step: fuse gaps + joinPQ ==> gapsPQ  .. IN ERROR yet ..
  355.  
  356. primesFPQW :: [Int]                                           
  357. primesFPQW = 2:3:5:7: gapsPQ 11 wheel mults
  358.  where                                                   
  359.    primes' = 11:13: gapsPQ 17 (drop 2 wheel) mults 
  360.    mults   = roll 11 wheel primes'                 
  361.    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:       
  362.            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 -- (*  PQ, 2x feed, w/ Wheel!   *)
  363.    roll k ws@(w:t) ps@ ~(p:u) | k==p  = (p,ws) : roll (k+w) t u 
  364.                               | True  = roll (k+w) t ps    
  365.  
  366. gapsPQ k0 ws0 ((p,ws) : sts@((p2,_):_) ) = go k0 ws0 (singletonPQ (p*p) (p,ws)) (p2*p2) sts
  367.   where
  368.     go k0 ws0@(w0:t0) pq q sts =                 -- pq is good up to q, the 1st one in sts
  369.        let  k1 = minKeyPQ pq                     -- next composite to produce
  370.        in  
  371.         if k0 < k1  
  372.          then -- k0 is the next prime; return it without advancing the PQ
  373.                   k0 : go (k0+w0) t0 pq q sts
  374.          else 
  375.           -- if k0 == k1  -- is composite; ignore and advance as needed
  376.           -- then
  377.             if k1 < q 
  378.              then  --   k1 : go (next k1 pq) q sts
  379.                    -- gaps k0 ws0 (k1 : go (next k1 pq) q sts)
  380.                       go (k0+w0) t0 (next k1 pq) q sts
  381.              else 
  382.                   let (a : t@((p2,_):_) ) = sts
  383.                   in  go (k0+w0) t0 (insertPQ q a pq) (p2*p2) t
  384.           -- else -- k0 > k1
  385.  
  386.     next k1 pq =                                 -- k1 is the top -- get all equals off 
  387.        let  (c,(p,d:ws)) = minKeyValuePQ pq      -- and reinsert with next values
  388.        in  if c==k1 
  389.              then next k1 $ deleteMinAndInsertPQ (c+p*d) (p,ws) pq
  390.              else pq
  391.  
  392.  
  393.  
  394. -- http://e...content-available-to-author-only...s.org/Priority_Queue_(Haskell) is almost TWICE SLOWER than O'Neill's PQ
  395.  
  396. -- O'Neill's PQ:
  397.  
  398. --------------------------------------------------------------------------------
  399. -- ======== the file ONeillPrimes.hs ========
  400.  
  401. -- Generate Primes using ideas from The Sieve of Eratosthenes
  402. --
  403. -- This code is intended to be a faithful reproduction of the
  404. -- Sieve of Eratosthenes, with the following change from the original
  405. --   - The list of primes is infinite
  406. -- (This change does have consequences for representing the number table
  407. -- from which composites are "crossed out".)
  408. --
  409. -- (c) 2006-2007 Melissa O'Neill.  Code may be used freely so long as
  410. -- this copyright message is retained and changed versions of the file
  411. -- are clearly marked.
  412.  
  413. ----- // next line is commented out for Ideone.com testing (willness-2011-02-13) // ---
  414. -- module ONeillPrimes (primes, sieve, calcPrimes, primesToNth, primesToLimit) where
  415.  
  416.  
  417. -- Priority Queues;  this is essentially a copy-and-paste-job of
  418. -- PriorityQ.hs.  By putting the code here, we allow the Haskell
  419. -- compiler to do some whole-program optimization.  (Based on ML
  420. -- code by L. Paulson in _ML for the Working Programmer_.)
  421.  
  422. data PriorityQ k v = Lf -- ***********
  423.                    | Br {-# UNPACK #-} !k v !(PriorityQ k v) !(PriorityQ k v)
  424.                deriving (Eq, Ord, Read, Show)
  425.  
  426. emptyPQ :: PriorityQ k v
  427. emptyPQ = Lf
  428.  
  429. isEmptyPQ :: PriorityQ k v -> Bool
  430. isEmptyPQ Lf  = True
  431. isEmptyPQ _   = False
  432.  
  433. minKeyValuePQ :: PriorityQ k v -> (k, v)
  434. minKeyValuePQ (Br k v _ _)    = (k,v)
  435. minKeyValuePQ _               = error "Empty heap!"
  436.  
  437. minKeyPQ :: PriorityQ k v -> k
  438. minKeyPQ (Br k v _ _)         = k
  439. minKeyPQ _                    = error "Empty heap!"
  440.  
  441. minValuePQ :: PriorityQ k v -> v
  442. minValuePQ (Br k v _ _)       = v
  443. minValuePQ _                  = error "Empty heap!"
  444.  
  445. insertPQ :: Ord k => k -> v -> PriorityQ k v -> PriorityQ k v
  446. insertPQ wk wv (Br vk vv t1 t2)
  447.                | wk <= vk   = Br wk wv (insertPQ vk vv t2) t1
  448.                | otherwise  = Br vk vv (insertPQ wk wv t2) t1
  449. insertPQ wk wv Lf             = Br wk wv Lf Lf
  450.  
  451. siftdown :: Ord k => k -> v -> PriorityQ k v -> PriorityQ k v -> PriorityQ k v
  452. siftdown wk wv Lf _             = Br wk wv Lf Lf
  453. siftdown wk wv (t @ (Br vk vv _ _)) Lf
  454.     | wk <= vk                  = Br wk wv t Lf
  455.     | otherwise                 = Br vk vv (Br wk wv Lf Lf) Lf
  456. siftdown wk wv (t1 @ (Br vk1 vv1 p1 q1)) (t2 @ (Br vk2 vv2 p2 q2))
  457.     | wk <= vk1 && wk <= vk2    = Br wk wv t1 t2
  458.     | vk1 <= vk2                = Br vk1 vv1 (siftdown wk wv p1 q1) t2
  459.     | otherwise                 = Br vk2 vv2 t1 (siftdown wk wv p2 q2)
  460.  
  461. deleteMinAndInsertPQ :: Ord k => k -> v -> PriorityQ k v -> PriorityQ k v
  462. deleteMinAndInsertPQ wk wv Lf             = error "Empty PriorityQ"
  463. deleteMinAndInsertPQ wk wv (Br _ _ t1 t2) = siftdown wk wv t1 t2
  464.  
  465. leftrem :: PriorityQ k v -> (k, v, PriorityQ k v)
  466. leftrem (Br vk vv Lf Lf) = (vk, vv, Lf)
  467. leftrem (Br vk vv t1 t2) = (wk, wv, Br vk vv t t2) where
  468.     (wk, wv, t) = leftrem t1
  469. leftrem _                = error "Empty heap!"
  470.  
  471. deleteMinPQ :: Ord k => PriorityQ k v -> PriorityQ k v
  472. deleteMinPQ (Br vk vv Lf _) = Lf
  473. deleteMinPQ (Br vk vv t1 t2) = siftdown wk wv t2 t where
  474.     (wk,wv,t) = leftrem t1
  475. deleteMinPQ _ = error "Empty heap!"
  476.  
  477. ----- // the rest of the file is omitted here (willness-2011-05-20) // ---
  478.  
  479. ----- an interface function (willness-2011-05-20) --
  480.  
  481. singletonPQ k v    = insertPQ k v emptyPQ
stdin
copy 
2000000 

 78499  1000003   150000 2015177  200000 2750159   400000 5800079  1000000 15485863   900000 13834103     

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 
 150000 2015177
 100000 1299709 
 78499  1000003 
 50000  611953 
 25000  287117 
 15000  163841 
 10000  104729 
 4000   37813 

2m 32452843
4m 67867967
5m 86028121
6m 104395301
7m 122949823
8m 141650939
9m 160481183   (13.56s lKxjW)
compilation info 
[1 of 1] Compiling Main             ( prog.hs, prog.o )
Linking prog ...
stdout
copy 
[32452657,32452681,32452687,32452727,32452759,32452781,32452789,32452837,32452841,32452843]
https://ideone.com/GMkfS
language:
Haskell (ghc 8.4.4)
created:
12 years ago
visibility:
secret

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!