# Prime numbers

## Simple Prime Sieve

The following is an elegant way to generate a list of all the prime numbers in the universe:

  primes :: [Integer]
primes = sieve [2..]
where sieve (p:xs) = p : sieve [x | x<-xs, x mod p /= 0]


While simple, this method is rather inefficient and not recommended for more than the first 1000 prime numbers. For every number x, it will test x mod p /= 0 for all prime numbers p that are smaller than the smallest prime factor of x.

Given an infinite list of prime numbers, we can implement primality tests and integer factorization:

  isPrime n = n > 1 && n == head (factorize 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


## Simple Prime Sieve II

The following method is slightly faster and works well for the first 5000 primes:

  primes :: [Integer]
primes = 2:filter isPrime [3,5..]
where
isPrime n   = all (not . divides n) $takeWhile (\p -> p*p <= n) primes divides n p = n mod p == 0  Compared to the previous sieve, it only tests odd numbers and avoids testing for prime factors of $n$ that are larger than $\sqrt{n}$. ## 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 that avoid the first 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 Int [Int]


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 Data.Array.ST
import Data.Array.Base
import System
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  ## 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)


## 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)