Rosetta Code

Miller–Rabin primality test: Difference between revisions

Page Discussion

Miller–Rabin primality test (view source)

Revision as of 23:02, 5 October 2013

225 bytes added ,  10 years ago
→‎{{header|Haskell}}
(No function doing the task)
Line 1,170:
 
=={{header|Haskell}}==
{{works with|Haskell|67.106.43}}
* 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>importmodule System.RandomPrimes 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
isMillerRabinPrimeisPrime :: Integer -> IO Bool▼
isPrime n = unsafePerformIO (isMillerRabinPrime n)
 
isMillerRabinPrime n:: |Integer even-> nIO = return (2 == n)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
e = uncurry (++evens,odds) . second(take 1) .= span even . (iterate (`div` 2) $ pnn_1)▼
powers = map (powerMod n w) (evens ++ take 1 odds)
 
witnswitnesses :: (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 ▼
then return $ take x $k (randomRs (2,yn-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 ::m (Integralx a, Integral b)n => ax^n -> a -> b ->`mod` am
powerMod m:: _Integral 0nat => nat 1-> nat -> nat -> nat
powerMod m x n | n > 0 = join (flip f (n - 1)) x x `rem` m where
where
f _ 0 y = y
f d a y = if d==0 then y =else g d a dy where
g i b iy | even i = g (b*bi `remquot` m2) (ib*b `quotrem` 2m) y
| otherwise = f b (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 ⟶ 1,226:
False
 
*~> filterMdropWhile isMillerRabinPrime(<900) $ filter isPrime [2..1000]>>= print. dropWhile(<=900)
[907,911,919,929,937,941,947,953,967,971,977,983,991,997]</pre>
 
Anonymous user
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/107833"