MonteCarlo simulation: calculate PI

gistfile1.hs

module Main where

import System.Environment
import System.Exit
import System.Random
import Control.Monad
import Control.Concurrent
import Control.Concurrent.MVar

workerFn :: Int -> Int -> MVar Double -> IO ()
workerFn numSamples id m = do
  g <- newStdGen
  samples <- genSamples g numSamples
  let cInside = filter (\(x, y) -> (x * x) + (y * y) <= 1.0) samples
  let nPoints      = fromIntegral (length samples) :: Double
  let pointsInside = fromIntegral (length cInside) :: Double
  let area = (4 * pointsInside / nPoints)
  area `seq` putMVar m area -- force the thread to do the calculations

runWorker :: Int -> Int -> IO (MVar Double)
runWorker numSamples id = do
  m <- newEmptyMVar
  forkIO $ workerFn numSamples id m
  return m

genSamples :: StdGen -> Int -> IO [(Double, Double)]
genSamples g numSamples = do
  return $ take numSamples $ transform (randoms g :: [Double])
  where transform (a:b:r) = (a,b) : transform r

calcPi :: Int -> Int -> IO Double
calcPi numSamples numWorkers = do
  let wrkNumSamples = numSamples `div` numWorkers
  mvars <- sequence [runWorker wrkNumSamples i | i <- [1 .. numWorkers]]
  nums  <- sequence $ map (\m ->
    do v <- takeMVar m
       return v) mvars
  return $ (sum nums) / (fromIntegral numWorkers)

showUsage :: IO ExitCode
showUsage = do
  progName <- getProgName
  putStrLn $ "USAGE: ./" ++ progName ++ " <num-samples> <num-workers>"
  exitFailure

main :: IO ExitCode
main = do
  args <- getArgs
  if length args /= 2
    then showUsage
    else do
      let samples = read (args !! 0) :: Int
      let numWorkers = read (args !! 1) :: Int
      if samples `mod` numWorkers /= 0
        then do
          putStrLn $ "Error: number of samples must be divisible "
                  ++ "by the number of workers"
          showUsage
        else do
          pi <- calcPi samples numWorkers
          putStrLn $ "PI ~= " ++ (show pi)
          exitSuccess