But you can remove sqrt from the C++ implementation as well, so it only improves the relative performance if the C++ implementation of sqrt is worse than its Haskell counterpart.

Oops, of course I mean, you only improve if Haskell's implementation is worse than C++'s implementation :)

From: haskell-cafe-bounces@haskell.org [mailto:haskell-cafe-bounces@haskell.org] On Behalf Of Roel van Dijk

I replaced the standard random number generated with the one from mersenne-random. On my system this makes the resulting program about 14 times faster than the original. I also made a change to accumulateHit because it doesn't need to count to total. That is already known.

I tried this too, but got a seg fault (!), so I stripped it back to a small test program. This is with mersenne-random, setup configured with -fuse_sse2: module Main (main) where import System.Random.Mersenne import System (getArgs) rands :: Int -> IO [Double] rands samples = do rng <- getStdGen rns <- randoms rng return (take (2*samples) rns) main = do args <- getArgs let samples = read (head args) randomNumbers <- rands samples print randomNumbers On WinXp with ghc-6.10.1 it always seg faults. If I rebuild without -fuse-sse2 then it's happy. I'm assuming that my machine has SSE2; it's a little old, but it's a P4 (3GHz), so I think it should. Is there an easy way to tell? Alistair ***************************************************************** Confidentiality Note: The information contained in this message, and any attachments, may contain confidential and/or privileged material. It is intended solely for the person(s) or entity to which it is addressed. Any review, retransmission, dissemination, or taking of any action in reliance upon this information by persons or entities other than the intended recipient(s) is prohibited. If you received this in error, please contact the sender and delete the material from any computer. *****************************************************************

Might also mean bad alignment of data structures, maybe when marshalling data between the GHC runtime and the C implementation? If I recall correctly SSE2 requires 16-byte alignment. On Thu, Feb 26, 2009 at 5:24 PM, Don Stewart wrote:

Alistair.Bayley:

...

...

From: haskell-cafe-bounces@haskell.org [mailto:haskell-cafe-bounces@haskell.org] On Behalf Of Roel van Dijk

I replaced the standard random number generated with the one from mersenne-random. On my system this makes the resulting program about 14 times faster than the original. I also made a change to accumulateHit because it doesn't need to count to total. That is already known.

I tried this too, but got a seg fault (!), so I stripped it back to a small test program. This is with mersenne-random, setup configured with -fuse_sse2:

This in the past has always meant: wrong architecture (or GCC can't handle sse2 on your system)

-- Don _______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe

On Thu, Feb 26, 2009 at 6:23 PM, Don Stewart wrote:

But note the lazy list of Double pairs, so the inner loop still looks like this though:

   $wlgo :: Int# -> [(Double, Double)] -> Int

   $wlgo =      \ (ww_s1pv :: Int#)        (w_s1px :: [(Double, Double)]) ->        case w_s1px of wild_aTl {          [] -> I# ww_s1pv;          : x_aTp xs_aTq ->            case x_aTp of wild1_B1 { (x1_ak3, y_ak5) ->            case x1_ak3 of wild2_aX8 { D# x2_aXa ->            case y_ak5 of wild3_XYs { D# x3_XYx ->            case <=##                   (sqrtDouble#                      (+##                         (*## x2_aXa x2_aXa) (*## x3_XYx x3_XYx)))                   1.0            of wild4_X1D {              False -> $wlgo ww_s1pv xs_aTq;              True -> $wlgo (+# ww_s1pv 1) xs_aTq            }

while we want to keep everything in registers with something like:

   Int# -> Double# -> Double# -> Int#

So we'll be paying a penalty to force the next elem of the list (instead of just calling the Double generator).  This definitely has an impact on performance.

   $ ghc-core B.hs -O2 -fvia-C -optc-O3 -fexcess-precision -optc-march=core2 -funbox-strict-fields

   $ time ./B 10000000    3.1407688    ./B 10000000  2.41s user 0.01s system 99% cpu 2.415 total

Now, what if we just rewrote that inner loop directly to avoid intermediate stuff? That'd give us a decent lower bound.

   {-# LANGUAGE BangPatterns #-}

   import System.Environment    import System.Random.Mersenne

   isInCircle :: Double -> Double -> Bool    isInCircle x y = sqrt (x*x + y*y) <= 1.0

   countHits :: Int -> IO Int    countHits lim = do        g <- newMTGen Nothing        let go :: Int -> Int -> IO Int            go !throws !hits                | throws >= lim  = return hits                | otherwise = do                    x <- random g   -- use mersenne-random-pure64 to stay pure!                    y <- random g                    if isInCircle x y                        then go (throws+1) (hits+1)                        else go (throws+1) hits        go 0 0

   monteCarloPi :: Int -> IO Double    monteCarloPi n = do        hits <- countHits n        return $ 4.0 * fromIntegral hits / fromIntegral n

   main = do        [n] <- getArgs        res <- monteCarloPi (read n)        print res

And now the inner loop looks like:

     $wa_s1yW :: Int#                  -> Int#                  -> State# RealWorld                  -> (# State# RealWorld, Int #)

Pretty good. Can't avoid the Int boxed return (and resulting heap check) due to use of IO monad. But at least does away with heap allocs in the inner loop!

How does it go:

   $ ghc-core A.hs -O2 -fvia-C -optc-O3 -fexcess-precision -optc-march=core2 -funbox-strict-fields

   $ time ./A 10000000    3.1412564    ./A 10000000  0.81s user 0.00s system 99% cpu 0.818 total

Ok. So 3 times faster. Now the goal is to recover the high level version. We have many tools to employ: switching to mersenne-random-pure64 might help here. And seeing if you can fuse filling a uvector with randoms, and folding over it... t

-- Don

Very nice! I also wrote a naive version which used uvector but it was about twice as slow as the original Haskell version. I wanted to write "lengthU . filterU isInCircle" because that clearly expresses the algorithm. Sadly I was at work and didn't have time for profile the program to see what was wrong. Still, I couldn't resist having a go at the problem :-)