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
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
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. -}
Success #stdin compilation error #stdout 1.21s 4816KB
stdin
1000000
compilation info
[1 of 1] Compiling Main             ( prog.hs, prog.o )