------------------------------------------------------------------------------- |-- Module : Math.Sieve.ONeill-- Copyright : (c) 2006-2010 Melissa O'Neill-- License : BSD3---- Maintainer : leon@melding-monads.com-- Stability : experimental-- Portability : portable---- Incremental primality sieve based on priority queues, described-- in the paper /The Genuine Sieve of Eratosthenes/ by Melissa O'Neill,-- /Journal of Functional Programming/, 19(1), pp95-106, Jan 2009.---- <http://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf>---- Code is unchanged, other than packaging, from that written by-- Melissa O'Neill. This version contains optimizations not described-- in the paper, primarily improving memory consumption.---- <http://www.cs.hmc.edu/~oneill/code/haskell-primes.zip>--------------------------------------------------------------------------------- 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.moduleMath.Sieve.ONeill(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_.)dataPriorityQkv=Lf|Br{-# UNPACK #-}!kv!(PriorityQkv)!(PriorityQkv)deriving(Eq,Ord,Read,Show)emptyPQ::PriorityQkvemptyPQ=LfisEmptyPQ::PriorityQkv->BoolisEmptyPQLf=TrueisEmptyPQ_=FalseminKeyValuePQ::PriorityQkv->(k,v)minKeyValuePQ(Brkv__)=(k,v)minKeyValuePQ_=error"Empty heap!"minKeyPQ::PriorityQkv->kminKeyPQ(Brkv__)=kminKeyPQ_=error"Empty heap!"minValuePQ::PriorityQkv->vminValuePQ(Brkv__)=vminValuePQ_=error"Empty heap!"insertPQ::Ordk=>k->v->PriorityQkv->PriorityQkvinsertPQwkwv(Brvkvvt1t2)|wk<=vk=Brwkwv(insertPQvkvvt2)t1|otherwise=Brvkvv(insertPQwkwvt2)t1insertPQwkwvLf=BrwkwvLfLfsiftdown::Ordk=>k->v->PriorityQkv->PriorityQkv->PriorityQkvsiftdownwkwvLf_=BrwkwvLfLfsiftdownwkwv(t@(Brvkvv__))Lf|wk<=vk=BrwkwvtLf|otherwise=Brvkvv(BrwkwvLfLf)Lfsiftdownwkwv(t1@(Brvk1vv1p1q1))(t2@(Brvk2vv2p2q2))|wk<=vk1&&wk<=vk2=Brwkwvt1t2|vk1<=vk2=Brvk1vv1(siftdownwkwvp1q1)t2|otherwise=Brvk2vv2t1(siftdownwkwvp2q2)deleteMinAndInsertPQ::Ordk=>k->v->PriorityQkv->PriorityQkvdeleteMinAndInsertPQwkwvLf=error"Empty PriorityQ"deleteMinAndInsertPQwkwv(Br__t1t2)=siftdownwkwvt1t2leftrem::PriorityQkv->(k,v,PriorityQkv)leftrem(BrvkvvLfLf)=(vk,vv,Lf)leftrem(Brvkvvt1t2)=(wk,wv,Brvkvvtt2)where(wk,wv,t)=leftremt1leftrem_=error"Empty heap!"deleteMinPQ::Ordk=>PriorityQkv->PriorityQkvdeleteMinPQ(BrvkvvLf_)=LfdeleteMinPQ(Brvkvvt1t2)=siftdownwkwvt2twhere(wk,wv,t)=leftremt1deleteMinPQ_=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.typeHybridQkv=(PriorityQkv,[(k,v)])initHQ::PriorityQkv->[(k,v)]->HybridQkvinitHQpqfeeder=(pq,feeder)insertHQ::(Ordk)=>k->v->HybridQkv->HybridQkvinsertHQkv(pq,q)=(insertPQkvpq,q)deleteMinAndInsertHQ::(Ordk)=>k->v->HybridQkv->HybridQkvdeleteMinAndInsertHQkv(pq,q)=postRemoveHQ(deleteMinAndInsertPQkvpq,q)wherepostRemoveHQmq@(pq,[])=mqpostRemoveHQmq@(pq,(qk,qv):qs)|qk<minKeyPQpq=(insertPQqkqvpq,qs)|otherwise=mqminKeyHQ::HybridQkv->kminKeyHQ(pq,q)=minKeyPQpqminKeyValueHQ::HybridQkv->(k,v)minKeyValueHQ(pq,q)=minKeyValuePQpq-- 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::Integrala=>[a]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] #-}-- | An infinite stream of primescalcPrimes::Integrala=>()->[a]calcPrimes()=2:3:5:7:sieve11wheel{-# SPECIALIZE primes :: [Int] #-}{-# SPECIALIZE primes :: [Integer] #-}-- | An infinite stream of primesprimes::Integrala=>[a]primes=calcPrimes(){-# SPECIALIZE primesToNth :: Int -> [Integer] #-}{-# SPECIALIZE primesToNth :: Int -> [Int] #-}-- | Returns the first @n@ primesprimesToNth::Integrala=>Int->[a]primesToNthn=taken(calcPrimes()){-# SPECIALIZE primesToLimit :: Integer -> [Integer] #-}{-# SPECIALIZE primesToLimit :: Int -> [Int] #-}-- | Returns primes up to some limitprimesToLimit::Integrala=>a->[a]primesToLimitlimit=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] #-}-- | The first argument specifies which number to start at,-- the second argument is a wheel of deltas for skipping composites.-- For example, @primes@ could be defined as---- > 2 : 3 : sieve 5 (cycle [2,4])sieve::Integrala=>a->[a]->[a]sieven[]=[]sievenwheel@(d:ds)=n:(map(\(p,wheel)->p)primes1)whereprimes1=sieve'(n+d)dsinitialTableprimes2=sieve'(n+d)dsinitialTableinitialTable=initHQ(insertPQ(n*n)(n,wheel)emptyPQ)(map(\(p,wheel)->(p*p,(p,wheel)))primes2)sieve'n[]table=[]sieve'nwheel@(d:ds)table|nextComposite<=n=sieve'(n+d)ds(adjusttable)|otherwise=(n,wheel):sieve'(n+d)dstablewherenextComposite=minKeyHQtableadjusttable|m<=n=adjust(deleteMinAndInsertHQm'(p,ds)table)|otherwise=tablewhere(m,(p,d:ds))=minKeyValueHQtablem'=m+p*d