[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
> 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.
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
import Control.Monad
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.
More information about the Haskell-Cafe
mailing list