The fastest way to generate a billion random doubles in Haskell

Question

I am doing simulation in Monte Carlo and am currently using System.Random .

 import System.Random main = do g <- newStdGen let xs = randoms g :: [Double] -- normally, I'd do other magic here putStrLn $ show $ length $ take 10^9 xs 

Unfortunately, this takes a very long time, at least 5 times slower than Python random.random() , not to mention calling C rand() .

With ghc -O2 -optc-ffast-math -optc-O3 -O3

 import System.Random main = do g <- newStdGen let xs = randoms h :: [Double] putStrLn $ show $ length $ take (10^7) xs 

takes ~ 8s against (in iPython)

 import random %timeit len([random.random() for _ in range(10 ** 7)]) 

takes ~ 1.3s. My goal is one billion, but Haskell cannot generate them for a reasonable amount of time.

I also have a C ++ program that generates a float with rand() . He makes 10^7 samples in 0.2s.

How can I quickly generate random doubles in the range [0-1) in Haskell?

Ideally, the GHC program only generates platoons of rand() POSIX calls and collects them into a list. It turns out the answer with the cleanest and fastest code. (No, 10x code to accelerate 1% is not worth it.)

Answer 1

Here is Mersenne, which unexpectedly seemed faster than MWC and beats C ++, although we are on different computers ;-). It is tempting to see how much parallelizes it, but I better get back to work.

 {-# LANGUAGE BangPatterns #-} {-# OPTIONS_GHC -Wall #-} {-# OPTIONS_GHC -fno-warn-name-shadowing #-} {-# OPTIONS_GHC -fno-warn-type-defaults #-} import System.Random.Mersenne.Pure64 testUniform :: Int -> Double -> PureMT -> Double testUniform 0 !x _ = x testUniform n !x gen = testUniform (n - 1) (x + y) gen' where (y, gen') = randomDouble gen n :: Int n = 10^7 total :: Double total = testUniform n 0 (pureMT $ fromIntegral arbSeed) arbSeed :: Int arbSeed = 8 mean :: Double mean = total / fromIntegral n main :: IO () main = print mean ~/Dropbox/Private/Stochastic $ ./MersennePure +RTS -s 0.4999607889729769 802,924,992 bytes allocated in the heap 164,240 bytes copied during GC 44,312 bytes maximum residency (2 sample(s)) 21,224 bytes maximum slop 1 MB total memory in use (0 MB lost due to fragmentation) Tot time (elapsed) Avg pause Max pause Gen 0 1634 colls, 0 par 0.00s 0.01s 0.0000s 0.0000s Gen 1 2 colls, 0 par 0.00s 0.00s 0.0001s 0.0002s INIT time 0.00s ( 0.00s elapsed) MUT time 0.11s ( 0.11s elapsed) GC time 0.00s ( 0.01s elapsed) EXIT time 0.00s ( 0.00s elapsed) Total time 0.12s ( 0.12s elapsed) %GC time 4.2% (5.4% elapsed) Alloc rate 7,336,065,126 bytes per MUT second Productivity 95.7% of total user, 93.5% of total elapsed