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

Don Stewart dons at galois.com
Mon Nov 5 15:40:21 EST 2007

```alex:
> Hi all,
>
> import Random
> import System.Environment
> import List
>
> 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
> 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

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.
```