https://blog.derbuihan.xyz/posts/haskell_primes/

fermat_test.hs

import Data.Numbers.Primes hiding (primes')
import Control.Parallel.Strategies

gcd' :: Integer -> Integer -> Integer
gcd' n 0 = n
gcd' n m = gcd' m (n `mod` m)

exp_mod :: Integer -> Integer -> Integer
exp_mod a n = iter a (n-1) n
  where
    iter a n base
      | n == 1 = a
      | even n = (iter a (n `div` 2) base) ^ 2 `mod` base
      | otherwise = a * (iter a ((n-1) `div` 2) base) ^ 2 `mod` base

primes' :: [Integer]
primes' = 2: filter is_prime' [3,5..]

is_prime' :: Integer -> Bool
is_prime' n = and $ map (\a -> test n a) [2..(n-1)]
  where
    test n a
      | gcd' n a == 1 = (1 == exp_mod a n)
      | otherwise = False

is_prime'' :: Integer -> Bool
is_prime'' n = and $ map (\a -> test n a) $ takeWhile (<n) [2,3,5]
  where
    test n a
      | gcd' n a == 1 = (1 == exp_mod a n)
      | otherwise = False

is_carmichael :: Integer -> Bool
is_carmichael n = and $ ((not . isPrime) n):(map (\a -> test n a) [2..(n-1)])
  where
    test n a
      | gcd' n a == 1 = (1 == exp_mod a n)
      | otherwise = True

is_carmichael' :: Integer -> Bool
is_carmichael' n = and $ map (\a -> test n a) $ takeWhile (< n) primes
  where
    test n a
      | gcd' n a == 1 = (1 == exp_mod a n)
      | otherwise = True

main :: IO ()
main = do
  let ps = takeWhile (< 5000) primes'
      mersenne = runEval $ parFilter (\n -> is_prime'' (2^n - 1)) ps
      carmichael = runEval $ parFilter is_carmichael' $ filter (not . isPrime) $ takeWhile (< 10 ^ 7) [2..]
  print carmichael

parFilter :: (a -> Bool) -> [a] -> Eval [a]
parFilter f [] = return []
parFilter f (a:as) = do
        b <- rpar (f a)
        bs <- parFilter f as
        return (if b then a:bs else bs)

コンパイル

ghc -O2 fermat_test.hs -threaded -rtsopts

実行

./fermat_test +RTS -N4