mwc-random 0.10.0.1 → 0.11.0.0

raw patch · 7 files changed

+469/−148 lines, 7 filesdep +HUnitdep +QuickCheckdep +mwc-randomdep ~basedep ~vectorPVP ok

version bump matches the API change (PVP)

Dependencies added: HUnit, QuickCheck, mwc-random, statistics, test-framework, test-framework-hunit, test-framework-quickcheck2

Dependency ranges changed: base, vector

API changes (from Hackage documentation)

Files

• System/Random/MWC.hs +106/−132
• System/Random/MWC/Distributions.hs +171/−0
• benchmarks/Benchmark.hs +42/−13
• mwc-random.cabal +29/−3
• test/run-dieharder-test.sh +26/−0
• test/visual.R +62/−0
• test/visual.hs +33/−0

@@ -2,7 +2,7 @@     MagicHash, Rank2Types, ScopedTypeVariables, TypeFamilies, UnboxedTuples #-} -- | -- Module    : System.Random.MWC--- Copyright : (c) 2009, 2010, 2011 Bryan O'Sullivan+-- Copyright : (c) 2009-2012 Bryan O'Sullivan -- License   : BSD3 -- -- Maintainer  : bos@serpentine.com@@ -10,34 +10,75 @@ -- Portability : portable -- -- Pseudo-random number generation.  This module contains code for--- generating high quality random numbers that follow either a uniform--- or normal distribution.+-- generating high quality random numbers that follow a uniform+-- distribution. --+-- For non-uniform distributions, see the+-- 'System.Random.MWC.Distributions' module.+-- -- The uniform PRNG uses Marsaglia's MWC256 (also known as MWC8222) -- multiply-with-carry generator, which has a period of 2^8222 and -- fares well in tests of randomness.  It is also extremely fast, -- between 2 and 3 times faster than the Mersenne Twister.+--+-- The generator state is stored in the 'Gen' data type. It can be+-- created in several ways:+--+--   1. Using the 'withSystemRandom' call, which creates a random state.+--+--   2. Supply your own seed to 'initialize' function.+--+--   3. Finally, 'create' makes a generator from a fixed seed.+--      Generators created in this way aren't really random.+--+-- For repeatability, the state of the generator can be snapshotted+-- and replayed using the 'save' and 'restore' functions.+--+-- The simplest use is to generate a vector of uniformly distributed values:+--+-- @+--   vs <- withSystemRandom (uniformVector 100)+-- @+--+-- These values can be of any type which is an instance of the class 'Variate'.+--+-- To generate random values on demand, first 'create' a random number generator.+--+-- @+--   gen <- create+-- @+--+-- Keep this generator and use it wherever random values are required. Get a random+-- value using 'uniform' or 'uniformR':+--+-- @+--   v <- uniform gen+-- @+--+-- @+--   v <- uniformR (1, 52) gen+-- @ module System.Random.MWC     (-    -- * Types+    -- * Gen: Pseudo-Random Number Generators       Gen     , GenIO     , GenST-    , Seed-    , fromSeed-    , toSeed-    , Variate(..)-    -- * Other distributions-    , normal-    -- * Creation     , create     , initialize     , withSystemRandom-    -- * State management++    -- * Variates: uniformly distributed values+    , Variate(..)+    , uniformVector++    -- * Seed: state management+    , Seed+    , fromSeed+    , toSeed     , save     , restore-    -- * Helper functions-    , uniformVector+     -- * References     -- $references     ) where@@ -56,24 +97,20 @@ import Data.Ratio              ((%), numerator) import Data.Time.Clock.POSIX   (getPOSIXTime) import Data.Typeable           (Typeable)-import Data.Vector.Generic     (Vector, unsafeFreeze)+import Data.Vector.Generic     (Vector) import Data.Word               (Word, Word8, Word16, Word32, Word64) import Foreign.Marshal.Alloc   (allocaBytes) import Foreign.Marshal.Array   (peekArray) import Prelude hiding (catch) import qualified Data.Vector.Generic         as G-import qualified Data.Vector.Generic.Mutable as GM import qualified Data.Vector.Unboxed         as I import qualified Data.Vector.Unboxed.Mutable as M import System.CPUTime   (cpuTimePrecision, getCPUTime) import System.IO        (IOMode(..), hGetBuf, hPutStrLn, stderr, withBinaryFile) import System.IO.Unsafe (unsafePerformIO) --- FIXME: removal of Unbox constraint leads to severe (~10x)---        performance drop with GHC 6.12. For details see bug #33 in the ---        vector bug tracker[1]--- [1] http://trac.haskell.org/vector/ticket/33 + -- | The class of types for which we can generate uniformly -- distributed random variates. --@@ -84,7 +121,7 @@ -- -- /Note/: Marsaglia's PRNG is not known to be cryptographically -- secure, so you should not use it for cryptographic operations.-class M.Unbox a => Variate a where+class Variate a where     -- | Generate a single uniformly distributed random variate.  The     -- range of values produced varies by type:     --@@ -111,49 +148,49 @@  instance Variate Int8 where     uniform  = uniform1 fromIntegral-    uniformR = uniformRange+    uniformR a b = uniformRange a b     {-# INLINE uniform  #-}     {-# INLINE uniformR #-}  instance Variate Int16 where     uniform  = uniform1 fromIntegral-    uniformR = uniformRange+    uniformR a b = uniformRange a b     {-# INLINE uniform  #-}     {-# INLINE uniformR #-}  instance Variate Int32 where     uniform  = uniform1 fromIntegral-    uniformR = uniformRange+    uniformR a b = uniformRange a b     {-# INLINE uniform  #-}     {-# INLINE uniformR #-}  instance Variate Int64 where     uniform  = uniform2 wordsTo64Bit-    uniformR = uniformRange+    uniformR a b = uniformRange a b     {-# INLINE uniform  #-}     {-# INLINE uniformR #-}  instance Variate Word8 where     uniform  = uniform1 fromIntegral-    uniformR = uniformRange+    uniformR a b = uniformRange a b     {-# INLINE uniform  #-}     {-# INLINE uniformR #-}  instance Variate Word16 where     uniform  = uniform1 fromIntegral-    uniformR = uniformRange+    uniformR a b = uniformRange a b     {-# INLINE uniform  #-}     {-# INLINE uniformR #-}  instance Variate Word32 where     uniform  = uniform1 fromIntegral-    uniformR = uniformRange+    uniformR a b = uniformRange a b     {-# INLINE uniform  #-}     {-# INLINE uniformR #-}  instance Variate Word64 where     uniform  = uniform2 wordsTo64Bit-    uniformR = uniformRange+    uniformR a b = uniformRange a b     {-# INLINE uniform  #-}     {-# INLINE uniformR #-} @@ -184,7 +221,7 @@ #else     uniform = uniform2 wordsTo64Bit #endif-    uniformR = uniformRange+    uniformR a b = uniformRange a b     {-# INLINE uniform  #-}     {-# INLINE uniformR #-} @@ -194,7 +231,7 @@ #else     uniform = uniform2 wordsTo64Bit #endif-    uniformR = uniformRange+    uniformR a b = uniformRange a b     {-# INLINE uniform  #-}     {-# INLINE uniformR #-} @@ -287,10 +324,12 @@ -- verbatim, then its elements are 'xor'ed against elements of the -- default seed until 256 elements are reached. ----- If a seed contains exactly 258 elements last two elements are used--- to set generator state. It's to ensure that @gen' == gen@+-- If a seed contains exactly 258 elements, then the last two elements+-- are used to set the generator's initial state. This allows for+-- complete generator reproducibility, so that e.g. @gen' == gen@ in+-- the following example: ----- > gen' <- initialize . fromSeed =<< save+-- @gen' <- 'initialize' . 'fromSeed' =<< 'save'@ initialize :: (PrimMonad m, Vector v Word32) =>               v Word32 -> m (Gen (PrimState m)) initialize seed = do@@ -330,29 +369,14 @@ toSeed :: (Vector v Word32) => v Word32 -> Seed toSeed v = Seed $ I.create $ do { Gen q <- initialize v; return q } --- Safe version of unsafeFreeze.--- NOTE: vector-0.7 will provide function `freeze' with same---       functionality. This function shall be removed when support for---       vector<=0.6 is dropped-safeFreeze :: (PrimMonad m, Vector v a) => G.Mutable v (PrimState m) a -> m (v a)-safeFreeze v = do-  v' <- GM.unsafeNew (GM.length v)-  GM.unsafeCopy v' v-  unsafeFreeze v'- -- | Save the state of a 'Gen', for later use by 'restore'. save :: PrimMonad m => Gen (PrimState m) -> m Seed-save (Gen q) = Seed `liftM` safeFreeze q+save (Gen q) = Seed `liftM` G.freeze q {-# INLINE save #-} --- NOTE: with vector-0.7 all code could be replaced with `clone' -- | Create a new 'Gen' that mirrors the state of a saved 'Seed'. restore :: PrimMonad m => Seed -> m (Gen (PrimState m))-restore (Seed s) = M.unsafeNew n >>= fill-  where fill q = go 0 where-          go !i | i >= n    = return $! Gen q-                | otherwise = M.unsafeWrite q i (I.unsafeIndex s i) >> go (i+1)-        n = I.length s+restore (Seed s) = Gen `liftM` G.thaw s {-# INLINE restore #-}  @@ -379,6 +403,10 @@ -- pseudo-random numbers (\"\/dev\/urandom\" on Unix-like systems), -- then run the given action. --+-- This is a heavyweight function, intended to be called only+-- occasionally (e.g. once per thread).  You should use the `Gen` it+-- creates to generate many random numbers.+-- -- /Note/: on Windows, this code does not yet use the native -- Cryptographic API as a source of random numbers (it uses the system -- clock instead). As a result, the sequences it generates may not be@@ -404,16 +432,16 @@     where j = fromIntegral (i+1) :: Word8 {-# INLINE nextIndex #-} -a :: Word64-a = 1540315826-{-# INLINE a #-}+aa :: Word64+aa = 1540315826+{-# INLINE aa #-}  uniformWord32 :: PrimMonad m => Gen (PrimState m) -> m Word32 uniformWord32 (Gen q) = do   i  <- nextIndex `liftM` M.unsafeRead q ioff   c  <- fromIntegral `liftM` M.unsafeRead q coff   qi <- fromIntegral `liftM` M.unsafeRead q i-  let t  = a * qi + c+  let t  = aa * qi + c       c' = fromIntegral (t `shiftR` 32)       x  = fromIntegral t + c'       (# x', c'' #)  | x < c'    = (# x + 1, c' + 1 #)@@ -437,12 +465,12 @@   c  <- fromIntegral `liftM` M.unsafeRead q coff   qi <- fromIntegral `liftM` M.unsafeRead q i   qj <- fromIntegral `liftM` M.unsafeRead q j-  let t   = a * qi + c+  let t   = aa * qi + c       c'  = fromIntegral (t `shiftR` 32)       x   = fromIntegral t + c'       (# x', c'' #)  | x < c'    = (# x + 1, c' + 1 #)                      | otherwise = (# x,     c' #)-      u   = a * qj + fromIntegral c''+      u   = aa * qj + fromIntegral c''       d'  = fromIntegral (u `shiftR` 32)       y   = fromIntegral u + d'       (# y', d'' #)  | y < d'    = (# y + 1, d' + 1 #)@@ -475,31 +503,36 @@ -- unsigned data type of same size sub :: (Integral a, Integral (Unsigned a)) => a -> a -> Unsigned a sub x y = fromIntegral x - fromIntegral y+{-# INLINE sub #-}  add :: (Integral a, Integral (Unsigned a)) => a -> Unsigned a -> a add m x = m + fromIntegral x---- Generate uniform value in the range [0,n). Values must be--- unsigned. Second parameter is random number generator-unsignedRange :: (PrimMonad m, Integral a, Bounded a) => a -> m a -> m a-unsignedRange n rnd = go-  where-    buckets = maxBound `div` n-    maxN    = buckets * n-    go = do x <- rnd-            if x < maxN then return (x `div` buckets)-                        else go-{-# INLINE unsignedRange #-}+{-# INLINE add #-} --- Generate unformly distributed value in inclusive range.+-- Generate uniformly distributed value in inclusive range.+--+-- NOTE: This function must be fully applied. Otherwise it won't be+--       inlined, which will cause a severe performance loss.+--+-- > uniformR     = uniformRange      -- won't be inlined+-- > uniformR a b = uniformRange a b  -- will be inlined uniformRange :: ( PrimMonad m                 , Integral a, Bounded a, Variate a                 , Integral (Unsigned a), Bounded (Unsigned a), Variate (Unsigned a))              => (a,a) -> Gen (PrimState m) -> m a uniformRange (x1,x2) g-  | x1 == minBound && x2 == maxBound = uniform g-  | otherwise                        = do x <- unsignedRange (sub x2 x1 + 1) (uniform g)-                                          return $! add x1 x+  | n == 0    = uniform g   -- Abuse overflow in unsigned types+  | otherwise = loop+  where+    -- Allow ranges where x2<x1+    (# i, j #) | x1 < x2   = (# x1, x2 #)+               | otherwise = (# x2, x1 #)+    n       = 1 + sub j i+    buckets = maxBound `div` n+    maxN    = buckets * n+    loop    = do x <- uniform g+                 if x < maxN then return $! add i (x `div` buckets)+                             else loop {-# INLINE uniformRange #-}  -- | Generate a vector of pseudo-random variates.  This is not@@ -510,57 +543,7 @@ uniformVector gen n = G.replicateM n (uniform gen) {-# INLINE uniformVector #-} -data T = T {-# UNPACK #-} !Double {-# UNPACK #-} !Double --- | Generate a normally distributed random variate.------ 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.-normal :: PrimMonad m => Gen (PrimState m) -> m Double-normal gen = loop-  where-    loop = do-      u  <- (subtract 1 . (*2)) `liftM` uniform gen-      ri <- uniform gen-      let i  = fromIntegral ((ri :: Word32) .&. 127)-          bi = I.unsafeIndex blocks i-          bj = I.unsafeIndex blocks (i+1)-      if abs u < I.unsafeIndex ratios i-        then return $! u * bi-        else if i == 0-        then normalTail (u < 0)-        else do-          let x  = u * bi-              xx = x * x-              d  = exp (-0.5 * (bi * bi - xx))-              e  = exp (-0.5 * (bj * bj - xx))-          c <- uniform gen-          if e + c * (d - e) < 1-            then return x-            else loop-    blocks = let f = exp (-0.5 * r * r)-             in (`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-    {-# 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-              if y * (-2) < x * x-                then tailing-                else return $! if neg then x - r else r - x-{-# INLINE normal #-}- defaultSeed :: I.Vector Word32 defaultSeed = I.fromList [   0x7042e8b3, 0x06f7f4c5, 0x789ea382, 0x6fb15ad8, 0x54f7a879, 0x0474b184,@@ -609,15 +592,6 @@ {-# NOINLINE defaultSeed #-}  -- $references------ * Doornik, J.A. (2005) An improved ziggurat method to generate---   normal random samples. Mimeo, Nuffield College, University of---   Oxford.  <http://www.doornik.com/research/ziggurat.pdf>------ * Doornik, J.A. (2007) Conversion of high-period random numbers to---   floating point.---   /ACM Transactions on Modeling and Computer Simulation/ 17(1).---   <http://www.doornik.com/research/randomdouble.pdf> -- -- * Marsaglia, G. (2003) Seeds for random number generators. --   /Communications of the ACM/ 46(5):90&#8211;93.

