# Testing primality

### From HaskellWiki

(Difference between revisions)

(→Primality Test and Integer Factorization) |
(→Primality Test and Integer Factorization: add one-liner 'factors') |
||

(5 intermediate revisions by 2 users not shown) | |||

Line 1: | Line 1: | ||

− | + | = Testing Primality = | |

(for a context to this see [[Prime_numbers | Prime numbers]]). | (for a context to this see [[Prime_numbers | Prime numbers]]). | ||

− | + | == 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 n = unfoldr (\(d, n) -> listToMaybe [(x, (x, div n x)) | ||

+ | | n > 1, x <- [d..isqrt n] ++ [n], rem n x == 0]) (2,n) | ||

+ | |||

+ | isPrime n = factors n == [n] | ||

+ | isqrt n = floor . sqrt . fromIntegral $ n | ||

+ | </haskell> | ||

+ | |||

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

<haskell> | <haskell> | ||

− | + | 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 -- or go n (2:[3,5..]) | |

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

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

− | go n ps@(p: | + | | p*p > n = [n] |

− | | p*p > n | + | | r == 0 = p : go q ps |

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

− | | otherwise | + | where |

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

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

− | + | When trying to factorize only one number or two, it will be faster to just use <code>(2:[3,5..])</code> 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 | |

<haskell> | <haskell> | ||

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

+ | </haskell> | ||

+ | |||

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

+ | |||

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

+ | <haskell> | ||

+ | -- (eq. to) find2km (2^k * n) = (k,n) | ||

find2km :: Integral a => a -> (a,a) | find2km :: Integral a => a -> (a,a) | ||

find2km n = f 0 n | find2km n = f 0 n | ||

Line 31: | Line 51: | ||

| otherwise = f (k+1) q | | otherwise = f (k+1) q | ||

where (q,r) = quotRem m 2 | where (q,r) = quotRem m 2 | ||

− | + | ||

+ | -- n is the number to test; a is the (presumably randomly chosen) witness | ||

millerRabinPrimality :: Integer -> Integer -> Bool | millerRabinPrimality :: Integer -> Integer -> Bool | ||

millerRabinPrimality n a | millerRabinPrimality n a | ||

Line 51: | Line 72: | ||

| x == n' = True | | x == n' = True | ||

| otherwise = iter xs | | 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' :: (Num a, Integral b) => (a->a->a) -> (a->a) -> a -> b -> a | ||

pow' _ _ _ 0 = 1 | pow' _ _ _ 0 = 1 | ||

Line 68: | Line 90: | ||

squareMod :: Integral a => a -> a -> a | squareMod :: Integral a => a -> a -> a | ||

squareMod a b = (b * b) `rem` 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 :: Integral a => a -> a -> a -> a | ||

powMod m = pow' (mulMod m) (squareMod m) | powMod m = pow' (mulMod m) (squareMod m) | ||

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

+ | Example: | ||

+ | |||

+ | <haskell>-- check if '1212121' is prime with several witnesses | ||

+ | > map (millerRabinPrimality 1212121) [5432,1265,87532,8765,26] | ||

+ | [True,True,True,True,True] | ||

+ | </haskell> | ||

[[Category:Mathematics]] | [[Category:Mathematics]] |

## Revision as of 09:22, 3 June 2014

# 1 Testing Primality

(for a context to this see Prime numbers).

## 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 n = unfoldr (\(d, n) -> listToMaybe [(x, (x, div n x)) | n > 1, x <- [d..isqrt n] ++ [n], rem n x == 0]) (2,n) isPrime n = factors n == [n] isqrt n = floor . sqrt . fromIntegral $ n

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

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.

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