fork download
  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
1000000
compilation info
[1 of 1] Compiling Main             ( prog.hs, prog.o )
Linking prog ...
stdout
(1000000,15485863)