[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