Miller–Rabin primality test: Difference between revisions

Content added Content deleted
Line 1,183: Line 1,183:
-- Miller-Rabin wrapped up as an (almost deterministic) pure function
-- Miller-Rabin wrapped up as an (almost deterministic) pure function
isPrime :: Integer -> Bool
isPrime :: Integer -> Bool
isPrime n = unsafePerformIO (isMillerRabinPrime n)
isPrime n = unsafePerformIO (isMillerRabinPrime 100 n)



⚫
isMillerRabinPrime :: Integer -> IO Bool
isMillerRabinPrime n
isMillerRabinPrime :: Int -> Integer -> IO Bool
⚫
isMillerRabinPrime k n
| even n = return (n==2)
| even n = return (n==2)
| n < 100 = return (n `elem` primesTo100)
| n < 100 = return (n `elem` primesTo100)
| otherwise = do ws <- witnesses 100 n
| otherwise = do ws <- witnesses k n
return $ and (map (test n (pred n)) ws)
return $ and [test n (pred n) evens d a | a <- ws]
where
⚫
(evens,odds) = span even (iterate (`div` 2) (pred n))
⚫
d = head odds



test :: Integral nat => nat -> nat -> nat -> Bool
test :: Integral nat => nat -> nat -> [nat] -> nat -> nat -> Bool
test n n_1 w = n_1 `elem` powers || last powers == 1
test n n_1 evens d a = x `elem` [1,n_1] || n_1 `elem` powers
where
where
x = powerMod n a d
⚫
(evens,odds) = span even (iterate (`div` 2) n_1)
powers = map (powerMod n w) (evens ++ take 1 odds)
powers = map (powerMod n a) evens


witnesses :: (Num a, Ord a, Random a) => Int -> a -> IO [a]
witnesses :: (Num a, Ord a, Random a) => Int -> a -> IO [a]
Line 1,209: Line 1,214:
primesTo100 :: [Integer]
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]
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 x n = x^n `mod` m
-- powerMod m x n = x^n `mod` m
Line 1,218: Line 1,224:
| otherwise = f (i-1) b (b*y `rem` m)
| otherwise = f (i-1) b (b*y `rem` m)


--------------------------------------


Testing in GHCi:
Testing in GHCi:
Line 1,228: Line 1,235:
*~> dropWhile (<900) $ filter isPrime [2..1000]
*~> dropWhile (<900) $ filter isPrime [2..1000]
[907,911,919,929,937,941,947,953,967,971,977,983,991,997]</pre>
[907,911,919,929,937,941,947,953,967,971,977,983,991,997]</pre>

=={{header|Icon}} and {{header|Unicon}}==
The following code works in both languages:
<lang unicon>procedure main()
every writes(primeTest(901 to 1000, 10)," ")
write()
end

procedure primeTest(n, k)
if n = 2 then return n
if n%2 = 0 then fail
s := 0
⚫
d := n-1
while (d%2 ~= 0, s+:=1, d/:=2)

every (1 to k, x := ((1+?(n-2))^d)%n) do {
if x = (1 | (n-1)) then next
every (1 to s-1, x := (x^2)%n) do {
if x = 1 then fail
if x = n-1 then break next
}
fail
}
return n
end</lang>
{{out|Sample run}}
<pre>->mrpt
907 911 919 929 937 941 947 953 967 971 977 983 991 997
-></pre>


=={{header|J}}==
=={{header|J}}==
Retrieved from "https://rosettacode.org/wiki/Miller–Rabin_primality_test"