# Testing primality

### From HaskellWiki

JaimeSoffer (Talk | contribs) (documentation of Miller-Rabin test) |
(→Primality Test and Integer Factorization) |
||

(3 intermediate revisions by one user not shown) | |||

Line 5: | Line 5: | ||

== Primality Test and Integer Factorization == | == Primality Test and Integer Factorization == | ||

− | Given an infinite list of | + | Simplest primality test and integer factorization is by trial division: |

+ | <haskell> | ||

+ | import Data.List (unfoldr) | ||

+ | import Data.Maybe (listToMaybe) | ||

+ | |||

+ | factors :: Integer -> [Integer] | ||

+ | factors n | ||

+ | = unfoldr (\n -> listToMaybe [(x, div n x) | x <- [2..n], mod n x==0]) n | ||

+ | = unfoldr (\(d,n) -> listToMaybe [(x, (x, div n x)) | x <- [d..n], mod n x==0]) (2,n) | ||

+ | = unfoldr (\(d,n) -> listToMaybe [(x, (x, div n x)) | x <- takeWhile ((<=n).(^2)) | ||

+ | [d..] ++ [n|n>1], mod n x==0]) (2,n) | ||

+ | isPrime n = n > 1 && head (factors n) == n | ||

+ | </haskell> | ||

+ | |||

+ | The factors produced by this code are all prime by construction, because we enumerate possible divisors in ascending order while dividing each found factor out of the number being tested. | ||

+ | |||

+ | Of course there's no need to try any even numbers above 2. Given an infinite list of primes we can avoid ''any'' composites, not just evens: | ||

+ | <haskell> | ||

+ | pfactors prs n = unfoldr (\(ds,n) -> listToMaybe | ||

+ | [(x, (dropWhile (< x) ds, div n x)) | x <- takeWhile ((<=n).(^2)) ds ++ | ||

+ | [n|n>1], mod n x==0]) (prs,n) | ||

+ | primes :: [Integer] | ||

+ | primes = 2 : 3 : [x | x <- [5,7..], head (pfactors (tail primes) x) == x] | ||

+ | </haskell> | ||

+ | |||

+ | Re-writing the above as a recursive code, we get: | ||

<haskell> | <haskell> | ||

− | |||

isPrime n = n > 1 && | isPrime n = n > 1 && | ||

foldr (\p r -> p*p > n || ((n `rem` p) /= 0 && r)) | foldr (\p r -> p*p > n || ((n `rem` p) /= 0 && r)) | ||

True primes | True primes | ||

− | primeFactors n | n > 1 = go n primes | + | primeFactors n | n > 1 = go n primes -- or go n (2:[3,5..]) |

− | where | + | where -- for one-off invocation |

− | go n ps@(p: | + | go n ps@(p:t) |

− | | p*p > n | + | | p*p > n = [n] |

− | | | + | | r == 0 = p : go q ps |

− | | otherwise | + | | otherwise = go n t |

+ | where | ||

+ | (q,r) = quotRem n p | ||

</haskell> | </haskell> | ||

− | When no other primes source is available, just use | + | |

+ | When trying to factorize only one number or two, it might be faster to just use <code>(2:[3,5..])</code> as a source of possible divisors instead of calculating the prime numbers first, depending on the speed of your primes generator. For more than a few factorizations, when no other primes source is available, just use | ||

<haskell> | <haskell> | ||

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

</haskell> | </haskell> | ||

+ | |||

+ | More at [[Prime numbers#Optimal trial division]]. | ||

== Miller-Rabin Primality Test == | == Miller-Rabin Primality Test == |

## Latest revision as of 02:53, 23 November 2015

# [edit] 1 Testing Primality

(for a context to this see Prime numbers).

## [edit] 1.1 Primality Test and Integer Factorization

Simplest primality test and integer factorization is by trial division:

import Data.List (unfoldr) import Data.Maybe (listToMaybe) factors :: Integer -> [Integer] factors n = unfoldr (\n -> listToMaybe [(x, div n x) | x <- [2..n], mod n x==0]) n = unfoldr (\(d,n) -> listToMaybe [(x, (x, div n x)) | x <- [d..n], mod n x==0]) (2,n) = unfoldr (\(d,n) -> listToMaybe [(x, (x, div n x)) | x <- takeWhile ((<=n).(^2)) [d..] ++ [n|n>1], mod n x==0]) (2,n) isPrime n = n > 1 && head (factors n) == n

The factors produced by this code are all prime by construction, because we enumerate possible divisors in ascending order while dividing each found factor out of the number being tested.

Of course there's no need to try any even numbers above 2. Given an infinite list of primes we can avoid *any* composites, not just evens:

pfactors prs n = unfoldr (\(ds,n) -> listToMaybe [(x, (dropWhile (< x) ds, div n x)) | x <- takeWhile ((<=n).(^2)) ds ++ [n|n>1], mod n x==0]) (prs,n) primes :: [Integer] primes = 2 : 3 : [x | x <- [5,7..], head (pfactors (tail primes) x) == x]

Re-writing the above as a recursive code, we get:

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 might be faster to just use `(2:[3,5..])`

as a source of possible divisors instead of calculating the prime numbers first, depending on the speed of your primes generator. 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.

## [edit] 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]