Today’s Programming Praxis problem is about checking whether or not a number is prime. We’re supposed to implement the Miller-Rabin primality test, which is probably the most well-known algorithm to do this. Let’s go.
First the imports:
import Control.Arrow import Data.Bits import Data.List import System.Random
The Miller-Rabin algorithm calls for random numbers. Since I didn’t want to make isPrime an IO function, we have to pass it a random number generator as well.
isPrime :: Integer -> StdGen -> Bool isPrime n g = let (s, d) = (length *** head) . span even $ iterate (flip div 2) (n - 1) xs = map (expm n d) . take 50 $ randomRs (2, n - 2) g in all (\x -> elem x [1, n - 1] || any (== n - 1) (take s $ iterate (expm n 2) x)) xs
You may have noticed the expm in there. In theory all you have to do is mod (random_number ^ d) n. However, in our test case of 2^89 – 1, both the random number and d can easily be 10 digits. And 10-digit-number ^ 10-digit-number is a very big number, which takes a while to calculate. As a result, your algorithm is going to be terribly slow, as I discovered in my initial attempt. The solution lies in using modular exponentiation. This is a much faster way of calculating mod (a ^ b) m. In Scheme this function is apparently in the standard Prelude, but I couldn’t find a Haskell library that has it defined. No matter, we’ll just make our own:
expm :: Integer -> Integer -> Integer -> Integer expm m e b = foldl' (\r (b', _) -> mod (r * b') m) 1 . filter (flip testBit 0 . snd) . zip (iterate (flip mod m . (^ 2)) b) $ takeWhile (> 0) $ iterate (flip shiftR 1) e
And there’s our primality checking code. To test it, use:
main :: IO () main = print . isPrime (2 ^ 89 - 1) =<< getStdGen