# Prime numbers

### From HaskellWiki

## Contents |

## 1 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]

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

## 2 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 .

## 3 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 6*k* + 1 or 6*k* + 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

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]

nextSize (Wheel n rs) p = Wheel (p*n) [r' | k <- [0..(p-1)], r <- rs, let r' = n*k+r, r' `mod` p /= 0]

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.

## 4 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 thedata 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

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

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

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

## 8 External links

- A collection of prime generators; the file "ONeillPrimes.hs" contains one of the fastest pure-Haskell prime generators.