{-# OPTIONS_GHC -O2 -fno-cse #-}
-- // the above specifies options common with our testing (willness-2011-02-13) // ---
--------------------------------------------------------------------------------
-- ======== the file ONeillPrimesTest.hs ========
-- Print out the nth prime, where n is the 1st argument
module Main where
----- // next two lines commented out for Ideone.com testing (willness-2011-02-13) // ---
-- import ONeillPrimes (primes)
-- import System (getArgs)
--- The largest prime we can safely calculate with Int arithmetic is the
--- 105,095,435th prime, which is 2,147,437,291 (prime closest to
--- 2**31 - sqrt(2**31), but our nThPrimeApprox function overestimates,
--- so we have to limit choosing an the Int specializaton over the Integer
--- one at 103,945,550.
printNthPrime
:: Int -> IO ()printNthPrime n
| n
< 103945550 = print (n
, (primes
!! (n
- 1)) :: Int)
main = do
---- // next two lines are replaced for Ideone.com testing (willness-2011-02-13) // ---
-- args <- getArgs
-- printNthPrime $ read $ head args
printNthPrime
$ read $ arg
--------------------------------------------------------------------------------
-- ======== the file ONeillPrimes.hs ========
-- Generate Primes using ideas from The Sieve of Eratosthenes
--
-- This code is intended to be a faithful reproduction of the
-- Sieve of Eratosthenes, with the following change from the original
-- - The list of primes is infinite
-- (This change does have consequences for representing the number table
-- from which composites are "crossed out".)
--
-- (c) 2006-2007 Melissa O'Neill. Code may be used freely so long as
-- this copyright message is retained and changed versions of the file
-- are clearly marked.
----- // next line is commented out for Ideone.com testing (willness-2011-02-13) // ---
-- module ONeillPrimes (primes, sieve, calcPrimes, primesToNth, primesToLimit) where
-- Priority Queues; this is essentially a copy-and-paste-job of
-- PriorityQ.hs. By putting the code here, we allow the Haskell
-- compiler to do some whole-program optimization. (Based on ML
-- code by L. Paulson in _ML for the Working Programmer_.)
data PriorityQ k v = Lf
| Br {-# UNPACK #-} !k v !(PriorityQ k v) !(PriorityQ k v)
emptyPQ :: PriorityQ k v
emptyPQ = Lf
isEmptyPQ
:: PriorityQ k v
-> BoolisEmptyPQ Lf = True
isEmptyPQ _ = False
minKeyValuePQ :: PriorityQ k v -> (k, v)
minKeyValuePQ (Br k v _ _) = (k,v)
minKeyValuePQ
_ = error "Empty heap!"
minKeyPQ :: PriorityQ k v -> k
minKeyPQ (Br k v _ _) = k
minKeyPQ
_ = error "Empty heap!"
minValuePQ :: PriorityQ k v -> v
minValuePQ (Br k v _ _) = v
minValuePQ
_ = error "Empty heap!"
insertPQ
:: Ord k
=> k
-> v
-> PriorityQ k v
-> PriorityQ k v
insertPQ wk wv (Br vk vv t1 t2)
| wk <= vk = Br wk wv (insertPQ vk vv t2) t1
| otherwise = Br vk vv
(insertPQ wk wv t2
) t1
insertPQ wk wv Lf = Br wk wv Lf Lf
siftdown
:: Ord k
=> k
-> v
-> PriorityQ k v
-> PriorityQ k v
-> PriorityQ k v
siftdown wk wv Lf _ = Br wk wv Lf Lf
siftdown wk wv (t @ (Br vk vv _ _)) Lf
| wk <= vk = Br wk wv t Lf
siftdown wk wv (t1 @ (Br vk1 vv1 p1 q1)) (t2 @ (Br vk2 vv2 p2 q2))
| wk <= vk1 && wk <= vk2 = Br wk wv t1 t2
| vk1 <= vk2 = Br vk1 vv1 (siftdown wk wv p1 q1) t2
| otherwise = Br vk2 vv2 t1
(siftdown wk wv p2 q2
)
deleteMinAndInsertPQ
:: Ord k
=> k
-> v
-> PriorityQ k v
-> PriorityQ k v
deleteMinAndInsertPQ wk wv Lf
= error "Empty PriorityQ"deleteMinAndInsertPQ wk wv (Br _ _ t1 t2) = siftdown wk wv t1 t2
leftrem :: PriorityQ k v -> (k, v, PriorityQ k v)
leftrem (Br vk vv Lf Lf) = (vk, vv, Lf)
leftrem (Br vk vv t1 t2) = (wk, wv, Br vk vv t t2) where
(wk, wv, t) = leftrem t1
leftrem
_ = error "Empty heap!"
deleteMinPQ
:: Ord k
=> PriorityQ k v
-> PriorityQ k v
deleteMinPQ (Br vk vv Lf _) = Lf
deleteMinPQ (Br vk vv t1 t2) = siftdown wk wv t2 t where
(wk,wv,t) = leftrem t1
deleteMinPQ
_ = error "Empty heap!"
-- A hybrid of Priority Queues and regular queues. It allows a priority
-- queue to have a feeder queue, filled with items that come in an
-- increasing order. By keeping the feed for the queue separate, we
-- avoid needlessly filling an O(log n) data structure with data that
-- it won't need for a long time.
type HybridQ k v = (PriorityQ k v, [(k,v)])
initHQ :: PriorityQ k v -> [(k,v)] -> HybridQ k v
initHQ pq feeder = (pq, feeder)
insertHQ
:: (Ord k
) => k
-> v
-> HybridQ k v
-> HybridQ k v
insertHQ k v (pq, q) = (insertPQ k v pq, q)
deleteMinAndInsertHQ
:: (Ord k
) => k
-> v
-> HybridQ k v
-> HybridQ k v
deleteMinAndInsertHQ k v (pq, q) = postRemoveHQ(deleteMinAndInsertPQ k v pq, q)
where
postRemoveHQ mq@(pq, []) = mq
postRemoveHQ mq@(pq, (qk,qv) : qs)
| qk < minKeyPQ pq = (insertPQ qk qv pq, qs)
minKeyHQ :: HybridQ k v -> k
minKeyHQ (pq, q) = minKeyPQ pq
minKeyValueHQ :: HybridQ k v -> (k, v)
minKeyValueHQ (pq, q) = minKeyValuePQ pq
-- Finally, we have acceptable queues, now on to finding ourselves some
-- primes.
-- Here we use a wheel to generate all the number that are not multiples
-- of 2, 3, 5, and 7. We use some hard-coded data for that.
{-# SPECIALIZE wheel :: [Int] #-}
{-# SPECIALIZE wheel :: [Integer] #-}
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
:2:6:6:4:2:4:6:2:6:4:2:4:2:10:2:10:wheel
-- Now generate the primes using that wheel
-- Sometimes, memoization isn't your friend. Maybe you don't actually want
-- to remember all the primes for the duration of your program and doing so
-- is just wasted space. For that situation, we provide calcPrimes which
-- calculates the infinite list of primes from scratch.
{-# SPECIALIZE calcPrimes :: () -> [Int] #-}
{-# SPECIALIZE calcPrimes :: () -> [Integer] #-}
calcPrimes () = 2 : 3 : 5 : 7 : sieve 11 wheel
{-# SPECIALIZE primes :: [Int] #-}
{-# SPECIALIZE primes :: [Integer] #-}
primes = calcPrimes ()
{-# SPECIALIZE primesToNth :: Int -> [Integer] #-}
{-# SPECIALIZE primesToNth :: Int -> [Int] #-}
primesToNth n
= take n
(calcPrimes
())
{-# SPECIALIZE primesToLimit :: Integer -> [Integer] #-}
{-# SPECIALIZE primesToLimit :: Int -> [Int] #-}
primesToLimit limit
= takeWhile (< limit
) (calcPrimes
())
-- This version of the sieve takes a wheel, not a list to be sieved.
-- primes1 and primes2 represent the same infinite list of, but they are
-- consumed at different speeds. By creating separate expressions, we
-- avoid retaining all the material between the two points. Sometimes
-- (when you care about space usage) memoization is not your friend.
{-# SPECIALIZE sieve :: Int -> [Int] -> [Int] #-}
{-# SPECIALIZE sieve :: Integer -> [Integer] -> [Integer] #-}
sieve n [] = []
sieve n wheel
@(d:ds
) = n :
(map (\
(p
,wheel
) -> p
) primes1
) where primes1 = sieve' (n+d) ds initialTable
primes2 = sieve' (n+d) ds initialTable
initialTable = initHQ (insertPQ (n*n) (n, wheel) emptyPQ)
(map (\
(p
,wheel
) -> (p
*p
,(p
,wheel
))) primes2
) sieve' n [] table = []
sieve' n wheel@(d:ds) table
| nextComposite <= n = sieve' (n+d) ds (adjust table)
| otherwise = (n,wheel) : sieve' (n+d) ds table
where
nextComposite = minKeyHQ table
adjust table
| m <= n = adjust (deleteMinAndInsertHQ m' (p, ds) table)
| otherwise = table
where
(m, (p, d:ds)) = minKeyValueHQ table
m' = m + p * d
--------------------------------------------------------------------------------
--- // end of ONeill code (willness-2011-02-13) // ---
{- Ojq01 R3TBr
ONeill gaps/roll vips/roll
(1000000,15485863) 1.19s 4.8MB x1.16 1.38-4.7 1.46-4.8
(2000000,32452843) 2.74s 4.8MB x1.19 3.25-4.7 3.36-4.8
(3000000,49979687) 4.41s 4.8MB x1.22 5.40-4.7 5.55-4.8
(4000000,67867967) 6.23s 4.8MB x1.23 7.66-4.7 7.92-4.8
(5000000,86028121) 8.07s 4.8MB x1.26 10.16-4.7 10.40-4.8
(6000000,104395301) 10.14s 4.8MB x1.25 12.71-4.7 12.99-4.8
(7000000,122949823) 12.06s 4.8MB
(8000000,141650939) 14.20s 4.8MB
empirical cpxty: n^1.20 n^1.24 n^1.22
-}