Testing primality
Revision as of 11:36, 1 June 2014 by WillNess (talk | contribs) (→Primality Test and Integer Factorization: mention one-off invocation, on odds.)
The printable version is no longer supported and may have rendering errors. Please update your browser bookmarks and please use the default browser print function instead.
Testing Primality
(for a context to this see Prime numbers).
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..]
More at Prime numbers#Optimal trial division.
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]