# Testing primality

(Difference between revisions)

# 1 Testing Primality

(for a context to this see Prime numbers).

## 1.1 Primality Test and Integer Factorization

Given an infinite list of prime numbers, we can implement primality test and integer factorization by trial division:

```-- isPrime n = head (primeFactors n) == n
isPrime n = n > 1 &&
foldr (\p r -> p*p > n || ((n `rem` p) /= 0 && r))
True primes

primeFactors n | n > 1 = go n primes   -- or go n (2:[3,5..])
where                               -- for one-off invocation
go n ps@(p:t)
| p*p > n    = [n]
| r == 0     =  p : go q ps
| otherwise  =      go n t
where
(q,r) = quotRem n p```

When trying to factorize only one number or two, it will be faster to just use `(2:[3,5..])` as a source of possible divisors instead of first finding prime numbers and then using them. For more than a few factorizations, when no other primes source is available, just use

`primes = 2 : filter isPrime [3,5..]`

## 1.2 Miller-Rabin Primality Test

```-- (eq. to) find2km (2^k * n) = (k,n)
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

-- n is the number to test; a is the (presumably randomly chosen) witness
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

-- (eq. to) pow' (*) (^2) n k = n^k
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

-- (eq. to) powMod m n k = n^k `mod` m
powMod :: Integral a => a -> a -> a -> a
powMod m = pow' (mulMod m) (squareMod m)```

Example:

```-- check if '1212121' is prime with several witnesses
> map (millerRabinPrimality 1212121) [5432,1265,87532,8765,26]
[True,True,True,True,True]```