Difference between revisions of "Testing primality"
(→Primality Test and Integer Factorization: reduce code width) 

Line 12:  Line 12:  
factors :: Integer > [Integer] 
factors :: Integer > [Integer] 

factors n 
factors n 

−  = unfoldr (\n > listToMaybe 
+  = unfoldr (\n > listToMaybe 
−  +  [(x, div n x )  x < [2..n], mod n x==0]) n 

−  = unfoldr (\(d,n) > listToMaybe 
+  = unfoldr (\(d,n) > listToMaybe 
−  +  [(x, (x, div n x))  x < [d..n], mod n x==0]) (2,n) 

+  = unfoldr (\(d,n) > listToMaybe 

+  [(x, (x, div n x))  x < takeWhile ((<=n).(^2)) 

+  [d..] ++ [nn>1], mod n x==0]) (2,n) 

isPrime n = n > 1 && head (factors n) == n 
isPrime n = n > 1 && head (factors n) == n 

</haskell> 
</haskell> 

Line 23:  Line 26:  
Of course there's no need to try any even numbers above 2. Given an infinite list of primes we can avoid ''any'' composites, not just evens: 
Of course there's no need to try any even numbers above 2. Given an infinite list of primes we can avoid ''any'' composites, not just evens: 

<haskell> 
<haskell> 

−  pfactors 
+  pfactors primes n = unfoldr (\(ds, n) > listToMaybe 
−  +  [(x, (dropWhile (< x) ds, div n x)) 

−  +   x < takeWhile ((<=n).(^2)) ds ++ [nn>1], mod n x==0]) 

+  (primes,n) 

primes :: [Integer] 
primes :: [Integer] 

−  primes = 2 : 3 : [x  x < [5,7.. 
+  primes = 2 : 3 : [x  x < [5,7..] 
+  , head (pfactors (tail primes) x) == x] 

</haskell> 
</haskell> 

Line 34:  Line 39:  
<haskell> 
<haskell> 

isPrime n = n > 1 && 
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..]) 
primeFactors n  n > 1 = go n primes  or go n (2:[3,5..]) 
Revision as of 17:01, 18 July 2022
Testing Primality
(for a context to this see Prime numbers).
Primality Test and Integer Factorization
Simplest primality test and integer factorization is by trial division:
import Data.List (unfoldr)
import Data.Maybe (listToMaybe)
factors :: Integer > [Integer]
factors n
= unfoldr (\n > listToMaybe
[(x, div n x )  x < [2..n], mod n x==0]) n
= unfoldr (\(d,n) > listToMaybe
[(x, (x, div n x))  x < [d..n], mod n x==0]) (2,n)
= unfoldr (\(d,n) > listToMaybe
[(x, (x, div n x))  x < takeWhile ((<=n).(^2))
[d..] ++ [nn>1], mod n x==0]) (2,n)
isPrime n = n > 1 && head (factors n) == n
The factors produced by this code are all prime by construction, because we enumerate possible divisors in ascending order while dividing each found factor out of the number being tested.
Of course there's no need to try any even numbers above 2. Given an infinite list of primes we can avoid any composites, not just evens:
pfactors primes n = unfoldr (\(ds, n) > listToMaybe
[(x, (dropWhile (< x) ds, div n x))
 x < takeWhile ((<=n).(^2)) ds ++ [nn>1], mod n x==0])
(primes,n)
primes :: [Integer]
primes = 2 : 3 : [x  x < [5,7..]
, head (pfactors (tail primes) x) == x]
Rewriting the above as a recursive code, we get:
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 oneoff 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 might be faster to just use (2:[3,5..])
as a source of possible divisors instead of calculating the prime numbers first, depending on the speed of your primes generator. 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.
MillerRabin 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 >= n1 =
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' = n1
(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]