@@ -0,0 +1,171 @@+{-# LANGUAGE BangPatterns #-}+-- |+-- Module    : System.Random.MWC.Distributions+-- Copyright : (c) 2012 Bryan O'Sullivan+-- License   : BSD3+--+-- Maintainer  : bos@serpentine.com+-- Stability   : experimental+-- Portability : portable+--+-- Pseudo-random number generation for non-uniform distributions.++module System.Random.MWC.Distributions +    (+    -- * Variates: non-uniformly distributed values+      normal+    , standard+    , exponential+    , gamma+    , chiSquare++    -- * References+    -- $references+    ) where++import Control.Monad (liftM)+import Control.Monad.Primitive (PrimMonad, PrimState)+import Data.Bits ((.&.))+import Data.Word (Word32)+import System.Random.MWC (Gen, uniform)+import qualified Data.Vector.Unboxed as I++-- Unboxed 2-tuple+data T = T {-# UNPACK #-} !Double {-# UNPACK #-} !Double+++-- | Generate a normally distributed random variate with given mean+-- and standard deviation.+normal :: PrimMonad m+       => Double                -- ^ Mean+       -> Double                -- ^ Standard deviation+       -> Gen (PrimState m)+       -> m Double+{-# INLINE normal #-}+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+      ri <- uniform gen+      let i  = fromIntegral ((ri :: Word32) .&. 127)+          bi = I.unsafeIndex blocks i+          bj = I.unsafeIndex blocks (i+1)+      case () of+        _| abs u < I.unsafeIndex ratios i -> return $! u * bi+         | i == 0                         -> normalTail (u < 0)+         | otherwise                      -> do+             let x  = u * bi+                 xx = x * x+                 d  = exp (-0.5 * (bi * bi - xx))+                 e  = exp (-0.5 * (bj * bj - xx))+             c <- uniform gen+             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+              if y * (-2) < x * x+                then tailing+                else return $! if neg then x - r else r - x+++-- | Generate exponentially distributed random variate.+exponential :: PrimMonad m+            => Double            -- ^ Scale parameter+            -> Gen (PrimState m) -- ^ Generator+            -> m Double+{-# INLINE exponential #-}+exponential beta gen = do+  x <- uniform gen+  return $! - log x / beta+++-- | Random variate generator for gamma distribution.+gamma :: PrimMonad m+      => Double                 -- ^ Shape parameter+      -> Double                 -- ^ Scale parameter+      -> Gen (PrimState m)      -- ^ Generator+      -> m Double+{-# INLINE gamma #-}+gamma a b gen+  | a <= 0    = pkgError "gamma" "negative alpha parameter"+  | otherwise = mainloop+    where+      mainloop = do+        T x v <- innerloop+        u     <- uniform gen+        let cont =  u > 1 - 0.331 * sqr (sqr x)+                 && log u > 0.5 * sqr x + a1 * (1 - v + log v) -- Rarely evaluated+        case () of+          _| cont      -> mainloop+           | a >= 1    -> return $! a1 * v * b+           | otherwise -> do y <- uniform gen+                             return $! y ** (1 / a) * a1 * v * b+      -- inner loop+      innerloop = do+        x <- standard gen+        case 1 + a2*x of+          v | v <= 0    -> innerloop+            | otherwise -> return $! T x (v*v*v)+      -- constants+      a' = if a < 1 then a + 1 else a+      a1 = a' - 1/3+      a2 = 1 / sqrt(9 * a1)+++-- | Random variate generator for chi square distribution.+chiSquare :: PrimMonad m+          => Int                -- ^ Number of degrees of freedom+          -> Gen (PrimState m)  -- ^ Generator+          -> m Double+{-# INLINE chiSquare #-}+chiSquare n gen+  | n <= 0    = pkgError "chiSquare" "number of degrees of freedom must be positive"+  | otherwise = do x <- gamma (0.5 * fromIntegral n) 1 gen+                   return $! 2 * x+++sqr :: Double -> Double+sqr x = x * x+{-# INLINE sqr #-}++pkgError :: String -> String -> a+pkgError func msg = error $ "System.Random.MWC.Distributions." ++ func +++                            ": " ++ msg++-- $references+--+-- * Doornik, J.A. (2005) An improved ziggurat method to generate+--   normal random samples. Mimeo, Nuffield College, University of+--   Oxford.  <http://www.doornik.com/research/ziggurat.pdf>+--+-- * Doornik, J.A. (2007) Conversion of high-period random numbers to+--   floating point.+--   /ACM Transactions on Modeling and Computer Simulation/ 17(1).+--   <http://www.doornik.com/research/randomdouble.pdf>

@@ -5,6 +5,7 @@ import Data.Word import qualified System.Random as R import System.Random.MWC+import System.Random.MWC.Distributions import qualified System.Random.Mersenne as M  main = do@@ -12,19 +13,47 @@   mtg <- M.newMTGen . Just =<< uniform mwc   defaultMain      [ bgroup "mwc"-      [ bench "Double"  (uniform mwc :: IO Double)-      , bench "Int"     (uniform mwc :: IO Int)-      , bench "Int8"    (uniform mwc :: IO Int8)-      , bench "Int16"   (uniform mwc :: IO Int16)-      , bench "Int32"   (uniform mwc :: IO Int32)-      , bench "Int64"   (uniform mwc :: IO Int64)-      , bench "Word"    (uniform mwc :: IO Word)-      , bench "Word8"   (uniform mwc :: IO Word8)-      , bench "Word16"  (uniform mwc :: IO Word16)-      , bench "Word32"  (uniform mwc :: IO Word32)-      , bench "Word64"  (uniform mwc :: IO Word64)-      , bench "Integer" (uniform mwc :: IO Word64)-      , bench "normal"  (normal mwc :: IO Double)+      -- One letter group names are used so they will fit on the plot.+      --+      --  U - uniform+      --  R - uniformR+      --  D - distribution+      [ bgroup "U"+        [ bench "Double"  (uniform mwc :: IO Double)+        , bench "Int"     (uniform mwc :: IO Int)+        , bench "Int8"    (uniform mwc :: IO Int8)+        , bench "Int16"   (uniform mwc :: IO Int16)+        , bench "Int32"   (uniform mwc :: IO Int32)+        , bench "Int64"   (uniform mwc :: IO Int64)+        , bench "Word"    (uniform mwc :: IO Word)+        , bench "Word8"   (uniform mwc :: IO Word8)+        , bench "Word16"  (uniform mwc :: IO Word16)+        , bench "Word32"  (uniform mwc :: IO Word32)+        , bench "Word64"  (uniform mwc :: IO Word64)+        ]+      , bgroup "R"+        -- I'm not entirely convinced that this is right way to test+        -- uniformR. /A.Khudyakov/+        [ bench "Double"  (uniformR (-3.21,26) mwc :: IO Double)+        , bench "Int"     (uniformR (-12,679)  mwc :: IO Int)+        , bench "Int8"    (uniformR (-12,4)    mwc :: IO Int8)+        , bench "Int16"   (uniformR (-12,679)  mwc :: IO Int16)+        , bench "Int32"   (uniformR (-12,679)  mwc :: IO Int32)+        , bench "Int64"   (uniformR (-12,679)  mwc :: IO Int64)+        , bench "Word"    (uniformR (34,633)   mwc :: IO Word)+        , bench "Word8"   (uniformR (34,63)    mwc :: IO Word8)+        , bench "Word16"  (uniformR (34,633)   mwc :: IO Word16)+        , bench "Word32"  (uniformR (34,633)   mwc :: IO Word32)+        , bench "Word64"  (uniformR (34,633)   mwc :: IO Word64)+        ]+      , bgroup "D"+        [ bench "standard"    (standard      mwc :: IO Double)+        , bench "normal"      (normal 1 3    mwc :: IO Double)+        , 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)+        , bench "chiSquare"   (chiSquare 4   mwc :: IO Double)+        ]       ]     , bgroup "random"       [

@@ -1,5 +1,5 @@ name:           mwc-random-version:        0.10.0.1+version:        0.11.0.0 synopsis:       Fast, high quality pseudo random number generation description:   This package contains code for generating high quality random@@ -23,21 +23,25 @@ copyright:      2009, 2010, 2011 Bryan O'Sullivan category:       Math, Statistics build-type:     Simple-cabal-version:  >= 1.6+cabal-version:  >= 1.8 extra-source-files:   README.markdown   benchmarks/*.hs   benchmarks/Quickie.hs   benchmarks/mwc-random-benchmarks.cabal+  test/*.R+  test/*.sh+  test/visual.hs  library   exposed-modules:     System.Random.MWC+    System.Random.MWC.Distributions   build-depends:     base < 5,     primitive,     time,-    vector >= 0.6.0.2+    vector >= 0.7   if impl(ghc >= 6.10)     build-depends:       base >= 4@@ -48,6 +52,28 @@   ghc-options: -Wall -funbox-strict-fields   if impl(ghc >= 6.8)     ghc-options: -fwarn-tabs++test-suite tests+  buildable:      False+  type:           exitcode-stdio-1.0+  hs-source-dirs: test+  main-is:        tests.hs+  other-modules:  KS+                  QC+                  Uniform++  ghc-options:+    -Wall -threaded -rtsopts++  build-depends:+    HUnit,+    QuickCheck,+    base,+    mwc-random,+    statistics >= 0.10.1.0,+    test-framework,+    test-framework-hunit,+    test-framework-quickcheck2  source-repository head   type:     git

@@ -0,0 +1,26 @@+#!/bin/sh+#+# Run dieharder set of tests for PRNG. All command line parameters are+# passed directly to the dieharder. If no parameters are given -a flag+# is passed which runs all available tests. Full list of dieharder+# options is available at dieharder manpage+#+# NOTE:+#   Full set of test require a lot of time to complete. From several+#   hours to a few days depending on CPU speed and thoroughness+#   settings.+#+# dieharder-source.hs is enthropy source for this test.+#+# This test require dieharder to be installed. It is available at:+#   http://www.phy.duke.edu/~rgb/General/dieharder.php++which dieharder > /dev/null || { echo "dieharder is not found. Aborting"; exit 1; }++ghc -fforce-recomp -O2 diehard-source+(+    date+    ./diehard-source | \+	if [ $# = 0 ]; then dieharder -a -g 200; else dieharder "$@" -g 200; fi+    date+) | tee diehard.log

@@ -0,0 +1,62 @@+# Ugly script for displaying distributions alogside with theoretical+# distribution.+++view.dumps <- function() {+  load.d <- function(name) read.table(name)[,1]+  plot.d <- function(name, dens, rng) {+    smp <- load.d( name )+    plot( density(smp), xlim=rng, main=name, col='blue', lwd=2)+    hist( smp, probability=TRUE, breaks=100, add=TRUE)+    plot( dens, xlim=rng, col='red', add=TRUE, lwd=2)+  }+  ################################################################+  # Normal+  plot.d ("distr/normal-0-1",+          function(x) dnorm( x, 0, 1 ),+          c(-4,4) )+  readline()+  # +  plot.d ("distr/normal-1-2",+          function(x) dnorm( x, 1, 2 ),+          c(-6,8) )+  readline();++  ################################################################+  # Gamma+  plot.d ("distr/gamma-1.0-1.0",+          function(x) dgamma( x, 1, 1 ),+          c(-1,8) )+  readline();+  #+  plot.d ("distr/gamma-0.3-0.4",+          function(x) dgamma( x, 0.3, scale=0.4 ),+          c(-0.25,2) )+  readline();+  #+  plot.d ("distr/gamma-0.3-3.0",+          function(x) dgamma( x, 0.3, scale=3.0 ),+          c(-1,5) )+  readline();+  #+  plot.d ("distr/gamma-3.0-0.4",+          function(x) dgamma( x, 3.0, scale=0.4 ),+          c(-1,6) )+  readline();+  #+  plot.d ("distr/gamma-3.0-3.0",+          function(x) dgamma( x, 3.0, scale=3.0 ),+          c(-1,32) )+  readline();+  ################################################################+  # Exponential+  plot.d ("distr/exponential-1",+          function(x) dexp(x,1),+          c(-0.5, 9) )+  readline()+  #+  plot.d ("distr/exponential-3",+          function(x) dexp(x,3),+          c(-0.5, 3) )+  readline()+}

@@ -0,0 +1,33 @@+-- Generates samples of value for display with visual.R+import Control.Monad++import System.Directory  (createDirectoryIfMissing,setCurrentDirectory)+import System.IO++import qualified System.Random.MWC               as MWC+import qualified System.Random.MWC.Distributions as MWC+++dumpSample :: Show a => Int -> FilePath -> IO a -> IO ()+dumpSample n fname gen =+  withFile fname WriteMode $ \h -> +    replicateM_ n (hPutStrLn h . show =<< gen)+  +main :: IO ()+main = MWC.withSystemRandom $ \g -> do+  let n   = 10000+      dir = "distr"+  createDirectoryIfMissing True dir+  setCurrentDirectory           dir+  -- Normal+  dumpSample n "normal-0-1" $ MWC.normal 0 1 g+  dumpSample n "normal-1-2" $ MWC.normal 1 2 g+  -- Gamma+  dumpSample n "gamma-1.0-1.0" $ MWC.gamma  1.0 1.0 g+  dumpSample n "gamma-0.3-0.4" $ MWC.gamma  0.3 0.4 g+  dumpSample n "gamma-0.3-3.0" $ MWC.gamma  0.3 3.0 g+  dumpSample n "gamma-3.0-0.4" $ MWC.gamma  3.0 0.4 g+  dumpSample n "gamma-3.0-3.0" $ MWC.gamma  3.0 3.0 g+  -- Exponential+  dumpSample n "exponential-1" $ MWC.exponential 1 g+  dumpSample n "exponential-3" $ MWC.exponential 3 g