Testing primality
Revision as of 07:04, 4 July 2011 by WillNess (talk | contribs) (→Primality Test and Integer Factorization: extra parentheses added for clarity)
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 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..]
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)