# Testing primality

### From HaskellWiki

(Difference between revisions)

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

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

Line 1: | Line 1: | ||

− | == Testing Primality == |
+ | = 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 === |
+ | == Primality Test and Integer Factorization == |

Given an infinite list of prime numbers, we can implement primality tests and integer factorization: |
Given an infinite list of prime numbers, we can implement primality tests and integer factorization: |
||

<haskell> |
<haskell> |
||

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

− | 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 1 = [] |
+ | primeFactors n | n > 1 = go n primes |

− | primeFactors n = go n primes |
+ | where |

− | where |
+ | go n ps@(p:ps') |

− | go n ps@(p:pt) |
||

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

− | | n `rem` p == 0 = p : go (n `quot` p) ps |
+ | | n `rem` p == 0 = p : go (n `quot` p) ps |

− | | otherwise = go n pt |
+ | | otherwise = go n ps' |

+ | </haskell> |
||

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

+ | <haskell> |
||

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

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

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

<haskell> |
<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 30: | Line 31: | ||

| 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 50: | Line 51: | ||

| 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 67: | Line 68: | ||

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 08:13, 17 August 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 == head (primeFactors 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 where go n ps@(p:ps') | p*p > n = [n] | n `rem` p == 0 = p : go (n `quot` p) ps | otherwise = go n ps'

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]