Generate large primes

README.md

Compile with:

ghc -O2 LargePrime.hs -main-is LargePrime -o large-prime -rtsopts -threaded

The only argument is the number of bits for the desired prime, so

./large-prime 4096 +RTS -N4

gives a 4096-bit prime and tries to use 4 cores.

LargePrime.hs

module LargePrime where

import Control.Parallel ( par, pseq )
import GHC.Conc ( numCapabilities )
import System.Environment ( getArgs )
import System.Random ( newStdGen, randomRs )

primes :: [Integer]
primes = 2 : filter (go primes) [3, 5..]
    where
      go [] _ = True
      go (p:ps) n
          | p*p > n        = True
          | n `mod` p == 0 = False
          | otherwise      = go ps n

pmap :: [a] -> [a]
pmap vs = foldl1 par vs `pseq` vs

fermatTest :: [Integer] -> Integer -> Bool
fermatTest [] _ = True
fermatTest as n = if all id (pmap $ map divTest cs)
                  then fermatTest as' n
                  else False
    where
      (cs,as') = splitAt numCapabilities as
      divTest a = modPow a (n-1) n == 1

modPow :: Integer -> Integer -> Integer -> Integer
modPow _ 0 _ = 1
modPow x y n
    | y `mod` 2 == 0 = (modPow x' (y `div` 2) n)^(2 :: Int) `mod` n
    | otherwise      = x' * modPow x' (y-1) n
    where
      x' = x `mod` n

randomNumbers :: Integer -> Integer -> IO [Integer]
randomNumbers lo hi = do
  gen <- newStdGen
  return $ randomRs (lo, hi) gen

largeRands :: Int -> IO [Integer]
largeRands sz = randomNumbers (2^sz) (2^(sz+1))

isPrime :: Integer -> IO Bool
isPrime n = do
  as <- randomNumbers 1 n
  let as' = take 100 primes ++ take 100 as
  return $ fermatTest as' n

findLargePrime :: Int -> IO Integer
findLargePrime sz = go =<< largeRands sz
    where
      go [] = fail "ran out of large numbers; weird"
      go (r:rs) = do
        ip <- isPrime r
        if ip then return r else go rs

main :: IO ()
main = do
  [sz] <- map read `fmap` getArgs
  print =<< findLargePrime sz