Rabin

Rabin.hs

isPrime :: Int -> Bool isPrime n = millerRabinTest n (factorizeN (n-1)) {- factorizeN finds a number s and odd number d such that n -1 = (2^s)d by succesively dividing n by two if it is even. -} factorizeN :: Int -> (Int, Int) factorizeN n = fN n 0 where fN n s | even n = fN (n `div` 2) (s + 1) | otherwise = (n,s) {- this is the main function. it takes w values from a set of witnesses and checks if n passes the test. If it doesn't, n is not prime, if it does for all w, it is probably prime. -} millerRabinTest :: Int -> (Int,Int) -> Bool millerRabinTest n (d,s) = and [test n (expmod w d n) s | w <- onesToCheck] {- this is the test that is used in the millerRabinTest function. it sees if w^d = 1 mod n or n-1 mod n, if not it multiplies by two and checks again for a total of s-1 times. If it is never true then the number is not prime -} test :: Int -> Int -> Int -> Bool test n w s | w `elem` [1,n-1] = True | otherwise = or [ (expmod w (2^k) n) `elem` [1,n-1]| k <- [1..s]] {- set of witnesses that should make the Miller Rabin test deterministic if n < 2^64. -} onesToCheck :: [Int] onesToCheck = [2,325,9375,28178,450775,9780504,1795265022] {- function that calculates a^e mod n. -} expmod :: Int -> Int -> Int -> Int expmod a e n | e == 1 = a `mod` n | (e `mod` 2) == 0 = (expmod ((a*a) `mod` n) (e `div` 2) n) | otherwise = (a*(expmod ((a*a) `mod` n) (e `div` 2) n)) `mod` n