|
*~> dropWhile (<900) $ filter isPrime [2..1000]
[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}}==
|