[Haskell-cafe] Monte Carlo Pi calculation (newbie learnings)

[Don Stewart]

alex:
> Hi all,
> 
> import Random
> import System.Environment
> import List
> import Monad
> 
> randMax = 32767
> unitRadius = randMax * randMax
> 
> rand :: IO Int
> rand = getStdRandom (randomR (0, randMax))
> 
> randListTail accum 0 = accum
> randListTail accum n = randListTail (rand : accum) (n - 1)
> 
> randList :: Int -> [IO Int]
> randList n = randListTail [] n
> 
> randPairs :: Int -> [(IO Int, IO Int)]
> randPairs n = zip (randList n) (randList n)
> 
> pairIsInside x y = if x*x + y*y < unitRadius then 1 else 0
> 
> doCountPair :: (IO Int, IO Int) -> IO Int
> doCountPair (a, b) = do
>   x <- a
>   y <- b
>   return (pairIsInside x y)
> 
> fSumListTail :: Int -> [(IO Int, IO Int)] -> IO Int
> fSumListTail total [] = do
>   return total
> fSumListTail total (x:xs) = do
>   y <- doCountPair x
>   fSumListTail (total+y) xs
> 
> fSumList :: [(IO Int, IO Int)] -> IO Int
> fSumList l = fSumListTail 0 l
> 
> piAsDouble :: Int -> Int -> Double
> piAsDouble a b =
>   (fromInteger (toInteger a)) / (fromInteger (toInteger b))
> 
> calculatePi total = do
>   count <- fSumList (randPairs total)
>   return (piAsDouble (4*count) total)
> 
> main = do
>   args <- getArgs
>   (calculatePi (read (args !! 0))) >>= (putStr . show)
> {--------------------------------------------------}
> 
> Now, I have two questions.  The easy question is, how can I make this 
> more idiomatic?  I seem to be jumping through hoops rather a lot to 
> achieve what should be rather simple operations.  The piAsDouble and 
> fSumListTail functions are perfect examples, but I'm not wildly keen on 
> doCountPair either.
> 
> The second question is performance-related.  The Haskell code above 
> overflows the stack when run with an argument of 1000000, so either 
> somewhere a list I'm intending to be lazily evaluated and chucked away 
> is being retained, or one of my recursive functions isn't 
> tail-optimising.  Is it obvious to a more trained eye where the problem 
> is?  If not, how should I start to diagnose this?  I'm not really sure 
> where I should look for more information.

You can replace most of your loops with Data.List functions, and
simplify the code overall by threading around a lazy list of randoms,
rather than calling into IO all the time:

    import System.Random
    import System.Environment
    import Data.List
    import Control.Monad

    randMax    = 32767
    unitRadius = randMax * randMax

    countPair :: (Int, Int) -> Int
    countPair (x, y) = fromEnum (x*x + y*y < unitRadius)

    calculatePi total g = fromIntegral (4*count) / fromIntegral total
     where
      count = sum
            . map countPair
            . take total
            $ zip (randomRs (0,randMax) a)
                  (randomRs (0,randMax) b)

      (a,b) = split g

    main = do
      [v] <- getArgs
      g <- newStdGen
      print $ calculatePi (read v) g

Compiled like so:

    $ ghc -O2 A.hs -o A

    $ time ./A 100000
    3.13548
    ./A 100000  0.08s user 0.02s system 98% cpu 0.101 total

We get no stack overflow. But note you're implementing a different
algorithm to the C program, with different data structures, so expect
different performance.