Personal tools

Testing primality

From HaskellWiki

(Difference between revisions)
Jump to: navigation, search
(→Primality Test and Integer Factorization)
Current revision (08:13, 17 August 2011) (edit) (undo)
(documentation of Miller-Rabin test)
 
(3 intermediate revisions 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 | n > 1 = 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 31: Line 35:
| 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 56:
| 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 74:
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]]

Current revision

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]
Retrieved from "http://www.haskell.org/haskellwiki/Testing_primality"

Category: Mathematics