diff --git a/System/Random/MWC/Distributions.hs b/System/Random/MWC/Distributions.hs --- a/System/Random/MWC/Distributions.hs +++ b/System/Random/MWC/Distributions.hs @@ -47,7 +47,20 @@ {-# INLINE normal #-} -- We express standard in terms of normal and not other way round -- because of bug in GHC. See bug #16 for more details. -normal m s gen = loop >>= (\x -> return $! m + s * x) +normal m s gen = do + x <- standard gen + return $! m + s * x + +-- | Generate a normally distributed random variate with zero mean and +-- unit variance. +-- +-- The implementation uses Doornik's modified ziggurat algorithm. +-- Compared to the ziggurat algorithm usually used, this is slower, +-- but generates more independent variates that pass stringent tests +-- of randomness. +standard :: PrimMonad m => Gen (PrimState m) -> m Double +{-# INLINE standard #-} +standard gen = loop where loop = do u <- (subtract 1 . (*2)) `liftM` uniform gen @@ -67,35 +80,35 @@ if e + c * (d - e) < 1 then return x else loop - blocks = (`I.snoc` 0) . I.cons (v/f) . I.cons r . I.unfoldrN 126 go $! T r f - where - go (T b g) = let !u = T h (exp (-0.5 * h * h)) - h = sqrt (-2 * log (v / b + g)) - in Just (h, u) - v = 9.91256303526217e-3 - f = exp (-0.5 * r * r) - {-# NOINLINE blocks #-} - r = 3.442619855899 - ratios = I.zipWith (/) (I.tail blocks) blocks - {-# NOINLINE ratios #-} normalTail neg = tailing where tailing = do - x <- ((/r) . log) `liftM` uniform gen - y <- log `liftM` uniform gen + x <- ((/rNorm) . log) `liftM` uniform gen + y <- log `liftM` uniform gen if y * (-2) < x * x then tailing - else return $! if neg then x - r else r - x + else return $! if neg then x - rNorm else rNorm - x --- | Generate a normally distributed random variate with zero mean and --- unit variance. --- --- The implementation uses Doornik's modified ziggurat algorithm. --- Compared to the ziggurat algorithm usually used, this is slower, --- but generates more independent variates that pass stringent tests --- of randomness. -standard :: PrimMonad m => Gen (PrimState m) -> m Double -{-# INLINE standard #-} -standard = normal 0 1 +-- Constants used by standard/normal. They are floated to the top +-- level to avoid performance regression (Bug #16) when blocks/ratios +-- are recalculated on each call to standard/normal. It's also +-- somewhat difficult to trigger reliably. +blocks :: I.Vector Double +blocks = (`I.snoc` 0) . I.cons (v/f) . I.cons rNorm . I.unfoldrN 126 go $! T rNorm f + where + go (T b g) = let !u = T h (exp (-0.5 * h * h)) + h = sqrt (-2 * log (v / b + g)) + in Just (h, u) + v = 9.91256303526217e-3 + f = exp (-0.5 * rNorm * rNorm) +{-# NOINLINE blocks #-} + +rNorm :: Double +rNorm = 3.442619855899 + +ratios :: I.Vector Double +ratios = I.zipWith (/) (I.tail blocks) blocks +{-# NOINLINE ratios #-} + -- | Generate an exponentially distributed random variate. diff --git a/benchmarks/Benchmark.hs b/benchmarks/Benchmark.hs --- a/benchmarks/Benchmark.hs +++ b/benchmarks/Benchmark.hs @@ -1,4 +1,6 @@ +{-# LANGUAGE BangPatterns #-} import Control.Exception +import Control.Monad import Control.Monad.ST import Criterion.Main import Data.Int @@ -57,6 +59,19 @@ , bgroup "D" [ bench "standard" (standard mwc :: IO Double) , bench "normal" (normal 1 3 mwc :: IO Double) + -- Regression tests for #16. These functions should take 10x + -- longer to execute. + -- + -- N.B. Bang patterns are necessary to trigger the bug with + -- GHC 7.6 + , bench "standard/N" (replicateM_ 10 $ do + !_ <- standard mwc :: IO Double + return () + ) + , bench "normal/N" (replicateM_ 10 $ do + !_ <- normal 1 3 mwc :: IO Double + return () + ) , bench "exponential" (exponential 3 mwc :: IO Double) , bench "gamma,a<1" (gamma 0.5 1 mwc :: IO Double) , bench "gamma,a>1" (gamma 2 1 mwc :: IO Double) diff --git a/benchmarks/tsts.hs b/benchmarks/tsts.hs new file mode 100644 --- /dev/null +++ b/benchmarks/tsts.hs @@ -0,0 +1,65 @@ +{-# LANGUAGE BangPatterns #-} +import Control.Monad +import System.Random.MWC +import System.Random.MWC.Distributions + +main = do + withSystemRandom $ \g -> replicateM_ (300*1000) $ do + -- + !n <- normal 0 1 g + !n <- normal 0 2 g + !n <- normal 3 3 g + !n <- normal 2 4 g + !n <- normal 2 5 g + !n <- normal 1 6 g + !n <- normal 3 7 g + !n <- normal 3 8 g + !n <- normal 3 9 g + !n <- normal 3 10 g + -- + !n <- normal 0 1 g + !n <- normal 0 2 g + !n <- normal 3 3 g + !n <- normal 2 4 g + !n <- normal 2 5 g + !n <- normal 1 6 g + !n <- normal 3 7 g + !n <- normal 3 8 g + !n <- normal 3 9 g + !n <- normal 3 10 g + -- + !n <- normal 0 1 g + !n <- normal 0 2 g + !n <- normal 3 3 g + !n <- normal 2 4 g + !n <- normal 2 5 g + !n <- normal 1 6 g + !n <- normal 3 7 g + !n <- normal 3 8 g + !n <- normal 3 9 g + !n <- normal 3 10 g + -- + !n <- normal 0 1 g + !n <- normal 0 2 g + !n <- normal 3 3 g + !n <- normal 2 4 g + !n <- normal 2 5 g + !n <- normal 1 6 g + !n <- normal 3 7 g + !n <- normal 3 8 g + !n <- normal 3 9 g + !n <- normal 3 10 g + -- + !n <- normal 0 1 g + !n <- normal 0 2 g + !n <- normal 3 3 g + !n <- normal 2 4 g + !n <- normal 2 5 g + !n <- normal 1 6 g + !n <- normal 3 7 g + !n <- normal 3 8 g + !n <- normal 3 9 g + !n <- normal 3 10 g + -- + return () :: IO () + diff --git a/mwc-random.cabal b/mwc-random.cabal --- a/mwc-random.cabal +++ b/mwc-random.cabal @@ -1,5 +1,5 @@ name: mwc-random -version: 0.13.1.0 +version: 0.13.1.1 synopsis: Fast, high quality pseudo random number generation description: This package contains code for generating high quality random