-
Notifications
You must be signed in to change notification settings - Fork 0
/
primeTest.hs
62 lines (51 loc) · 2.17 KB
/
primeTest.hs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
module Primes where
import System.Random
import System.IO.Unsafe
import Control.Monad.Par
import Data.List.Split
-- You need cabal install monad-par to use this!
-- Nice wrapper functions
-- Takes list, returns list of primes in list
-- list list list list list list
filterPrimes :: [Integer] -> [Integer]
filterPrimes a = filter isPrime a
-- Essentially filterPrimes, but in parallel.
-- Uses n threads to do filtering
findPrimes :: [Integer] -> Int -> [[Integer]]
findPrimes a n = runPar $ parMap filterPrimes $ chunksOf (length a `div` n) a
-- The following primality test was coppied shamelessly from RossettaCode --
-- http://rosettacode.org/wiki/Miller-Rabin_primality_test#Haskell --
-- Miller-Rabin wrapped up as an (almost deterministic) pure function
-- We use 40 iterations as it's good enough for our purposes
isPrime :: Integer -> Bool
isPrime n = unsafePerformIO (isMillerRabinPrime 40 n)
isMillerRabinPrime :: Int -> Integer -> IO Bool
isMillerRabinPrime k n
| even n = return (n==2)
| n < 100 = return (n `elem` primesTo100)
| otherwise = do ws <- witnesses k n
return $ and [test n (pred n) evens (head odds) a | a <- ws]
where
(evens,odds) = span even (iterate (`div` 2) (pred n))
test :: Integral nat => nat -> nat -> [nat] -> nat -> nat -> Bool
test n n_1 evens d a = x `elem` [1,n_1] || n_1 `elem` powers
where
x = powerMod n a d
powers = map (powerMod n a) evens
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 m x n = x^n `mod` m
powerMod :: Integral nat => nat -> nat -> nat -> nat
powerMod m x n = f (n - 1) x x `rem` m
where
f d a y = if d==0 then y else g d a y
g i b y | even i = g (i `quot` 2) (b*b `rem` m) y
| otherwise = f (i-1) b (b*y `rem` m)