# Difference between revisions of "Prime numbers"

(Added section "External links", with a link to haskell-primes.zip) |
m |
||

Line 13: | Line 13: | ||

<haskell> | <haskell> | ||

− | isPrime n = n > 1 && n == head ( | + | isPrime n = n > 1 && n == head (primeFactors n) |

primeFactors 1 = [] | primeFactors 1 = [] |

## Revision as of 17:36, 6 December 2008

## Contents

## 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 (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
```

## 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 that are larger than .

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

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

## External links

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