I implemented a prime number finder using the Miller-Rabin prime test in Haskell and it seems to be working. You enter a number and it finds the next prime following that number.

I also tried to make it parallel and I'm looking for some advice for how I can speed up the parallel version. Right now it runs much slower than the not parallel version.

Compiled with ghc -threaded -O2 miller_rabin.hs -o miller_rabin

Run with ./miller_rabin startnumber npar +RTS -N4 (I have a 4-core processor). npar is the parameter for parBuffer. If npar is 1 then the program runs single-threaded, otherwise is runs with multiple threads.

module Main where

import System.Random (StdGen, getStdGen, randomRs)
import System.Environment (getArgs)
import Control.Parallel.Strategies

main :: IO ()
main = do
    g <- getStdGen
    args <- getArgs 
    let number = read . head $ args :: Integer
        npar = read $ args !! 1 :: Int
    p <- case npar of
            1 -> do putStrLn "Run single thread"
                    return $ find_prime g 20 number
            _ -> do putStrLn "Run multi-threaded"
                    return $ head $ primeList npar g 20 number
    putStrLn . show $ p

miller_rabin :: StdGen -> Int -> Integer -> Bool
miller_rabin g k n = all (not . trial_composite n d r) 
                   $ take k (randomRs (2, n-2) g)    
        (d, r) = decompose (n-1) 0

trial_composite :: Integer -> Integer -> Integer -> Integer -> Bool
trial_composite n d r a = let x = fastPow a d n in
    if (x == 1) || (x==n-1) 
    then False
    else all ((/=) (n-1)) $ map (\i -> fastPow a (d*(2^i)) n) [0..r-1]

decompose :: Integer -> Integer -> (Integer, Integer)
decompose d r 
    | d `mod` 2 == 0 = decompose (d `div` 2) (r+1)
    | otherwise = (d, r) 

fastPow :: Integer -> Integer -> Integer -> Integer
fastPow base 1 m = mod base m
fastPow base pow m | even pow = mod ((fastPow base (div pow 2) m) ^ 2) m
                   | odd  pow = mod ((fastPow base (div (pow-1) 2) m) ^ 2 * base) m

primeList :: Int -> StdGen -> Int -> Integer -> [Integer]
primeList npar g k n = 
    let primes = filter (miller_rabin g k) [m,m+2..]
        primesP = primes `using` parBuffer npar rdeepseq
        m = if even n then n+1 else n
    in primesP

find_prime :: StdGen -> Int -> Integer -> Integer
find_prime g k n 
    | even n    = find_prime g k (n+1)
    | otherwise = case miller_rabin g k n of
                    True  -> n
                    False -> find_prime g k (n+2)
  • \$\begingroup\$ It seems to me like the parallel version is slower because it computes the list elements in parallel, rather than only needing to compute the first list element like the serial version. Make find_prime and primeList calculate and output the same and find_prime shouldn't complete faster anymore. \$\endgroup\$
    – Gurkenglas
    Jun 18, 2019 at 18:44
  • \$\begingroup\$ @Gurkenglas, I’ll try to do a rewrite, but I thought that the parallel list version would allow the list to be computed quickly and there would only be a penalty at the final chunk. \$\endgroup\$
    – user668074
    Jun 18, 2019 at 21:00


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.