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

Jonathan Cast jonathanccast at fastmail.fm
Mon Nov 5 15:30:00 EST 2007

```On Mon, 2007-11-05 at 20:11 +0000, Alex Young wrote:
> 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
>
> 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?

> main = do

Get two standard generators (one per dimension)

>   g0 <- newStdGen
>   g1 <- newStdGen

Get an infinite list of pairs

>   let pairs = [ (x, y) | x <- randoms (-1, 1) g0,
>                          y <- randoms (-1, 1) g1 ]

Get a finite list

>       consideredPairs = take total pairs

Count how many are in the unit circle

>       circleArea = length \$ filter (\ (x, y) -> x^2 + y^2 < 1)
>                      consideredPairs

Divide by total
>       ratio = fromInteger circleArea / fromIntegral total

Now, pi is approximated by ratio.

>   putStr \$ show ratio

jcc

```