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

Alex Young alex at blackkettle.org
Mon Nov 5 15:11:10 EST 2007


Hi all,

I'm new to Haskell, and don't quite understand how IO and lazy 
evaluation work together yet.  In order to improve my understanding, I 
thought I'd try to knock together a translation of a fairly simple 
problem I can code up in C, that of calculating pi by throwing random 
numbers about.  The C code I started with is as follows:

//////////////////////////////////////////////////
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int main(int argc, char *argv[])
{
   if(argc != 2)
   {
     printf("usage: %s iteration_count", argv[0]);
     return 1;
   }

   int total;
   if(sscanf(argv[1], "%d", &total) != 1)
     return 1;

   int count = 0;
   int unit_radius = RAND_MAX * RAND_MAX;
   int i;

   srand(time(NULL));

   for(i = 0; i < total; i++)
   {
     int x, y;
     x = rand();
     y = rand();
     count += x*x + y*y < unit_radius;
   }

   double pi = (count << 2) / (double)total;

   printf("Count: %d\nTotal: %d\n", count, total);
   printf("%f", pi);

   return 0;
}
///////////////////////////////////////////////////

All very simple - the key is that I'm counting the result of a 
conditional evaluated inside a loop.  Ignore whether it's correct, or 
accurate, for now; it happily runs and gives an accurate-looking result 
when run with an argument of 10000000.  My (thoroughly ugly, not 
currently working) Haskell version looks like this:

{--------------------------------------------------}
module Main where

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.

More information, if it's relevant:
   Haskell: Glasgow Haskell Compiler, Version 6.8.0.20070917, for 
Haskell 98, stage 3 booted by GHC version 6.6
   OS: Windows Vista, 32bit

Any guidance gratefully appreciated...

Thanks,
-- 
Alex


More information about the Haskell-Cafe mailing list