Miller–Rabin primality test: Difference between revisions

Line 1,170:

=={{header|Haskell}}==

* Ideas taken from [http://primes.utm.edu/prove/prove2_3.html Primality proving]

* Functions witns and isMillerRabinPrime follow closely the code outlined in [http://www.jsoftware.com/jwiki/Essays/Primality%20Tests#Miller-Rabin J/Essays]

* A useful powerMod function is taken from [[Multiplicative order#Haskell]]

* Original Rosetta code has been simplified to be easier to follow

Another Miller Rabin test can be found in D. Amos's Haskell for Math module [http://www.polyomino.f2s.com/david/haskell/numbertheory.html Primes.hs]

<lang Haskell>~~import~~ ~~System.Random~~

<lang Haskell>module Primes where

import Data.List

import Control.Monad

import Control.Arrow

import System.Random

import System.IO.Unsafe

-- Miller-Rabin wrapped up as an (almost deterministic) pure function

⚫

isPrime :: Integer -> Bool

isPrime n = unsafePerformIO (isMillerRabinPrime n)

⚫

isMillerRabinPrime :: Integer -> IO Bool

isMillerRabinPrime n

| even n = return (n==2)

| n < 100 = return (n `elem` primesTo100)

⚫

| otherwise = do ws <- witnesses 100 n

return $ and (map (test n (pred n)) ws)

test :: Integral nat => nat -> nat -> nat -> Bool

test n n_1 w = n_1 `elem` powers || last powers == 1

where

⚫

(evens,odds) = span even (iterate (`div` 2) n_1)

powers = map (powerMod n w) (evens ++ take 1 odds)

⚫

witnesses :: (Num a, Ord a, Random a) => Int -> a -> IO [a]

witnesses k n

| n < 9080191 = return [31,73]

| n < 4759123141 = return [2,7,61]

| n < 3474749660383 = return [2,3,5,7,11,13]

| n < 341550071728321 = return [2,3,5,7,11,13,17]

⚫

| otherwise = do g <- newStdGen

⚫

return $ take k (randomRs (2,n-1) g)

primesTo100 :: [Integer]

primesTo100 = [2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97]

powerMod :: ~~(Integral~~ ~~a, Integral b)~~ => a ~~-> a -> b ->~~ a

-- powerMod m x n = x^n `mod` m

powerMod m _ 0 = 1

powerMod :: Integral nat => nat -> nat -> nat -> nat

powerMod m x n ~~| n > 0~~ = ~~join (flip~~ f (n - 1)) x `rem` m ~~where~~

powerMod m x n = f (n - 1) x x `rem` m

where

f _ 0 y = y

f a d y = g a d ~~where~~

f d a y = if d==0 then y else g d a y

g b i | even i = g (~~b*b~~ `~~rem~~` m) (i `~~quot~~` 2)

g i b y | even i = g (i `quot` 2) (b*b `rem` m) y

| otherwise = f b (i-1) (b*y `rem` m)

| otherwise = f (i-1) b (b*y `rem` m)

⚫

~~witns~~ :: (Num a, Ord a, Random a) => Int -> a -> IO [a]

witns x y = do

⚫

g <- newStdGen

let r = [9080191, 4759123141, 2152302898747, 3474749660383, 341550071728321]

fs = [[31,73],[2,7,61],[2,3,5,7,11],[2,3,5,7,11,13],[2,3,5,7,11,13,17]]

if y >= 341550071728321

⚫

~~then~~ return $ take ~~x $~~ randomRs (2,y-1) g

else return $ snd.head.dropWhile ((<= y).fst) $ zip r fs

⚫

~~isMillerRabinPrime~~ :: Integer -> IO Bool

⚫

isMillerRabinPrime n | ~~even~~ n ~~= return (2 == n)~~

| n `elem` primesTo100 = return True

⚫

| otherwise = do

let pn = pred n

⚫

~~e = uncurry~~ (++) ~~. second(take 1) .~~ span even . iterate (`div` 2) ~~$ pn~~

try = return . all (\a -> let c = map (powerMod n a) e in

pn `elem` c || last c == 1)

witns 100 n >>= try</lang>

Testing in GHCi:

<pre>*~> isMillerRabinPrime 4547337172376300111955330758342147474062293202868155909489

Line 1,215:

Line 1,226:

False

*~> ~~filterM~~ ~~isMillerRabinPrime~~ [2..1000]~~>>= print. dropWhile(<=900)~~

*~> dropWhile (<900) $ filter isPrime [2..1000]

[907,911,919,929,937,941,947,953,967,971,977,983,991,997]</pre>