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. {-# LANGUAGE BangPatterns #-} import System( getArgs ) import Data.List( foldl' ) import System.Random.Mersenne pairs :: [a] -> [(a,a)] pairs [] = [] pairs (x:[]) = [] pairs (x:y:rest) = (x, y) : pairs rest isInCircle :: (Double, Double) -> Bool isInCircle (x,y) = sqrt (x*x + y*y) <= 1.0 accumulateHit :: Int -> (Double, Double) -> Int accumulateHit (!hits) pair | isInCircle pair = hits + 1 | otherwise = hits countHits :: [(Double, Double)] -> Int countHits ps = foldl' accumulateHit 0 ps monteCarloPi :: Int -> [(Double, Double)] -> Double monteCarloPi n xs = 4.0 * fromIntegral hits / fromIntegral n where hits = countHits $ take n xs main = do args <- getArgs let samples = read $ head args randomNumberGenerator <- getStdGen randomNumbers <- randoms randomNumberGenerator let res = monteCarloPi samples $ pairs randomNumbers putStrLn $ show $ res