Prime numbers

Prime Number Resources

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.

Prime Numbers at Wikipedia.

Sieve of Eratosthenes 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 then 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 Simple Primes Sieve

Attributed to David Turner (SASL Language Manual, 1975), the following is a direct translation of that idea, generating a list of all the 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

While appearing simple, this method is extremely inefficient and not recommended for more than a few thousand prime numbers, even when compiled. Here, every number is tested for divisibility by each prime up to its smallest factor, so any prime is checked against all its preceding primes, when in fact 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 (upto 89's filter only).

So this code in effect creates a cascade of nested filters in front of the infinite numbers supply, and in an extremely premature fashion at that. One way of fixing that would be to postpone the creation of filters until the right moment, as in

Postponed Filters Sieve

  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 only tests odd numbers, and by only 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 the first version exhibits near O(${n^{2}}$) behavior, this one exhibits near O(${n^{1.5}}$) behavior, and is good for creating a few 100,000s of primes (compiled). It takes as much time to generate 100,000 primes with it as it takes for the first one to generate 5500, GHC-compiled.

Postponed Multiples Removal i.e. Euler's Sieve

Instead of testing each number for divisibility by a prime we can just remove the prime's multiples in advance to achieve the same effect. We gain in that we now get the primes for free, after all the multiples are removed on a particular span, without performing any divisibility tests at all:

  primes :: [Integer]
  primes = 2: 3: sieve (tail primes) [5,7..]  
   where
    sieve (p:ps) xs = h ++ sieve ps (t `minus` [p*p+2*p,p*p+4*p..])
                            -- was: [x | x<-t, x `rem` p /= 0]
                      where (h,~(_:t)) = span (< p*p) xs 
    minus xs@(x:xt) ys@(y:yt)
            | x<y   = x: xt `minus` ys
            | x>y   =    xs `minus` yt
            | True  =    xt `minus` yt

This is, in fact, Euler's sieve.

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

More Filtering Sieves

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").

Odd numbers, by Trial Division

This is also good for generating a few 100,000s primes (when GHC-compiled 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 re-fetches 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

This version is a bit faster still, creating 158,000 primes (again, GHC-compiled) in the same time span (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*p-2] 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 filters, 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 one-call testing by the explicit list of primes, and direct generation of odds between the successive primes squares, this leads to:

  primes :: [Integer]
  primes = 2: 3: sieve 0 (tail primes) 5
  sieve k (p:ps) x = [x | x<-[x,x+2..p*p-2], 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 time span, 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).

Sieve of Eratosthenes

All the divisibility-testing versions thus far try to keep the primes among all numbers by testing each number in isolation. Testing composites is cheap (as most of them will have small factors, so the test is soon aborted), but testing prime numbers is costly. The sieve of Eratosthenes gets rid of composites from the spans of numbers and thus gets its primes "for free", as gaps between the generated/crossed-off/marked composites:

  primes :: [Integer]
  primes = 2: 3: 5: sieve [] (tail primes) 7
   where
    sieve fs (p:ps) x = gaps x comps
                        ++ sieve ((2*p,q+2*p):fs') ps (q+2)
     where
      q                    = p*p
      comps                = foldl mrg [] mults
      (mults,fs')          = un_zip $ map mark fs
      mark (s,x)           = (h,(s,x')) 
                               where (h,~(x':_)) = span(<q)[x,x+s..]
      gaps x (y:ys) | x==y = gaps (x+2) ys 
                    | True = x: gaps (x+2) (y:ys)
      gaps x []            = if x==q then [] else [x]   
      mrg u@(a:t) w@(b:r) 
                    | a<b  = a: t `mrg` w
                    | b<a  = b: u `mrg` r  
                    | True = a: t `mrg` r
      mrg u       w        = if null u then w else u

This modifies the preceding sieve to "mark" the odd composites in a given range (instead of testing their 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 - and then look for gaps in the composites list.

Compared with the previous version, if compiled with GHC, it is steadily 1.2x slower and takes 3.5x more memory. But interpreted under both GHCi and WinHugs, it runs faster, takes less memory, and has better asymptotic behavior.

Merged Multiples Sieve

The explicit partitioning into spans between the successive primes squares in both the Sieve of Eratosthenes and Euler's Sieve versions above is actually superfluous and is implicitly achieved when looking for gaps in the directly merged infinite primes multiples streams (as if turning (...((s-a)-b)-...) into (s-(a+b+...)) in the Euler's Sieve definition):

  primes :: [Integer]
  primes = [2,3] ++ gaps 5 
             (foldr (\p-> let q=p*p in (q:).merge [q+2*p,q+4*p..]) 
                    [] $ tail primes)
   where
    gaps x s@(y:ys) | x==y  = gaps (x+2) ys
                    | True  = x: gaps (x+2) s
    merge a@(x:xs) b@(y:ys)
                    | x > y = y: merge a  ys
                    | x < y = x: merge xs b
                    | True  = x: merge xs ys

This code is the fastest of all presented so far. Here the specific clause order used in merge achieves additional 30% - 50% speedup.

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 $6k+1$ or $6k+5$. 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..(p-1)], 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 the functional pearl titled Lazy wheel sieves and spirals of primes for more.

Implicit Heap

The following is a more efficient prime generator, implementing the sieve of Eratosthenes.

See also the message threads Re: "no-coding" 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) = f x . 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*k | k <- [p,p+2..]]

    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.

Bitwise prime sieve

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 -optc-O -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

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

Miller-Rabin 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 >= n-1 = 
        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' = n-1
        (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

• http://www.cs.hmc.edu/~oneill/code/haskell-primes.zip

A collection of prime generators; the file "ONeillPrimes.hs" contains one of the fastest pure-Haskell 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.