fork(1) download
copy 
  1. {-# OPTIONS_GHC -O2 -fno-cse #-}
  2.  
  3. -- // the above specifies options common with our testing (willness-2011-02-13) // ---
  4.  
  5. --------------------------------------------------------------------------------
  6. -- ======== the file ONeillPrimesTest.hs ========
  7.  
  8. -- Print out the nth prime, where n is the 1st argument
  9.  
  10. module Main where
  11.  
  12. ----- // next two lines commented out for Ideone.com testing (willness-2011-02-13) // ---
  13. -- import ONeillPrimes (primes)
  14. -- import System (getArgs)
  15.  
  16. --- The largest prime we can safely calculate with Int arithmetic is the
  17. --- 105,095,435th prime, which is 2,147,437,291 (prime closest to
  18. --- 2**31 - sqrt(2**31), but our nThPrimeApprox function overestimates,
  19. --- so we have to limit choosing an the Int specializaton over the Integer
  20. --- one at 103,945,550.
  21. printNthPrime :: Int -> IO ()
  22. printNthPrime n
  23.     | n < 103945550 = print (n, (primes !! (n - 1)) :: Int)
  24.     | otherwise     = print (n, (primes !! (n - 1)) :: Integer)
  25.  
  26. main = do
  27.     ---- // next two lines are replaced for Ideone.com testing (willness-2011-02-13) // ---
  28.     -- args <- getArgs
  29.     -- printNthPrime $ read $ head args
  30.     arg <- getLine
  31.     printNthPrime $ read $ arg
  32.  
  33.  
  34. --------------------------------------------------------------------------------
  35. -- ======== the file ONeillPrimes.hs ========
  36.  
  37. -- Generate Primes using ideas from The Sieve of Eratosthenes
  38. --
  39. -- This code is intended to be a faithful reproduction of the
  40. -- Sieve of Eratosthenes, with the following change from the original
  41. --   - The list of primes is infinite
  42. -- (This change does have consequences for representing the number table
  43. -- from which composites are "crossed out".)
  44. --
  45. -- (c) 2006-2007 Melissa O'Neill.  Code may be used freely so long as
  46. -- this copyright message is retained and changed versions of the file
  47. -- are clearly marked.
  48.  
  49. ----- // next line is commented out for Ideone.com testing (willness-2011-02-13) // ---
  50. -- module ONeillPrimes (primes, sieve, calcPrimes, primesToNth, primesToLimit) where
  51.  
  52.  
  53. -- Priority Queues;  this is essentially a copy-and-paste-job of
  54. -- PriorityQ.hs.  By putting the code here, we allow the Haskell
  55. -- compiler to do some whole-program optimization.  (Based on ML
  56. -- code by L. Paulson in _ML for the Working Programmer_.)
  57.  
  58. data PriorityQ k v = Lf
  59.                    | Br {-# UNPACK #-} !k v !(PriorityQ k v) !(PriorityQ k v)
  60.                deriving (Eq, Ord, Read, Show)
  61.  
  62. emptyPQ :: PriorityQ k v
  63. emptyPQ = Lf
  64.  
  65. isEmptyPQ :: PriorityQ k v -> Bool
  66. isEmptyPQ Lf  = True
  67. isEmptyPQ _   = False
  68.  
  69. minKeyValuePQ :: PriorityQ k v -> (k, v)
  70. minKeyValuePQ (Br k v _ _)    = (k,v)
  71. minKeyValuePQ _               = error "Empty heap!"
  72.  
  73. minKeyPQ :: PriorityQ k v -> k
  74. minKeyPQ (Br k v _ _)         = k
  75. minKeyPQ _                    = error "Empty heap!"
  76.  
  77. minValuePQ :: PriorityQ k v -> v
  78. minValuePQ (Br k v _ _)       = v
  79. minValuePQ _                  = error "Empty heap!"
  80.  
  81. insertPQ :: Ord k => k -> v -> PriorityQ k v -> PriorityQ k v
  82. insertPQ wk wv (Br vk vv t1 t2)
  83.                | wk <= vk   = Br wk wv (insertPQ vk vv t2) t1
  84.                | otherwise  = Br vk vv (insertPQ wk wv t2) t1
  85. insertPQ wk wv Lf             = Br wk wv Lf Lf
  86.  
  87. siftdown :: Ord k => k -> v -> PriorityQ k v -> PriorityQ k v -> PriorityQ k v
  88. siftdown wk wv Lf _             = Br wk wv Lf Lf
  89. siftdown wk wv (t @ (Br vk vv _ _)) Lf
  90.     | wk <= vk                  = Br wk wv t Lf
  91.     | otherwise                 = Br vk vv (Br wk wv Lf Lf) Lf
  92. siftdown wk wv (t1 @ (Br vk1 vv1 p1 q1)) (t2 @ (Br vk2 vv2 p2 q2))
  93.     | wk <= vk1 && wk <= vk2    = Br wk wv t1 t2
  94.     | vk1 <= vk2                = Br vk1 vv1 (siftdown wk wv p1 q1) t2
  95.     | otherwise                 = Br vk2 vv2 t1 (siftdown wk wv p2 q2)
  96.  
  97. deleteMinAndInsertPQ :: Ord k => k -> v -> PriorityQ k v -> PriorityQ k v
  98. deleteMinAndInsertPQ wk wv Lf             = error "Empty PriorityQ"
  99. deleteMinAndInsertPQ wk wv (Br _ _ t1 t2) = siftdown wk wv t1 t2
  100.  
  101. leftrem :: PriorityQ k v -> (k, v, PriorityQ k v)
  102. leftrem (Br vk vv Lf Lf) = (vk, vv, Lf)
  103. leftrem (Br vk vv t1 t2) = (wk, wv, Br vk vv t t2) where
  104.     (wk, wv, t) = leftrem t1
  105. leftrem _                = error "Empty heap!"
  106.  
  107. deleteMinPQ :: Ord k => PriorityQ k v -> PriorityQ k v
  108. deleteMinPQ (Br vk vv Lf _) = Lf
  109. deleteMinPQ (Br vk vv t1 t2) = siftdown wk wv t2 t where
  110.     (wk,wv,t) = leftrem t1
  111. deleteMinPQ _ = error "Empty heap!"
  112.  
  113. -- A hybrid of Priority Queues and regular queues.  It allows a priority
  114. -- queue to have a feeder queue, filled with items that come in an
  115. -- increasing order.  By keeping the feed for the queue separate, we
  116. -- avoid needlessly filling an O(log n) data structure with data that
  117. -- it won't need for a long time.
  118.  
  119. type HybridQ k v = (PriorityQ k v, [(k,v)])
  120.  
  121. initHQ :: PriorityQ k v -> [(k,v)] -> HybridQ k v
  122. initHQ pq feeder = (pq, feeder)
  123.  
  124. insertHQ :: (Ord k) => k -> v -> HybridQ k v -> HybridQ k v
  125. insertHQ k v (pq, q) = (insertPQ k v pq, q)
  126.  
  127. deleteMinAndInsertHQ :: (Ord k) => k -> v -> HybridQ k v -> HybridQ k v
  128. deleteMinAndInsertHQ k v (pq, q) = postRemoveHQ(deleteMinAndInsertPQ k v pq, q)
  129.     where
  130.         postRemoveHQ mq@(pq, []) = mq
  131.         postRemoveHQ mq@(pq, (qk,qv) : qs)
  132.             | qk < minKeyPQ pq = (insertPQ qk qv pq, qs)
  133.             | otherwise        = mq
  134.  
  135. minKeyHQ      :: HybridQ k v -> k
  136. minKeyHQ (pq, q) = minKeyPQ pq
  137.  
  138. minKeyValueHQ :: HybridQ k v -> (k, v)
  139. minKeyValueHQ (pq, q) = minKeyValuePQ pq
  140.  
  141.  
  142. -- Finally, we have acceptable queues, now on to finding ourselves some
  143. -- primes.
  144.  
  145.  
  146. -- Here we use a wheel to generate all the number that are not multiples
  147. -- of 2, 3, 5, and 7.  We use some hard-coded data for that.
  148.  
  149. {-# SPECIALIZE wheel :: [Int] #-}
  150. {-# SPECIALIZE wheel :: [Integer] #-}
  151. wheel :: Integral a => [a]
  152. 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
  153. 	:2:6:6:4:2:4:6:2:6:4:2:4:2:10:2:10:wheel
  154.  
  155. -- Now generate the primes using that wheel
  156.  
  157. -- Sometimes, memoization isn't your friend.  Maybe you don't actually want
  158. -- to remember all the primes for the duration of your program and doing so
  159. -- is just wasted space.  For that situation, we provide calcPrimes which
  160. -- calculates the infinite list of primes from scratch.
  161.  
  162. {-# SPECIALIZE calcPrimes :: () -> [Int] #-}
  163. {-# SPECIALIZE calcPrimes :: () -> [Integer] #-}
  164. calcPrimes :: Integral a => () -> [a]
  165. calcPrimes () = 2 : 3 : 5 : 7 : sieve 11 wheel
  166.  
  167. {-# SPECIALIZE primes :: [Int] #-}
  168. {-# SPECIALIZE primes :: [Integer] #-}
  169. primes :: Integral a => [a]
  170. primes = calcPrimes ()
  171.  
  172. {-# SPECIALIZE primesToNth :: Int -> [Integer] #-}
  173. {-# SPECIALIZE primesToNth :: Int -> [Int] #-}
  174. primesToNth :: Integral a => Int -> [a]
  175. primesToNth n = take n (calcPrimes ())
  176.  
  177. {-# SPECIALIZE primesToLimit :: Integer -> [Integer] #-}
  178. {-# SPECIALIZE primesToLimit :: Int -> [Int] #-}
  179. primesToLimit :: Integral a => a -> [a]
  180. primesToLimit limit = takeWhile (< limit) (calcPrimes ())
  181.  
  182. -- This version of the sieve takes a wheel, not a list to be sieved.
  183. -- primes1 and primes2 represent the same infinite list of, but they are
  184. -- consumed at different speeds.  By creating separate expressions, we
  185. -- avoid retaining all the material between the two points.  Sometimes
  186. -- (when you care about space usage) memoization is not your friend.
  187.  
  188. {-# SPECIALIZE sieve :: Int -> [Int] -> [Int] #-}
  189. {-# SPECIALIZE sieve :: Integer -> [Integer] -> [Integer] #-}
  190. sieve :: Integral a => a -> [a] -> [a]
  191. sieve n [] = []
  192. sieve n wheel@(d:ds) = n : (map (\(p,wheel) -> p) primes1) where
  193.     primes1 = sieve' (n+d) ds initialTable
  194.     primes2 = sieve' (n+d) ds initialTable
  195.     initialTable = initHQ (insertPQ (n*n) (n, wheel) emptyPQ)
  196.                    (map (\(p,wheel) -> (p*p,(p,wheel))) primes2)
  197.     sieve' n []     table = []
  198.     sieve' n wheel@(d:ds) table
  199.         | nextComposite <= n = sieve' (n+d) ds (adjust table)
  200.         | otherwise	     = (n,wheel) : sieve' (n+d) ds table
  201.         where
  202.             nextComposite = minKeyHQ table
  203.             adjust table
  204.                 | m <= n    = adjust (deleteMinAndInsertHQ m' (p, ds) table)
  205.                 | otherwise = table
  206.               where
  207. 		(m, (p, d:ds)) = minKeyValueHQ table
  208. 		m' = m + p * d
  209. --------------------------------------------------------------------------------
  210. --- // end of ONeill code (willness-2011-02-13) // ---
  211. {-                                           Ojq01         R3TBr
  212.                        ONeill              gaps/roll     vips/roll
  213. (1000000,15485863)   1.19s 4.8MB   x1.16    1.38-4.7      1.46-4.8
  214. (2000000,32452843)   2.74s 4.8MB   x1.19    3.25-4.7      3.36-4.8
  215. (3000000,49979687)   4.41s 4.8MB   x1.22    5.40-4.7      5.55-4.8
  216. (4000000,67867967)   6.23s 4.8MB   x1.23    7.66-4.7      7.92-4.8
  217. (5000000,86028121)   8.07s 4.8MB   x1.26   10.16-4.7     10.40-4.8
  218. (6000000,104395301) 10.14s 4.8MB   x1.25   12.71-4.7     12.99-4.8
  219. (7000000,122949823) 12.06s 4.8MB 
  220. (8000000,141650939) 14.20s 4.8MB 
  221.  
  222. empirical cpxty:     n^1.20                 n^1.24        n^1.22
  223. -}
stdin
copy 
1000000
compilation info 
[1 of 1] Compiling Main             ( prog.hs, prog.o )
Linking prog ...
stdout
copy 
(1000000,15485863)
https://ideone.com/H8Y5V
testing the ONeillPrimes*.hs files by Melissa O'Neill from her http://w...content-available-to-author-only...c.edu/~oneill/code/haskell-primes.zip
language:
Haskell (ghc 8.4.4)
created:
13 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!