# Testing primality

### From HaskellWiki

(Difference between revisions)

m (added "mathematics"-category) |
(→Primality Test and Integer Factorization: faster isPrime) |
||

Line 8: | Line 8: | ||

<haskell> |
<haskell> |
||

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

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

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

+ | True primes |
||

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

## Revision as of 11:27, 5 June 2011

## 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 tests and integer factorization:

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

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