Haskell prime number generator according to bits for very large numbers

Question

I just got into haskell. I'm trying to write a prime number generator that would return me a prime number depending on the number of bits I give it. Do I use some sort of Sieve of Eratosthenes? Would this be the quickest way? Currently, I already have a prime number checker for Miller-Rabin. Is there a correct way and a wrong way to do this? Also, I want to be able to generate large numbers very quickly.

Ex. generate a prime number that is of 32 bits

genp 32

Code so far.

import System.IO
import System.Random
import Data.List
import Control.Monad
import Control.Arrow


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 :: (Integral a, Integral b) => a -> a -> b -> a
powerMod m _ 0 =  1
powerMod m x n | n > 0 = join (flip f (n - 1)) x `rem` m where
  f _ 0 y = y
  f a d y = g a d where
    g b i | even i    = g (b*b `rem` m) (i `quot` 2)
          | otherwise = f b (i-1) (b*y `rem` m)

witns :: (Num a, Ord a, Random a) => Int -> a -> IO [a]
witns x y = do
     g <- newStdGen 
     let r = [9080191, 4759123141, 2152302898747, 3474749600383, 341550071728321]
         fs = [[31,73],[2,7,61],[2,3,5,7,11],[2,3,5,7,11,13],[2,3,5,7,11,13,17]]
     if  y >= 341550071728321
      then return $ take x $ randomRs (2,y-1) g
       else return $ snd.head.dropWhile ((<= y).fst) $ zip r fs

isMillerRabinPrime :: Integer -> IO Bool
isMillerRabinPrime n | n `elem` primesTo100 = return True
                     | otherwise = do
    let pn = pred n
        e = uncurry (++) . second(take 1) . span even . iterate (`div` 2) $ pn
        try = return . all (\a -> let c = map (powerMod n a) e in
                                  pn `elem` c || last c == 1)
    witns 100 n >>= try



type Prime = Integer

isProbablyPrime :: Prime -> Bool
isProbablyPrime n = isMillerRabinPrime n

pickFirstFrom :: Integer -> Prime
pickFirstFrom n = head $ filter isProbablyPrime [n..]

numBits = 1024
constantStdGen = mkStdGen 1234567 -- a random number

randomByBits n = fst $ randomR (2^(n-1), 2^n-1) constantStdGen

answer = pickFirstFrom (randomByBits numBits)

Using Library functions with RSA:

import Control.Monad.Fix 
import Math.NumberTheory.Primes

rndPrime :: Int -> IO Integer
rndPrime bits =
    fix $ \again -> do
        x <- fmap (.|. 1) $ randomRIO (2^(bits - 1), 2^bits - 1)
        if isPrime x then return x else again

rndPrimes :: Int -> IO (Integer, Integer)
rndPrimes bits = do
    p <- rndPrime bits
    fix $ \again -> do
        q <- rndPrime bits
        if p /= q then return (p, q) else again

Thanks guys, I really appreciate your help.

Answer 1

For generating an odd random prime number of a specific bit size there is this de facto way of doing it, assuming a (usually probabilistic) primality test isPrime:

rndPrime :: Int -> IO Integer
rndPrime bits =
    fix $ \again -> do
        x <- fmap (.|. 1) $ randomRIO (2^(bits - 1), 2^bits - 1)
        if isPrime x then return x else again

For RSA:

rndPrimes :: Int -> IO (Integer, Integer)
rndPrimes bits = do
    p <- rndPrime bits
    fix $ \again -> do
        q <- rndPrime bits
        if p /= q then return (p, q) else again

You can find ready-made fast primality tests in the arithmoi package, including the current state of the art, the Baillie PSW test. Then just import the Math.NumberTheory.Primes module and the code above should work out of the box.

Answer 2

You're in luck mister. The prime numbers are appearing so frequently that you just can try to enumerate all numbers from a given start number and pick the first one that is prime.

The solution would be something like this

import System.Random

type Prime = Integer

isProbablyPrime :: Prime -> Bool
isProbablyPrime n = error "insert miller rabin here"

pickFirstFrom :: Integer -> Prime
pickFirstFrom n = head $ filter isProbablyPrime [n..]

numBits = 1024
constantStdGen = mkStdGen 1234567 -- a random number

randomByBits n = fst $ randomR (2^(n-1), 2^n-1) constantStdGen

answer = pickFirstFrom (randomByBits numBits)

Here answer will be as you specified in your comments. Note however that you have to modify this code to allow more than one constant seed!