Difference between revisions of "Prime numbers"
m 
(mergeSP simplified, renamed unionSP) 

Line 181:  Line 181:  
primes = 2:primes' 
primes = 2:primes' 

where 
where 

−  primes' = [3] ++ [5,7..] `minus` foldr 
+  primes' = [3] ++ [5,7..] `minus` foldr union' [] mults 
mults = map (\p>let q=p*p in (q,tail [q,q+2*p..])) $ primes' 
mults = map (\p>let q=p*p in (q,tail [q,q+2*p..])) $ primes' 

−  +  union' (q,qs) xs = q : union qs xs 

</haskell> 
</haskell> 

This code is yet faster. Its main deficiency still is that it creates a linear nested merging structure, instead of a treelike structure. Each multiple produced by a prime has to percolate to the top eventually, so it's better to make its path shorter. It'll have to go through fewer merge nodes this way. 
This code is yet faster. Its main deficiency still is that it creates a linear nested merging structure, instead of a treelike structure. Each multiple produced by a prime has to percolate to the top eventually, so it's better to make its path shorter. It'll have to go through fewer merge nodes this way. 

−  The linearity is imposed by type asymmetry of our <code>( 
+  The linearity is imposed by type asymmetry of our <code>(union' :: a > b > b)</code> function, forcing us into the <code>1+(1+(1+(1+...)))</code> pattern, <code>+</code> standing for <code>union'</code> (which was defined that way to prevent the runahead when folding over the infinite list of lists of multiples). 
We need to turn it into an associative operation of uniform type <code>(:: a > a > a)</code> to be able to freely rearrange the combinations into arbitrary treelike patterns, as in e.g. ((1+1)+(2+2))+(...) etc. The type uniformity is what makes compositionality possible. 
We need to turn it into an associative operation of uniform type <code>(:: a > a > a)</code> to be able to freely rearrange the combinations into arbitrary treelike patterns, as in e.g. ((1+1)+(2+2))+(...) etc. The type uniformity is what makes compositionality possible. 

Line 201:  Line 201:  
primes' = [3,5,7,11,13,17] ++ [19,21..] `minus` comps 
primes' = [3,5,7,11,13,17] ++ [19,21..] `minus` comps 

mults = map (\p> let q=p*p in ([q],tail [q,q+2*p..])) $ primes' 
mults = map (\p> let q=p*p in ([q],tail [q,q+2*p..])) $ primes' 

−  comps = fst $ tfold 
+  comps = fst $ tfold unionSP (pairwise unionSP mults) 
+  
+  pairwise f (x:y:ys) = f x y : pairwise f ys 

+  tfold f (a:b:c:xs) = (a `f` (b `f` c)) `f` tfold f (pairwise f xs) 

+  unionSP (a:as,bs) v = (a:x,y) where (x,y)=unionSP (as,bs) v 

+  unionSP u@([],b:bs) v@(c:cs,ds) = case compare b c of 

+  LT > (b:x,y) where (x,y)=unionSP ([],bs) v 

+  EQ > (b:x,y) where (x,y)=unionSP ([],bs) (cs,ds) 

+  GT > (c:x,y) where (x,y)=unionSP u (cs,ds) 

+  unionSP ([],bs) ([],ds) = ([],union bs ds) 

−  pairwise f (x:y:ys) = f x y : pairwise f ys 

−  tfold f (a:b:c:xs) = (a `f` (b `f` c)) `f` tfold f (pairwise f xs) 

−  mergeSP (a,b) ~(c,d) = let (bc,bd) = spMerge b c d 

−  in (a ++ bc, bd) 

−  where 

−  spMerge b [] d = ([], union b d) 

−  spMerge b@(x:xs) c@(y:ys) d = case compare x y of 

−  LT > (x:u,v) where (u,v) = spMerge xs c d 

−  EQ > (x:u,v) where (u,v) = spMerge xs ys d 

−  GT > (y:u,v) where (u,v) = spMerge b ys d 

</haskell> 
</haskell> 

The fold used here creates a <code>(2+(2+2))+( (4+(4+4)) + ( (8+(8+8)) + ... ))</code> structure, better adjusted for primes multiples production than <code>1+(2+(4+(8+...)))</code>, used by the "Implicit Heap" code, giving it additional 10% speedup. 
The fold used here creates a <code>(2+(2+2))+( (4+(4+4)) + ( (8+(8+8)) + ... ))</code> structure, better adjusted for primes multiples production than <code>1+(2+(4+(8+...)))</code>, used by the "Implicit Heap" code, giving it additional 10% speedup. 

−  <code> 
+  <code> unionSP </code> is an associative operation, preserving of the invariant such that for a list of multiples <code>[(a,b),(c,d), ...]</code>, it's always the case that <code> last a < head b && last a < head c</code>. These "split pairs" represent ordered list as a pair of its known and (currently) finite prefix, and the rest of it. Such pairs under <code>unionSP</code> operation form a <i>monoid</i>, and if we were to declare a <hask>newtype SplitPair a = SP ([a],[a])</hask> a <hask>Monoid</hask> instance, with <code>unionSP</code> its <hask>mappend</hask> (and <code>tfold</code> its <hask>mconcat</hask>), the above code for <code>comps</code> would just become <hask>SP (comps,_) = mconcat mults</hask>. 
This code exhibits approximately O(<math>{n^{1.20}}</math>)..O(<math>{n^{1.15}}</math>) local asymptotic behavior (tested interpreted, in GHCi, for 10,000 to 300,000 primes produced). When compared with Melissa O'Neill's PQ code from the ZIP package which was modified to work on odds only as well, it is 3.2x times faster, with used memory reported about 2.5x times smaller. 
This code exhibits approximately O(<math>{n^{1.20}}</math>)..O(<math>{n^{1.15}}</math>) local asymptotic behavior (tested interpreted, in GHCi, for 10,000 to 300,000 primes produced). When compared with Melissa O'Neill's PQ code from the ZIP package which was modified to work on odds only as well, it is 3.2x times faster, with used memory reported about 2.5x times smaller. 

Line 229:  Line 219:  
primes' = [11,13,17,19,23,29] ++ (rollFrom 31) `minus` comps 
primes' = [11,13,17,19,23,29] ++ (rollFrom 31) `minus` comps 

mults = map (\p> fromList $ map (p*) $ rollFrom p) $ primes' 
mults = map (\p> fromList $ map (p*) $ rollFrom p) $ primes' 

−  comps = fst $ tfold 
+  comps = fst $ tfold unionSP (pairwise unionSP mults) 
fromList (x:xs) = ([x],xs) 
fromList (x:xs) = ([x],xs) 

rollFrom n = let x = (n11) `mod` 210 
rollFrom n = let x = (n11) `mod` 210 
Revision as of 03:54, 21 February 2010
In mathematics, a prime number (or a prime) is a natural number which has exactly two distinct natural number divisors: 1 and itself. The smallest prime is thus 2.
Contents
 1 Prime Number Resources
 2 Finding Primes
 3 Testing Primality
 4 External links
Prime Number Resources
 At Wikipedia:
 HackageDB packages:
 Numbers: An assortment of number theoretic functions.
 NumberSieves: Number Theoretic Sieves: primes, factorization, and Euler's Totient.
 primes: Efficient, purely functional generation of prime numbers.
Finding Primes
Any natural number is representable as a product of powers of its prime factors, and so a prime has no prime divisors other than itself. That means that starting with 2, for each newly found prime we can eliminate from the rest of the numbers all such numbers that have this prime as their divisor, giving us the next available number as next prime. This is known as sieving the natural numbers, so that in the end what we are left with are just primes.
The Classic Turner's Sieve
Attributed to David Turner (SASL Language Manual, 1975), the following is a direct translation of that idea, generating a list of all prime numbers:
primes :: [Integer]
primes = sieve [2..]
where
sieve (p:xs) = p : sieve [x  x<xs, x `mod` p /= 0]
 or: filter ((/=0).(`mod`p)) xs
This should only be considered a specification, not a code. When run as is, it is extremely inefficient because it starts up the filters prematurely, immediately after each prime, instead of only after the prime's square has been reached. To be admitted as prime, each number will be tested for divisibility here by all its preceding primes, while just those not greater than its square root would suffice. This means that e.g. to find the 1001st prime (7927
), 1000 filters are used, when in fact just the first 24 are needed (up to 89
's filter only).
So this in effect creates a cascade of nested filters in front of the infinite numbers supply, and in extremely premature fashion at that. One way of fixing that would be to postpone the creation of filters until the right moment, by decoupling the primes supply from the numbers supply.
Filtering To Keep the Prime Numbers In
The primes list creation with divisibility testing can be reformulated in a few more ways, using the list of primes as it is being built (a la "circular programming").
Postponed Filters Sieve
Here each filter's creation is postponed until the very moment it's needed.
primes :: [Integer]
primes = 2: 3: sieve (tail primes) [5,7..]
where
sieve (p:ps) xs = h ++ sieve ps [x  x<t, x `rem` p /= 0]
 or: filter ((/=0).(`rem`p)) t
where (h,~(_:t)) = span (< p*p) xs
This can be seen as essential framework for all the code to come. It only tests odd numbers, and only by the primes that are needed, for each numbers span between successive squares of primes. To find the 1001st prime, the divisibility test is performed by only 24 nested filters corresponding to the first 24 odd primes.
Whereas Turner's sieve exhibits near O() behavior, this one exhibits near O() behavior, with an ordersofmagnitude speedup.
Odd numbers, by Trial Division
This is also good for generating a few 100,000s primes (when GHCcompiled as well):
primes :: [Integer]
primes = 2: 3: filter isPrime [5,7..]
where
isPrime n = all (notDivs n)
$ takeWhile (\p> p*p <= n) (tail primes)
notDivs n p = n `mod` p /= 0
Instead of relying on nested filters, it tests each odd number by an explicit list of all the needed prime factors. But for each number tested it refetches this list anew which will be the same for the increasing spans of numbers between the successive squares of primes.
Generated Spans, by Nested Filters
The other way to go instead of concentrating on the numbers supply, is to directly work on the successive spans between the primes squares.
This version is a bit faster still, creating 158,000 primes (again, GHCcompiled) in the same time as the postponed filters does 100,000 primes:
primes :: [Integer]
primes = 2: 3: sieve [] (tail primes) 5
where
notDivsBy d n = n `mod` d /= 0
sieve ds (p:ps) x = foldr (filter . notDivsBy) [x, x+2..p*p2] ds
++ sieve (p:ds) ps (p*p+2)
This one explicitly maintains the list of primes needed for testing each odds span between successive primes squares, which it also explicitly generates. But it tests with nested filter
s, which it repeatedly recreates.
Generated Spans, by List of Primes
The list of primes needed to test each range of odds is actually just the prefix of the primes list itself, of known length, and need not be specifically generated at all. Combined with onecall testing by the explicit list of primes, and direct generation of odds between the successive primes squares, this calculates only the bare essentials and makes no extra recalculations:
primes :: [Integer]
primes = 2: 3: sieve 0 (tail primes) 5
sieve k (p:ps) x = [x  x<[x,x+2..p*p2], and [x`rem`p/=0  p<fs]]
 or: all ((>0).(x`rem`)) fs
++ sieve (k+1) ps (p*p+2)
where fs = take k (tail primes)
It produces about 222,000 primes in the same amount of time, and is good for creating about a million first primes, compiled.
The reason to have sieve
function available separately too is that it can also be used to produce primes above a given number, as in
primesFrom m = sieve (length h) ps $ m`div`2*2+1
where
(h,(_:ps)) = span (<= (floor.sqrt.fromIntegral) m) primes
It can thus produce a few primes e.g. above 239812076741689
, which is a square of the millionth odd prime, without having to compute all the preceding primes (which would number in trillions).
Getting the Composite Numbers Out
The divisibility testing too should be considered a specification (as in no multiples of p), and not a code per se, because although testing composites is cheap (as most of them will have small factors, so the test is soon aborted), testing prime numbers is costly, and is to be avoided.
Euler's Sieve
With each found prime Euler's sieve removes all its multiples in advance so that at each step the list to process is guaranteed to not have any multiples of all the preceding primes in it (consists only of numbers coprime with all the preceding primes) and thus starts with the next prime:
import Data.OrdList (minus)
euler xs@(p:xt) = p : euler (xt `minus` map (*p) xs)
primes = euler [2..]
This code is extremely inefficient, running at about O() complexity, and should be regarded a specification only. It works very hard trying to avoid double hits on multiples, which are relatively few and cheap to deal with in the first place. Turner's sieve could also be seen as its rendition, substituting the builtin filter
for minus
 computationally different but having the same effect.
But the situation can be improved using a different list representation, generating them not from a last element and an increment, but rather a last span and an increment, which entails a set of helpful equivalences:
fromSpan (xs,i) = concat $ iterate (map (+ i)) xs
{ map (*p) (fromSpan (xs,i))
=== fromSpan (map (*p) xs, i*p)
fromSpan (xs,i) `minus` fromSpan (ys,i)
=== fromSpan (xs `minus` ys, i)
fromSpan (p:xt,i) === p : fromSpan (xt ++ [p+i], i)
fromSpan (xs,i) === forall (p > 0).
fromSpan (concat $ take p $ iterate (map (+ i)) xs, i*p) }
eulerS () = iterate eulerStep ([2],1)
eulerStep w@(xs@(p:_), i) = w'
where
(_:ys) = concat $ take p $ iterate (map (+ i)) xs
w' = ( (ys ++ [p+i*p]) `minus` map (*p) xs, i*p)
{ > mapM_ print $ take 4 $ eulerS ()
([2],1)
([3],2)
([5,7],6)
([7,11,13,17,19,23,29,31],30) }
This is where the notion of wheel comes from, naturally and elegantly. For each wheel specification w@((p:_),_)
produced by eulerStep
, the numbers in (fromSpan w)
up to are all primes too, so that
eulerPrimesTo m = if m > 1 then go ([2],1) else []
where
go w@((p:_), _)
 m < p*p = takeWhile (<= m) (fromSpan w)
 True = p : go (eulerStep w)
This runs at about O() complexity, for n
primes produced.
Postponed Multiples Removal
All the filtering versions thus far try to keep the primes among all numbers by testing each number in isolation. Instead, we can just remove the prime's multiples in advance. We gain in speed because we now get the primes for free, after all the multiples are removed on a particular span, without performing any divisibility tests at all, with this simple variation of the Postponed Filters sieve:
primes :: [Integer]
primes = 2: 3: sieve (tail primes) [5,7..]
where
sieve (p:ps) xs = h ++ sieve ps (t `minus` tail [q,q+2*p..])
where (h,~(_:t)) = span (< q) xs
q = p*p
This is similar to Euler's sieve in that it removes multiples of each prime in advance, but will generate some of the multiples more than once, like e.g. 3*15
and 5*9
, unlike the true Euler's sieve.
Merged Multiples Removal Sieve
The construct (...((sa)b)...)
in Postponed Multiples Removal is the same as (s(a+b+...))
, and so we can just remove the merged infinite primes multiples, each starting at its prime's square, from the initial numbers supply. This way we needn't explicitly jump to a prime's square because it's where its multiples start anyway:
import Data.OrdList (minus,union)
primes :: [Integer]
primes = 2:primes'
where
primes' = [3] ++ [5,7..] `minus` foldr union' [] mults
mults = map (\p>let q=p*p in (q,tail [q,q+2*p..])) $ primes'
union' (q,qs) xs = q : union qs xs
This code is yet faster. Its main deficiency still is that it creates a linear nested merging structure, instead of a treelike structure. Each multiple produced by a prime has to percolate to the top eventually, so it's better to make its path shorter. It'll have to go through fewer merge nodes this way.
The linearity is imposed by type asymmetry of our (union' :: a > b > b)
function, forcing us into the 1+(1+(1+(1+...)))
pattern, +
standing for union'
(which was defined that way to prevent the runahead when folding over the infinite list of lists of multiples).
We need to turn it into an associative operation of uniform type (:: a > a > a)
to be able to freely rearrange the combinations into arbitrary treelike patterns, as in e.g. ((1+1)+(2+2))+(...) etc. The type uniformity is what makes compositionality possible.
The code in the "Implicit Heap" section below improves on that, and is essentially equivalent to using a treefold instead of a standard linear foldr
, as in:
Treefold Merged Multiples Removal
primes :: [Integer]
primes = 2:primes'
where
primes' = [3,5,7,11,13,17] ++ [19,21..] `minus` comps
mults = map (\p> let q=p*p in ([q],tail [q,q+2*p..])) $ primes'
comps = fst $ tfold unionSP (pairwise unionSP mults)
pairwise f (x:y:ys) = f x y : pairwise f ys
tfold f (a:b:c:xs) = (a `f` (b `f` c)) `f` tfold f (pairwise f xs)
unionSP (a:as,bs) v = (a:x,y) where (x,y)=unionSP (as,bs) v
unionSP u@([],b:bs) v@(c:cs,ds) = case compare b c of
LT > (b:x,y) where (x,y)=unionSP ([],bs) v
EQ > (b:x,y) where (x,y)=unionSP ([],bs) (cs,ds)
GT > (c:x,y) where (x,y)=unionSP u (cs,ds)
unionSP ([],bs) ([],ds) = ([],union bs ds)
The fold used here creates a (2+(2+2))+( (4+(4+4)) + ( (8+(8+8)) + ... ))
structure, better adjusted for primes multiples production than 1+(2+(4+(8+...)))
, used by the "Implicit Heap" code, giving it additional 10% speedup.
unionSP
is an associative operation, preserving of the invariant such that for a list of multiples [(a,b),(c,d), ...]
, it's always the case that last a < head b && last a < head c
. These "split pairs" represent ordered list as a pair of its known and (currently) finite prefix, and the rest of it. Such pairs under unionSP
operation form a monoid, and if we were to declare a newtype SplitPair a = SP ([a],[a])
a Monoid
instance, with unionSP
its mappend
(and tfold
its mconcat
), the above code for comps
would just become SP (comps,_) = mconcat mults
.
This code exhibits approximately O()..O() local asymptotic behavior (tested interpreted, in GHCi, for 10,000 to 300,000 primes produced). When compared with Melissa O'Neill's PQ code from the ZIP package which was modified to work on odds only as well, it is 3.2x times faster, with used memory reported about 2.5x times smaller.
It can be further improved with the wheel optimization (as seen in Euler's Sieve and described below in Prime Wheels section):
Treefold Merged Multiples, with Wheel
primes :: [Integer]
primes = 2:3:5:7:primes'
where
primes' = [11,13,17,19,23,29] ++ (rollFrom 31) `minus` comps
mults = map (\p> fromList $ map (p*) $ rollFrom p) $ primes'
comps = fst $ tfold unionSP (pairwise unionSP mults)
fromList (x:xs) = ([x],xs)
rollFrom n = let x = (n11) `mod` 210
(y,_) = span (< x) wheelNums
in roll n $ drop (length y) wheel
wheelNums = roll 0 wheel
roll = scanl (+)
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
This runs about 2.5x times faster than the Priority Queuebased code as present in Melissa O'Neill's ZIP package, with similar local asymptotic behavior of about O()..O() (tested interpreted, in GHCi, for 10,000 to 300,000 primes produced), with used memory reported about twice less as well.
Multiples Removal on Generated Spans, or Sieve of Eratosthenes
Potponed multiples removal, with work done segmentwise, on the spans of numbers between the successive primes squares, becomes
import Data.OrdList (minus,union)
primes :: [Integer]
primes = 2: 3: sieve [] (tail primes) 5
where
sieve fs (p:ps) x = [x,x+2..q2] `minus` foldl union [] mults
++ sieve ((2*p,q):fs') ps (q+2)
where
q = p*p
mults = [ [y+s,y+2*s..q]  (s,y)< fs]
fs' = [ (s,last ms)  ((s,_),ms)< zip fs mults]
This modifies the Generated Spans sieve to "mark" the odd composites in a given range (instead of testing all the odds' divisibility) by generating them  just as a person performing the original sieve of Eratosthenes would do, marking one by one the multiples of the relevant primes. These composites are independently generated so some will be generated multiple times.
Interpreted under both GHCi and WinHugs, it runs faster, takes less memory, and has better asymptotic behavior, its performance approximately the same as in the Merged Multiples Removal sieve. The advantage in working with spans explicitly is that this code is easily amendable to using arrays for the composites marking and removal on each finite span.
Using Immutable Arrays
Generating Segments of Primes
The removal of multiples on each segment of odds in the Sieve of Eratosthenes can be done by actually marking them in an array, instead of manipulating the ordered lists:
primes :: [Integer]
primes = 2: 3: sieve [] (tail primes) 5
where
sieve fs (p:ps) x = [i  i< [x,x+2..q2], a!i]
++ sieve ((2*p,q):fs') ps (q+2)
where
q = p*p
mults = [ [y+s,y+2*s..q]  (s,y)< fs]
fs' = [ (s,last ms)  ((s,_),ms)< zip fs mults]
a = accumArray (\a b>False) True (x,q2)
[(i,())  ms< mults, i< ms]
Apparently, arrays are fast too. This code, compared with Treefold Merged with Wheel version (itself 2.5x times faster than Melissa O'Neill's PQ version), runs at about the same time and memory usage, but improving slightly with slightly better local asymptotics.
Calculating Primes Upto a Given Value
primesToN n = 2: [i  i<[3,5..n], ar!i]
where
ar = f 5 $ accumArray (\a b>False) True (3,n)
[(i,())  i< [9,15..n]]
f p a  q > n = a
 True = if null x then a' else f (head x) a'
where q = p*p
a'= a//[(i,False)i<[q,q+2*p..n]]
x = [i  i<[p+2,p+4..n], a' !i]
Calculating Primes in a Given Range
primesFromTo a b = (if a<3 then [2] else [])
++ [i  i<[o,o+2..b], ar!i]
where
o = max (if even a then a+1 else a) 3
r = floor.sqrt.fromInteger $ b+1
ar = accumArray (\a b> False) True (o,b)
[(i,())  p < [3,5..r]
, let q = p*p
s = 2*p
(n,x) = quotRem (o  q) s
q' = if o <= q then q
else q + (n + signum x)*s
, i< [q',q'+s..b] ]
Although testing by odds instead of by primes, the array generation is so fast that it is very much feasible and even preferable for quick generation of some short spans of relatively big primes.
Using Mutable Arrays
Using mutable arrays is the fastest and most memory efficient way to calculate prime numbers in Haskell.
Using ST Array
This method implements the Sieve of Eratosthenes, similar to how you might do it in C.
import Control.Monad
import Control.Monad.ST
import Data.Array.IArray
import Data.Array.MArray
import Data.Array.ST
import Data.Array.Unboxed
primesToNA :: Int > UArray Int Bool
primesToNA n = runSTUArray sieve where
sieve = do
a < newArray (2,n) True :: ST s (STUArray s Int Bool)
let sr = floor . (sqrt::Double>Double) . fromIntegral $ n+1
forM_ [4,6..n] $ \j > writeArray a j False
forM_ [3,5..sr] $ \i > do
si < readArray a i
when si $
forM_ [i*i,i*i+i+i..n] $ \j > writeArray a j False
return a
primesToN :: Int > [Int]
primesToN n = [i  (i,e) < assocs . primesToNA $n, e]
Bitwise prime sieve with Template Haskell
Count the number of prime below a given 'n'. Shows fast bitwise arrays, and an example of Template Haskell to defeat your enemies.
{# OPTIONS O2 optcO XBangPatterns #}
module Primes (nthPrime) where
import Control.Monad.ST
import Data.Array.ST
import Data.Array.Base
import System
import Control.Monad
import Data.Bits
nthPrime :: Int > Int
nthPrime n = runST (sieve n)
sieve n = do
a < newArray (3,n) True :: ST s (STUArray s Int Bool)
let cutoff = truncate (sqrt $ fromIntegral n) + 1
go a n cutoff 3 1
go !a !m cutoff !n !c
 n >= m = return c
 otherwise = do
e < unsafeRead a n
if e then
if n < cutoff then
let loop !j
 j < m = do
x < unsafeRead a j
when x $ unsafeWrite a j False
loop (j+n)
 otherwise = go a m cutoff (n+2) (c+1)
in loop ( if n < 46340 then n * n else n `shiftL` 1)
else go a m cutoff (n+2) (c+1)
else go a m cutoff (n+2) c
And places in a module:
{# OPTIONS fth #}
import Primes
main = print $( let x = nthPrime 10000000 in [ x ] )
Run as:
$ ghc make o primes Main.hs
$ time ./primes
664579
./primes 0.00s user 0.01s system 228% cpu 0.003 total
Implicit Heap
The following is a more efficient prime generator, implementing the sieve of Eratosthenes.
See also the message threads Re: "nocoding" functional data structures via lazyness for more about how merging ordered lists amounts to creating an implicit heap and Re: Code and Perf. Data for Prime Finders for an explanation of the People a
structure that makes it work when tying a knot.
data People a = VIP a (People a)  Crowd [a]
mergeP :: Ord a => People a > People a > People a
mergeP (VIP x xt) ys = VIP x $ mergeP xt ys
mergeP (Crowd xs) (Crowd ys) = Crowd $ merge xs ys
mergeP xs@(Crowd (x:xt)) ys@(VIP y yt) = case compare x y of
LT > VIP x $ mergeP (Crowd xt) ys
EQ > VIP x $ mergeP (Crowd xt) yt
GT > VIP y $ mergeP xs yt
merge :: Ord a => [a] > [a] > [a]
merge xs@(x:xt) ys@(y:yt) = case compare x y of
LT > x : merge xt ys
EQ > x : merge xt yt
GT > y : merge xs yt
diff xs@(x:xt) ys@(y:yt) = case compare x y of
LT > x : diff xt ys
EQ > diff xt yt
GT > diff xs yt
foldTree :: (a > a > a) > [a] > a
foldTree f ~(x:xs) = x `f` foldTree f (pairs xs)
where pairs ~(x: ~(y:ys)) = f x y : pairs ys
primes, nonprimes :: [Integer]
primes = 2:3:diff [5,7..] nonprimes
nonprimes = serve . foldTree mergeP . map multiples $ tail primes
where
multiples p = vip [p*p,p*p+2*p..]
vip (x:xs) = VIP x $ Crowd xs
serve (VIP x xs) = x:serve xs
serve (Crowd xs) = xs
nonprimes
effectively implements a heap, exploiting lazy evaluation.
Prime Wheels
The idea of only testing odd numbers can be extended further. For instance, it is a useful fact that every prime number other than 2 and 3 must be of the form or . Thus, we only need to test these numbers:
primes :: [Integer]
primes = 2:3:primes'
where
1:p:candidates = [6*k+r  k < [0..], r < [1,5]]
primes' = p : filter isPrime candidates
isPrime n = all (not . divides n)
$ takeWhile (\p > p*p <= n) primes'
divides n p = n `mod` p == 0
Here, primes'
is the list of primes greater than 3 and isPrime
does not test for divisibility by 2 or 3 because the candidates
by construction don't have these numbers as factors. We also need to exclude 1 from the candidates and mark the next one as prime to start the recursion.
Such a scheme to generate candidate numbers first that avoid a given set of primes as divisors is called a prime wheel. Imagine that you had a wheel of circumference 6 to be rolled along the number line. With spikes positioned 1 and 5 units around the circumference, rolling the wheel will prick holes exactly in those positions on the line whose numbers are not divisible by 2 and 3.
A wheel can be represented by its circumference and the spiked positions.
data Wheel = Wheel Integer [Integer]
We prick out numbers by rolling the wheel.
roll (Wheel n rs) = [n*k+r  k < [0..], r < rs]
The smallest wheel is the unit wheel with one spike, it will prick out every number.
w0 = Wheel 1 [1]
We can create a larger wheel by rolling a smaller wheel of circumference n
along a rim of circumference p*n
while excluding spike positions at multiples of p
.
nextSize (Wheel n rs) p =
Wheel (p*n) [r'  k < [0..(p1)], r < rs,
let r' = n*k+r, r' `mod` p /= 0]
Combining both, we can make wheels that prick out numbers that avoid a given list ds
of divisors.
mkWheel ds = foldl nextSize w0 ds
Now, we can generate prime numbers with a wheel that for instance avoids all multiples of 2, 3, 5 and 7.
primes :: [Integer]
primes = small ++ large
where
1:p:candidates = roll $ mkWheel small
small = [2,3,5,7]
large = p : filter isPrime candidates
isPrime n = all (not . divides n)
$ takeWhile (\p > p*p <= n) large
divides n p = n `mod` p == 0
It's a pretty big wheel with a circumference of 210 and allows us to calculate the first 10000 primes in convenient time.
A fixed size wheel is fine, but how about adapting the wheel size while generating prime numbers? See Euler's Sieve above, or the functional pearl titled Lazy wheel sieves and spirals of primes for more.
Using IntSet for a traditional sieve
module Sieve where
import qualified Data.IntSet as I
 findNext  finds the next member of an IntSet.
findNext c is  I.member c is = c
 c > I.findMax is = error "Ooops. No next number in set."
 otherwise = findNext (c+1) is
 mark  delete all multiples of n from n*n to the end of the set
mark n is = is I.\\ (I.fromAscList (takeWhile (<=end) (map (n*) [n..])))
where
end = I.findMax is
 primes  gives all primes up to n
primes n = worker 2 (I.fromAscList [2..n])
where
worker x is
 (x*x) > n = is
 otherwise = worker (findNext (x+1) is) (mark x is)
Testing Primality
Primality Test and Integer Factorization
Given an infinite list of prime numbers, we can implement primality tests and integer factorization:
isPrime n = n > 1 && n == head (primeFactors n)
primeFactors 1 = []
primeFactors n = go n primes
where
go n ps@(p:pt)
 p*p > n = [n]
 n `rem` p == 0 = p : go (n `quot` p) ps
 otherwise = go n pt
MillerRabin Primality Test
find2km :: Integral a => a > (a,a)
find2km n = f 0 n
where
f k m
 r == 1 = (k,m)
 otherwise = f (k+1) q
where (q,r) = quotRem m 2
millerRabinPrimality :: Integer > Integer > Bool
millerRabinPrimality n a
 a <= 1  a >= n1 =
error $ "millerRabinPrimality: a out of range ("
++ show a ++ " for "++ show n ++ ")"
 n < 2 = False
 even n = False
 b0 == 1  b0 == n' = True
 otherwise = iter (tail b)
where
n' = n1
(k,m) = find2km n'
b0 = powMod n a m
b = take (fromIntegral k) $ iterate (squareMod n) b0
iter [] = False
iter (x:xs)
 x == 1 = False
 x == n' = True
 otherwise = iter xs
pow' :: (Num a, Integral b) => (a > a > a) > (a > a) > a > b > a
pow' _ _ _ 0 = 1
pow' mul sq x' n' = f x' n' 1
where
f x n y
 n == 1 = x `mul` y
 r == 0 = f x2 q y
 otherwise = f x2 q (x `mul` y)
where
(q,r) = quotRem n 2
x2 = sq x
mulMod :: Integral a => a > a > a > a
mulMod a b c = (b * c) `mod` a
squareMod :: Integral a => a > a > a
squareMod a b = (b * b) `rem` a
powMod :: Integral a => a > a > a > a
powMod m = pow' (mulMod m) (squareMod m)
External links
 A collection of prime generators; the file "ONeillPrimes.hs" contains one of the fastest pureHaskell prime generators. WARNING: Don't use the priority queue from that file for your projects: it's broken and works for primes only by a lucky chance.