/
Primes.hs
161 lines (142 loc) · 5.67 KB
/
Primes.hs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
-- |
-- Module : Data.Numbers.Primes
-- Copyright : Sebastian Fischer
-- License : BSD3
--
-- Maintainer : Sebastian Fischer (sebf@informatik.uni-kiel.de)
-- Stability : experimental
-- Portability : portable
--
-- This Haskell library provides an efficient lazy wheel sieve for
-- prime generation inspired by /Lazy wheel sieves and spirals of/
-- /primes/ by Colin Runciman
-- (<http://www.cs.york.ac.uk/ftpdir/pub/colin/jfp97lw.ps.gz>) and
-- /The Genuine Sieve of Eratosthenes/ by Melissa O'Neil
-- (<http://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf>).
--
module Data.Numbers.Primes ( primes, wheelSieve ) where
-- |
-- This global constant is an infinite list of prime numbers. It is
-- generated by a lazy wheel sieve and shared across the whole program
-- run. If you are concerned about the memory requirements of sharing
-- many primes you can call the function @wheelSieve@ directly.
--
primes :: [Integer]
primes = wheelSieve 6
-- |
-- This function returns an infinite list of prime numbers by sieving
-- with a wheel that cancels the multiples of the first @n@ primes
-- where @n@ is the argument given to @wheelSieve@. Don't use too
-- large wheels. The number @6@ is a good value to pass to this
-- function. Larger wheels improve the run time at the cost of higher
-- memory requirements.
--
wheelSieve :: Int -- ^ number of primes canceled by the wheel
-> [Integer] -- ^ infinite list of primes
wheelSieve k = reverse ps ++ map head (sieve p (cycle ns))
where (p:ps,ns) = wheel k
-- Auxiliary Definitions
------------------------------------------------------------------------------
-- Sieves prime candidates by computing composites from the result of
-- a recursive call with identical arguments. We could use sharing
-- instead of a recursive call with identical arguments but that would
-- lead to much higher memory requirements. The results of the
-- different calls are consumed at different speeds and we want to
-- avoid multiple far apart pointers into the result list to avoid
-- retaining everything in between.
--
-- Each list in the result starts with a prime. To obtain composites
-- that need to be cancelled, one can multiply all elements of the
-- list with its head.
--
sieve :: Integer -> [Integer] -> [[Integer]]
sieve p ns@(m:ms) = spin p ns : sieveComps (p+m) ms (composites p ns)
-- Composites are stored in increasing order in a priority queue. The
-- queue has an associated feeder which is used to avoid filling it
-- with entries that will only be used again much later.
--
type Composites = (Queue,[[Integer]])
-- The feeder is computed from the result of a call to 'sieve'.
--
composites :: Integer -> [Integer] -> Composites
composites p ns = (Empty, map comps (spin p ns : sieve p ns))
where comps xs@(x:_) = map (x*) xs
-- We can split all composites into the next and remaining
-- composites. We use the feeder when appropriate and discard equal
-- entries to not return a composite twice.
--
splitComposites :: Composites -> (Integer,Composites)
splitComposites (Empty, xs:xss) = splitComposites (Fork xs [], xss)
splitComposites (queue, xss@((x:xs):yss))
| x < z = (x, discard x (enqueue xs queue, yss))
| otherwise = (z, discard z (enqueue zs queue', xss))
where (z:zs,queue') = dequeue queue
-- Drops all occurrences of the given element.
--
discard :: Integer -> Composites -> Composites
discard n ns | n == m = discard n ms
| otherwise = ns
where (m,ms) = splitComposites ns
-- This is the actual sieve. It discards candidates that are
-- composites and yields lists which start with a prime and contain
-- all factors of the composites that need to be dropped.
--
sieveComps :: Integer -> [Integer] -> Composites -> [[Integer]]
sieveComps cand ns@(m:ms) xs
| cand == comp = sieveComps (cand+m) ms ys
| cand < comp = spin cand ns : sieveComps (cand+m) ms xs
| otherwise = sieveComps cand ns ys
where (comp,ys) = splitComposites xs
-- This function computes factors of composites of primes by spinning
-- a wheel.
--
spin :: Integer -> [Integer] -> [Integer]
spin x (y:ys) = x : spin (x+y) ys
-- A wheel consists of a list of primes whose multiples are canceled
-- and the actual wheel that is rolled for canceling.
--
type Wheel = ([Integer],[Integer])
-- Computes a wheel that cancels the multiples of the given number
-- (plus 1) of primes.
--
-- For example:
--
-- wheel 0 = ([2],[1])
-- wheel 1 = ([3,2],[2])
-- wheel 2 = ([5,3,2],[2,4])
-- wheel 3 = ([7,5,3,2],[4,2,4,2,4,6,2,6])
--
wheel :: Int -> Wheel
wheel n = iterate next ([2],[1]) !! n
next :: Wheel -> Wheel
next (ps@(p:_),xs) = (py:ps,cancel (product ps) p py ys)
where (y:ys) = cycle xs
py = p + y
cancel :: Integer -> Integer -> Integer -> [Integer] -> [Integer]
cancel 0 _ _ _ = []
cancel m p n (x:ys@(y:zs))
| nx `mod` p > 0 = x : cancel (m-x) p nx ys
| otherwise = cancel m p n (x+y:zs)
where nx = n + x
-- We use a special version of priority queues implemented as /pairing/
-- /heaps/ (see /Purely Functional Data Structures/ by Chris Okasaki).
--
-- The queue stores non-empty lists of composites; the first element
-- is used as priority.
--
data Queue = Empty | Fork [Integer] [Queue]
enqueue :: [Integer] -> Queue -> Queue
enqueue ns = merge (Fork ns [])
merge :: Queue -> Queue -> Queue
merge Empty y = y
merge x Empty = x
merge x y | prio x <= prio y = join x y
| otherwise = join y x
where prio (Fork (n:_) _) = n
join (Fork ns qs) q = Fork ns (q:qs)
dequeue :: Queue -> ([Integer], Queue)
dequeue (Fork ns qs) = (ns,mergeAll qs)
mergeAll :: [Queue] -> Queue
mergeAll [] = Empty
mergeAll [x] = x
mergeAll (x:y:qs) = merge (merge x y) (mergeAll qs)