parallel PI

{-# LANGUAGE BangPatterns #-}
import System.Environment
import System.Random
import Data.Word
import System.IO
import Data.Bits
import Text.Printf
import Debug.Trace
import Control.Monad
import Data.List
import Control.Concurrent

upBound :: Word32
upBound = 0xffff

type IPair = (Integer, Integer)
type WPair = (Word64, Word64)

accum :: WPair -> IPair -> IPair
accum (a,b) (p, all) = r
  where
    a2 = a `shiftR` 33
    b2 = b `shiftR` 33
    s2 =
      --traceShow ("a2:",printf "%x" a2 :: String,"b2:",printf "%x" b2 :: String) $
      a2 * a2 + b2 * b2
    p' = if s2 >= 0x4000000000000000 then p else p + 1
    all' = all + 1
    r =
      --traceShow ("s2:",printf "%016x" s2 :: String) $
      p' `seq` all' `seq` (p', all')

ncount :: Word64 -> IO IPair
ncount n = do
  (p,all) <- randIter n (0,0)
  return (p * 4, all)
  where
    randIter :: Word64 -> IPair -> IO IPair
    randIter 0 ip = return ip
    randIter n ip = do
      ra <- randomIO
      rb <- randomIO
      let !ip' = accum (ra, rb) ip
      randIter (n - 1) ip'

runPI :: Word64 -> IO Double
runPI n = do
  (p,all) <- ncount n
  return $ fromIntegral p / fromIntegral all


runNPI :: Int -> Word64 -> IO ()
runNPI n pn = do
  mList <- sequence $ take n $ repeat newEmptyMVar
  mapM_ (\m -> forkIO $ runPI pn >>= putMVar m) mList
  mapM_ (\m -> takeMVar m >>= printf "%f") mList

main :: IO ()
main = do
  (a:b:_) <- getArgs
  let n = read a :: Int
  let pn = read b :: Word64
  runNPI n pn

-- using 2cpus * 4cores * 2threads = 16 virtual cores

-- this: all 16 cpu consumes about 30%
-- +RTS -N16 --RTS 1 100000

-- +RTS -N16 --RTS 16 100000
-- +RTS -N2 --RTS 1 100000
-- +RTS -N2 --RTS 2 100000