module Arithmetic.Prime.Sieve
where
import OpenTheory.Primitive.Natural
import qualified Arithmetic.Utility.Heap as Heap
newtype Sieve = Sieve { unSieve :: Heap.Heap (Natural,Natural) }
instance Show Sieve where
show s = show (unSieve s)
initial :: Sieve
initial =
Sieve (Heap.empty lep)
where
lep (kp1,_) (kp2,_) = kp1 <= kp2
add :: Natural -> Sieve -> (Natural,Sieve)
add m (Sieve ps) =
(p, Sieve (Heap.add (m',p) ps))
where
p = 2 * m + 1
m' = (p + 1) * m
bump :: Sieve -> (Natural,Sieve)
bump (Sieve ps) =
case Heap.remove ps of
Nothing -> error "GenuineSieve.bump"
Just ((kp,p),ps') -> (kp, Sieve (Heap.add (kp + p, p) ps'))
advance :: Natural -> Natural -> Sieve -> [Natural]
advance m n s =
if m < n
then let (p,s') = add m s in p : advance m' n s'
else let (n',s') = bump s in advance (if m == n then m' else m) n' s'
where
m' = m + 1