Miller–Rabin primality test: Difference between revisions
Content added Content deleted
Line 1,303: | Line 1,303: | ||
*~> 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> |
||
== Perhaps a slightly clearer (less monadic) version. Transcription of pseudocode. == |
|||
* The code above likely has better complexity. |
|||
<lang Haskell> |
|||
import Control.Monad (liftM) |
|||
import Data.Bits (Bits, testBit, shiftR) |
|||
import System.Random (Random, getStdGen, randomRs) |
|||
import System.IO.Unsafe (unsafePerformIO) |
|||
import Prelude hiding (even, odd) |
|||
odd :: (Integral a, Bits a) => a -> Bool |
|||
odd = (`testBit` 0) |
|||
even :: (Integral a, Bits a) => a -> Bool |
|||
even = not . odd |
|||
-- modPow - Recursive modular exponentiation by taking successive powers of two |
|||
modPow :: (Integral a, Bits a) => a -> a -> a -> a |
|||
modPow _ 0 _ = 1 |
|||
modPow base ex m = let term |
|||
| testBit ex 0 = base `mod` m |
|||
| otherwise = 1 |
|||
in (term * modPow (base^2 `mod` m) (ex `shiftR` 1) m) `mod` m |
|||
isPrime :: (Integral a, Bits a, Random a) => a -> a -> Bool |
|||
isPrime n k |
|||
| n < 4 = if n > 1 then True else False -- Deal with 0-3. |
|||
| even n = False |
|||
| otherwise = let randPool = unsafePerformIO $ randNums (n - 2) |
|||
in witness k randPool |
|||
where |
|||
randNums upper = do |
|||
g <- getStdGen |
|||
return (randomRs (2, upper) g) |
|||
(d, r) = let decompose d r |
|||
| odd d = (d, r) |
|||
| otherwise = decompose (d `shiftR` 1) (r + 1) |
|||
in decompose (n - 1) 0 |
|||
witness 0 _ = True |
|||
witness k (a:rands) |
|||
| x == 1 || x == n - 1 = witness (k - 1) rands |
|||
| otherwise = check x (r - 1) |
|||
where |
|||
x = modPow a d n |
|||
check _ 0 = False |
|||
check x count |
|||
| x' == 1 = False |
|||
| x' == n - 1 = witness (k - 1) rands |
|||
| otherwise = check x' (count - 1) |
|||
where x' = modPow x 2 n |
|||
-- main function for testing |
|||
main :: IO() |
|||
main = do |
|||
[n,k] <- liftM (map (\x -> read x :: Integer) . words) getLine |
|||
print $ isPrime n k |
|||
</lang> |
|||
{{out|Sample Output}} |
|||
<pre> |
|||
4547337172376300111955330758342147474062293202868155909489 10 |
|||
True |
|||
4547337172376300111955330758342147474062293202868155909393 10 |
|||
False |
|||
226801 7 |
|||
False |
|||
</pre> |
|||
=={{header|Icon}} and {{header|Unicon}}== |
=={{header|Icon}} and {{header|Unicon}}== |