[Haskell-beginners] parallelizing a function for generating prime numbers

[mukesh tiwari]

You can use Rabin-Miller[1] primality testing.  The idea is, divide the
range  1 to 2000000 in chunk of 10000 numbers and evaluate the all the
chunks in parallel.

import Data.Bits
import Control.Parallel.Strategies

powM :: Integer -> Integer -> Integer -> Integer
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 powM a d n of
          1 -> True
          x -> any ( == pred n ) . take ( fromIntegral s )
                      . iterate (\e -> mod ( e^2 ) n ) $ x


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



sum' :: Integer -> Integer -> Integer
sum' m n = sum . map ( \( x, y ) -> if y then x else 0 ) .  zip [ m .. n ]
. primeRange m  $ n


main = print ( sum' 1  2000000 )


Mukeshs-MacBook-Pro:Puzzles mukeshtiwari$ time ./Proj10 +RTS -N2
142913828922

real    0m6.301s
user    0m11.937s
sys    0m0.609s
Mukeshs-MacBook-Pro:Puzzles mukeshtiwari$ time ./Proj10 +RTS -N1
142913828922

real    0m8.202s
user    0m8.026s
sys    0m0.174s


-Mukesh
[1] http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test


On Fri, May 16, 2014 at 9:23 PM, Norbert Melzer <timmelzer at gmail.com> wrote:

> Hi there!
>
> I am trying to enhence the speed of my Project Euler solutions…
>
> My original function is this:
>
> ```haskell
> problem10' ::  Integer
> problem10' = sum $ takeWhile (<=2000000) primes
>   where
>     primes                  = filter isPrime possiblePrimes
>     isPrime n               = n == head (primeFactors n)
>     possiblePrimes          = (2:3:concat [ [6*pp-1, 6*pp+1] | pp <- [1..]
> ])
>     primeFactors m          = pf 2 m
>     pf n m | n*n > m        = [m]
>            | n*n       == m  = [n,n]
>            | m `mod` n == 0  = n:pf n (m `div` n)
>            | otherwise      = pf (n+1) m
> ```
>
> Even if the generation of primes is relatively slow and could be much
> better, I want to focus on parallelization, so I tried the following:
>
> ```haskell
> parFilter :: (a -> Bool) -> [a] -> [a]
> parFilter _ [] = []
> parFilter f (x:xs) =
>   let px = f x
>       pxs = parFilter f xs
>   in par px $ par pxs $ case px of True -> x : pxs
>                                    False -> pxs
>
> problem10' ::  Integer
> problem10' = sum $ takeWhile (<=2000000) primes
>   where
>     primes                  = parFilter isPrime possiblePrimes
>     isPrime n               = n == head (primeFactors n)
>     possiblePrimes          = (2:3:concat [ [6*pp-1, 6*pp+1] | pp <- [1..]
> ])
>     primeFactors m          = pf 2 m
>     pf n m | n*n > m        = [m]
>            | n*n       == m  = [n,n]
>            | m `mod` n == 0  = n:pf n (m `div` n)
>            | otherwise      = pf (n+1) m
> ```
>
> This approach was about half as slow as the first solution (~15 seconds
> old, ~30 the new one!).
>
> Trying to use `Control.Parallel.Strategies.evalList` for `possiblePrimes`
> resulted in a huge waste of memory, since it forced to generate an endless
> list, and does not stop…
>
> Trying the same for `primeFactors` did not gain any speed, but was not
> much slower at least, but I did not expect much, since I look at its head
> only…
>
> Only thing I could imagine to parallelize any further would be the
> takeWhile, but then I don't get how I should do it…
>
> Any ideas how to do it?
>
> TIA
> Norbert
>
>
> _______________________________________________
> Beginners mailing list
> Beginners at haskell.org
> http://www.haskell.org/mailman/listinfo/beginners
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.haskell.org/pipermail/beginners/attachments/20140518/3d5b5d2b/attachment.html>