My Weblog

Blog about programming and math

Generating primes in parallel

We can use Miller-Rabin primality testing to test a number is prime or not with certain probability. We can use this method on chunk of data to get the primes in parallel. I wrote this code to see how it goes and it’s quite interesting. My other post about data parallelism and Don Stewart answer regarding different options to improve parallelism. The crux of this code is

primeRange :: Integer -> Integer -> [ Bool ]
primeRange m n = ( map isPrime [ m .. n ] ) `using` parListChunk 10000 rdeepseq 

For parallelism in Haskell , see Haskell-for-multicores and Real World Haskell. Code on github.

import Data.Bits
import Control.Parallel.Strategies 

powM :: Integer -> Integer -> Integer ->  Integer
powM a d n = powM' a d n where 
  powM' a d n 
        | d == 0 = 1
        | d == 1 = mod a n
        | otherwise = mod q n  where
	       p = powM'  ( mod ( a^2 ) n ) ( shiftR d 1 ) n
	       q = if (.&.) d 1 == 1 then mod ( a * p ) n else p

calSd :: Integer -> ( Integer , Integer ) 
calSd n = ( s , d ) where 
      s = until ( \x -> testBit ( n - 1) ( fromIntegral x ) )  ( +1 ) 0
      d = div ( n - 1 ) (  shiftL 1 ( fromIntegral s )  )


rabinMiller::Integer->Integer->Integer->Integer-> Bool
rabinMiller  n s d a 
   | n == a = True 
   | otherwise =  case x == 1 of 
          True -> True 
          _ ->   any ( == pred n ) . take ( fromIntegral s ) 
                      . iterate (\e -> mod ( e^2 ) n ) $ x  
        where 
              x = powM a d n 
   

isPrime::Integer-> Bool
isPrime n 
   | n <= 1 = False
   | n == 2 = True
   | even n = False
   | otherwise	= all ( == True ) . map ( rabinMiller n s d ) $  [ 2 , 3 , 5 , 7 , 11 , 13 , 17 ] where 
                ( s , d ) = calSd n 



primeRange :: Integer -> Integer -> [ Bool ]
primeRange m n = ( map isPrime [ m .. n ] ) `using` parListChunk 10000 rdeepseq 


main = do 
     let t = filter ( == True ) . primeRange 2 $ 10 ^ 6 
     print . length $ t 

Mukeshs-MacBook-Pro:Haskell mukeshtiwari$ ghc -O2 -rtsopts --make RabinMiller.hs -threaded  -fforce-recomp
[1 of 1] Compiling Main             ( RabinMiller.hs, RabinMiller.o )
Linking RabinMiller ...
Mukeshs-MacBook-Pro:Haskell mukeshtiwari$ time ./RabinMiller +RTS -N1
78498

real	0m4.841s
user	0m4.695s
sys	0m0.075s
Mukeshs-MacBook-Pro:Haskell mukeshtiwari$ time ./RabinMiller +RTS -N2
78498

real	0m3.115s
user	0m5.969s
sys	0m0.215s
Mukeshs-MacBook-Pro:Haskell mukeshtiwari$ time ./RabinMiller +RTS -N3
78498

real	0m2.687s
user	0m7.556s
sys	0m0.415s

I am no expert in parallelism or Haskell so if you have any comment or suggestion then please let me know.

Advertisements

February 5, 2013 - Posted by | Haskell, Programming | , , , , , ,

3 Comments »

  1. can u pls give the equivalent c++ or c code…??…so that many would be benefited from this….

    Comment by jayanth | February 5, 2013 | Reply

    • @Jayanth Writing the same code will take lot of time , efforts and debugging 🙂 . I have given the link of data parallelism so you can try take the array the divide it into chunks of 10000 elements and run each chunk in separate thread. This way you can learn threading in C because for this problem your threads are not communicating. Also why not learn Haskell itself to explore more parallelism 🙂

      Comment by tiwari_mukesh | February 5, 2013 | Reply

      • @tiwari:okay…ll try learning haskell quite soon…:D…and thanks for the post….:)

        Comment by jayanth | February 5, 2013


Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: