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

Share this: Twitter

Facebook

Like this: Like Loading... Related

Tags: prime primality checking programming praxis kata