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

[Jonathan Cast]

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
> 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?

> 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