summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAlexeyKhudyakov <>2018-07-11 16:24:00 (GMT)
committerhdiff <hdiff@hdiff.luite.com>2018-07-11 16:24:00 (GMT)
commit5bc8cc08c05484ac89a9b5109cf66cd5cac3e66a (patch)
tree69fe702136a99f6fbcb48a59983e2faa132e5479
parent1d28f0ce0867538527b4d337ad7871e3eead1103 (diff)
version 0.14.0.0HEAD0.14.0.0master
-rw-r--r--README.markdown2
-rw-r--r--System/Random/MWC.hs139
-rw-r--r--System/Random/MWC/SeedSource.hs97
-rw-r--r--benchmarks/Benchmark.hs113
-rw-r--r--benchmarks/Quickie.hs13
-rw-r--r--benchmarks/mwc-random-benchmarks.cabal18
-rw-r--r--changelog.md11
-rw-r--r--mwc-random.cabal60
-rw-r--r--test/KS.hs54
-rw-r--r--test/QC.hs56
-rw-r--r--test/run-dieharder-test.sh26
-rw-r--r--test/tests.hs16
-rw-r--r--test/visual.R105
-rw-r--r--test/visual.hs44
14 files changed, 167 insertions, 587 deletions
diff --git a/README.markdown b/README.markdown
index 9cd884d..975c9f9 100644
--- a/README.markdown
+++ b/README.markdown
@@ -1,4 +1,6 @@
# Efficient, general purpose pseudo-random number generation
+[![Build Status](https://travis-ci.org/Shimuuar/mwc-random.png?branch=master)](https://travis-ci.org/Shimuuar/mwc-random)
+[![Build status](https://ci.appveyor.com/api/projects/status/4228vkxje4as3nhw/branch/master)](https://ci.appveyor.com/project/Shimuuar/mwc-random)
This package provides the System.Random.MWC module, a Haskell library
for generating high-quality pseudo-random numbers in a space- and
diff --git a/System/Random/MWC.hs b/System/Random/MWC.hs
index 9453705..33e102e 100644
--- a/System/Random/MWC.hs
+++ b/System/Random/MWC.hs
@@ -1,6 +1,6 @@
{-# LANGUAGE BangPatterns, CPP, DeriveDataTypeable, FlexibleContexts,
- MagicHash, Rank2Types, ScopedTypeVariables, TypeFamilies, UnboxedTuples,
- ForeignFunctionInterface #-}
+ MagicHash, Rank2Types, ScopedTypeVariables, TypeFamilies, UnboxedTuples
+ #-}
-- |
-- Module : System.Random.MWC
-- Copyright : (c) 2009-2012 Bryan O'Sullivan
@@ -99,35 +99,25 @@ module System.Random.MWC
#endif
import Control.Monad (ap, liftM, unless)
-import Control.Monad.Primitive (PrimMonad, PrimState, unsafePrimToIO)
-#if MIN_VERSION_primitive(0,6,0)
-import Control.Monad.Primitive (PrimBase)
-#endif
+import Control.Monad.Primitive (PrimMonad, PrimBase, PrimState, unsafePrimToIO)
import Control.Monad.ST (ST)
import Data.Bits ((.&.), (.|.), shiftL, shiftR, xor)
import Data.Int (Int8, Int16, Int32, Int64)
import Data.IORef (atomicModifyIORef, newIORef)
-import Data.Ratio ((%), numerator)
-import Data.Time.Clock.POSIX (getPOSIXTime)
import Data.Typeable (Typeable)
import Data.Vector.Generic (Vector)
-import Data.Word (Word8, Word16, Word32, Word64)
-#if !MIN_VERSION_base(4,8,0)
-import Data.Word (Word)
-#endif
-import Foreign.Marshal.Alloc (allocaBytes)
-import Foreign.Marshal.Array (peekArray)
+import Data.Word
import qualified Data.Vector.Generic as G
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 (hPutStrLn, stderr)
import System.IO.Unsafe (unsafePerformIO)
import qualified Control.Exception as E
#if defined(mingw32_HOST_OS)
import Foreign.Ptr
import Foreign.C.Types
#endif
+import System.Random.MWC.SeedSource
-- | The class of types for which we can generate uniformly
@@ -202,7 +192,7 @@ instance Variate Word16 where
{-# INLINE uniformR #-}
instance Variate Word32 where
- uniform = uniform1 fromIntegral
+ uniform = uniform1 id
uniformR a b = uniformRange a b
{-# INLINE uniform #-}
{-# INLINE uniformR #-}
@@ -359,6 +349,17 @@ create = initialize defaultSeed
-- the following example:
--
-- @gen' <- 'initialize' . 'fromSeed' =<< 'save'@
+--
+-- In the MWC algorithm, the /carry/ value must be strictly smaller than the
+-- multiplicator (see https://en.wikipedia.org/wiki/Multiply-with-carry).
+-- Hence, if a seed contains exactly 258 elements, the /carry/ value, which is
+-- the last of the 258 values, is moduloed by the multiplicator.
+--
+-- Note that if the /first/ carry value is strictly smaller than the multiplicator,
+-- all subsequent carry values are also strictly smaller than the multiplicator
+-- (a proof of this is in the comments of the code of 'uniformWord32'), hence
+-- when restoring a saved state, we have the guarantee that moduloing the saved
+-- carry won't modify its value.
initialize :: (PrimMonad m, Vector v Word32) =>
v Word32 -> m (Gen (PrimState m))
initialize seed = do
@@ -367,7 +368,7 @@ initialize seed = do
if fini == 258
then do
M.unsafeWrite q ioff $ G.unsafeIndex seed ioff .&. 255
- M.unsafeWrite q coff $ G.unsafeIndex seed coff
+ M.unsafeWrite q coff $ G.unsafeIndex seed coff `mod` fromIntegral aa
else do
M.unsafeWrite q ioff 255
M.unsafeWrite q coff 362436
@@ -409,70 +410,6 @@ restore (Seed s) = Gen `liftM` G.thaw s
{-# INLINE restore #-}
--- Aquire seed from current time. This is horrible fallback for
--- Windows system.
-acquireSeedTime :: IO [Word32]
-acquireSeedTime = do
- c <- (numerator . (%cpuTimePrecision)) `liftM` getCPUTime
- t <- toRational `liftM` getPOSIXTime
- let n = fromIntegral (numerator t) :: Word64
- return [fromIntegral c, fromIntegral n, fromIntegral (n `shiftR` 32)]
-
--- | Acquire seed from the system entropy source. On Unix machines,
--- this will attempt to use @/dev/urandom@. On Windows, it will internally
--- use @RtlGenRandom@.
-acquireSeedSystem :: IO [Word32]
-acquireSeedSystem = do
-#if !defined(mingw32_HOST_OS)
- -- Read 256 random Word32s from /dev/urandom
- let nbytes = 1024
- random = "/dev/urandom"
- allocaBytes nbytes $ \buf -> do
- nread <- withBinaryFile random ReadMode $
- \h -> hGetBuf h buf nbytes
- peekArray (nread `div` 4) buf
-#else
- let nbytes = 1024
- -- Generate 256 random Word32s from RtlGenRandom
- allocaBytes nbytes $ \buf -> do
- ok <- c_RtlGenRandom buf (fromIntegral nbytes)
- if ok then return () else fail "Couldn't use RtlGenRandom"
- peekArray (nbytes `div` 4) buf
-
--- Note: on 64-bit Windows, the 'stdcall' calling convention
--- isn't supported, so we use 'ccall' instead.
-#if defined(i386_HOST_ARCH)
-# define WINDOWS_CCONV stdcall
-#elif defined(x86_64_HOST_ARCH)
-# define WINDOWS_CCONV ccall
-#else
-# error Unknown mingw32 architecture!
-#endif
-
--- Note: On Windows, the typical convention would be to use
--- the CryptoGenRandom API in order to generate random data.
--- However, here we use 'SystemFunction036', AKA RtlGenRandom.
---
--- This is a commonly used API for this purpose; one bonus is
--- that it avoids having to bring in the CryptoAPI library,
--- and completely sidesteps the initialization cost of CryptoAPI.
---
--- While this function is technically "subject to change" that is
--- extremely unlikely in practice: rand_s in the Microsoft CRT uses
--- this, and they can't change it easily without also breaking
--- backwards compatibility with e.g. statically linked applications.
---
--- The name 'SystemFunction036' is the actual link-time name; the
--- display name is just for giggles, I guess.
---
--- See also:
--- - http://blogs.msdn.com/b/michael_howard/archive/2005/01/14/353379.aspx
--- - https://bugzilla.mozilla.org/show_bug.cgi?id=504270
---
-foreign import WINDOWS_CCONV unsafe "SystemFunction036"
- c_RtlGenRandom :: Ptr a -> CULong -> IO Bool
-#endif
-
-- | Seed a PRNG with data from the system's fast source of
-- pseudo-random numbers (\"@\/dev\/urandom@\" on Unix-like systems or
-- @RtlGenRandom@ on Windows), then run the given action.
@@ -480,22 +417,13 @@ foreign import WINDOWS_CCONV unsafe "SystemFunction036"
-- This is a somewhat expensive function, and is intended to be called
-- only occasionally (e.g. once per thread). You should use the `Gen`
-- it creates to generate many random numbers.
-withSystemRandom ::
-#if MIN_VERSION_primitive(0,6,0)
- PrimBase m
-#else
- PrimMonad m
-#endif
+withSystemRandom :: PrimBase m
=> (Gen (PrimState m) -> m a) -> IO a
withSystemRandom act = do
- seed <- acquireSeedSystem `E.catch` \(_::E.IOException) -> do
+ seed <- acquireSeedSystem 256 `E.catch` \(_::E.IOException) -> do
seen <- atomicModifyIORef warned ((,) True)
unless seen $ E.handle (\(_::E.IOException) -> return ()) $ do
-#if !defined(mingw32_HOST_OS)
- hPutStrLn stderr ("Warning: Couldn't open /dev/urandom")
-#else
- hPutStrLn stderr ("Warning: Couldn't use RtlGenRandom")
-#endif
+ hPutStrLn stderr $ "Warning: Couldn't use randomness source " ++ randomSourceName
hPutStrLn stderr ("Warning: using system clock for seed instead " ++
"(quality will be lower)")
acquireSeedTime
@@ -516,6 +444,10 @@ nextIndex i = fromIntegral j
where j = fromIntegral (i+1) :: Word8
{-# INLINE nextIndex #-}
+-- The multiplicator : 0x5BCF5AB2
+--
+-- Eventhough it is a 'Word64', it is important for the correctness of the proof
+-- on carry value that it is /not/ greater than maxBound 'Word32'.
aa :: Word64
aa = 1540315826
{-# INLINE aa #-}
@@ -526,10 +458,29 @@ uniformWord32 (Gen q) = do
c <- fromIntegral `liftM` M.unsafeRead q coff
qi <- fromIntegral `liftM` M.unsafeRead q i
let t = aa * qi + c
+ -- The comments in this function are a proof that:
+ -- "if the carry value is strictly smaller than the multiplicator,
+ -- the next carry value is also strictly smaller than the multiplicator."
+ -- Eventhough the proof is written in terms of the actual value of the multiplicator,
+ -- it holds for any multiplicator value /not/ greater than maxBound 'Word32'
+ --
+ -- (In the code, the multiplicator is aa, the carry value is c,
+ -- the next carry value is c''.)
+ --
+ -- So we'll assume that c < aa, and show that c'' < aa :
+ --
+ -- by definition, aa = 0x5BCF5AB2, qi <= 0xFFFFFFFF (because it is a 'Word32')
+ -- hence aa*qi <= 0x5BCF5AB200000000 - 0x5BCF5AB2.
+ --
+ -- hence t < 0x5BCF5AB200000000 (because t = aa * qi + c and c < 0x5BCF5AB2)
+ -- hence t <= 0x5BCF5AB1FFFFFFFF
c' = fromIntegral (t `shiftR` 32)
+ -- c' < 0x5BCF5AB1
x = fromIntegral t + c'
(# x', c'' #) | x < c' = (# x + 1, c' + 1 #)
| otherwise = (# x, c' #)
+ -- hence c'' < 0x5BCF5AB2,
+ -- hence c'' < aa, which is what we wanted to prove.
M.unsafeWrite q i x'
M.unsafeWrite q ioff (fromIntegral i)
M.unsafeWrite q coff (fromIntegral c'')
diff --git a/System/Random/MWC/SeedSource.hs b/System/Random/MWC/SeedSource.hs
new file mode 100644
index 0000000..6018ba0
--- /dev/null
+++ b/System/Random/MWC/SeedSource.hs
@@ -0,0 +1,97 @@
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE ForeignFunctionInterface #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+-- |
+-- Low level source of random values for seeds. It should work on both
+-- unices and windows
+module System.Random.MWC.SeedSource (
+ acquireSeedSystem
+ , acquireSeedTime
+ , randomSourceName
+ ) where
+
+import Control.Monad (liftM)
+import Data.Word (Word32,Word64)
+import Data.Bits (shiftR)
+import Data.Ratio ((%), numerator)
+import Data.Time.Clock.POSIX (getPOSIXTime)
+
+import Foreign.Storable
+import Foreign.Marshal.Alloc (allocaBytes)
+import Foreign.Marshal.Array (peekArray)
+#if defined(mingw32_HOST_OS)
+import Foreign.Ptr
+import Foreign.C.Types
+#endif
+import System.CPUTime (cpuTimePrecision, getCPUTime)
+import System.IO (IOMode(..), hGetBuf, withBinaryFile)
+
+-- Acquire seed from current time. This is horrible fallback for
+-- Windows system.
+acquireSeedTime :: IO [Word32]
+acquireSeedTime = do
+ c <- (numerator . (%cpuTimePrecision)) `liftM` getCPUTime
+ t <- toRational `liftM` getPOSIXTime
+ let n = fromIntegral (numerator t) :: Word64
+ return [fromIntegral c, fromIntegral n, fromIntegral (n `shiftR` 32)]
+
+-- | Acquire seed from the system entropy source. On Unix machines,
+-- this will attempt to use @/dev/urandom@. On Windows, it will internally
+-- use @RtlGenRandom@.
+acquireSeedSystem :: forall a. Storable a => Int -> IO [a]
+acquireSeedSystem nElts = do
+ let eltSize = sizeOf (undefined :: a)
+ nbytes = nElts * eltSize
+#if !defined(mingw32_HOST_OS)
+ allocaBytes nbytes $ \buf -> do
+ nread <- withBinaryFile "/dev/urandom" ReadMode $ \h -> hGetBuf h buf nbytes
+ peekArray (nread `div` eltSize) buf
+#else
+ -- Generate 256 random Word32s from RtlGenRandom
+ allocaBytes nbytes $ \buf -> do
+ ok <- c_RtlGenRandom buf (fromIntegral nbytes)
+ if ok then return () else fail "Couldn't use RtlGenRandom"
+ peekArray nElts buf
+
+-- Note: on 64-bit Windows, the 'stdcall' calling convention
+-- isn't supported, so we use 'ccall' instead.
+#if defined(i386_HOST_ARCH)
+# define WINDOWS_CCONV stdcall
+#elif defined(x86_64_HOST_ARCH)
+# define WINDOWS_CCONV ccall
+#else
+# error Unknown mingw32 architecture!
+#endif
+
+-- Note: On Windows, the typical convention would be to use
+-- the CryptoGenRandom API in order to generate random data.
+-- However, here we use 'SystemFunction036', AKA RtlGenRandom.
+--
+-- This is a commonly used API for this purpose; one bonus is
+-- that it avoids having to bring in the CryptoAPI library,
+-- and completely sidesteps the initialization cost of CryptoAPI.
+--
+-- While this function is technically "subject to change" that is
+-- extremely unlikely in practice: rand_s in the Microsoft CRT uses
+-- this, and they can't change it easily without also breaking
+-- backwards compatibility with e.g. statically linked applications.
+--
+-- The name 'SystemFunction036' is the actual link-time name; the
+-- display name is just for giggles, I guess.
+--
+-- See also:
+-- - http://blogs.msdn.com/b/michael_howard/archive/2005/01/14/353379.aspx
+-- - https://bugzilla.mozilla.org/show_bug.cgi?id=504270
+--
+foreign import WINDOWS_CCONV unsafe "SystemFunction036"
+ c_RtlGenRandom :: Ptr a -> CULong -> IO Bool
+#endif
+
+
+-- | Name of source of randomness. It should be used in error messages
+randomSourceName :: String
+#if !defined(mingw32_HOST_OS)
+randomSourceName = "/dev/urandom"
+#else
+randomSourceName = "RtlGenRandom"
+#endif
diff --git a/benchmarks/Benchmark.hs b/benchmarks/Benchmark.hs
deleted file mode 100644
index 9810595..0000000
--- a/benchmarks/Benchmark.hs
+++ /dev/null
@@ -1,113 +0,0 @@
-{-# LANGUAGE BangPatterns #-}
-import Control.Exception
-import Control.Monad
-import Control.Monad.ST
-import Criterion.Main
-import Data.Int
-import Data.Word
-import qualified Data.Vector.Unboxed as U
-import qualified System.Random as R
-import System.Random.MWC
-import System.Random.MWC.Distributions
-import System.Random.MWC.CondensedTable
-import qualified System.Random.Mersenne as M
-
-makeTableUniform :: Int -> CondensedTable U.Vector Int
-makeTableUniform n =
- tableFromProbabilities $ U.zip (U.enumFromN 0 n) (U.replicate n (1 / fromIntegral n))
-{-# INLINE makeTableUniform #-}
-
-
-main = do
- mwc <- create
- mtg <- M.newMTGen . Just =<< uniform mwc
- defaultMain
- [ bgroup "mwc"
- -- 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)
- -- 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)
- , bench "chiSquare" (chiSquare 4 mwc :: IO Double)
- ]
- , bgroup "CT/gen" $ concat
- [ [ bench ("uniform "++show i) (genFromTable (makeTableUniform i) mwc :: IO Int)
- | i <- [2..10]
- ]
- , [ bench ("poisson " ++ show l) (genFromTable (tablePoisson l) mwc :: IO Int)
- | l <- [0.01, 0.2, 0.8, 1.3, 2.4, 8, 12, 100, 1000]
- ]
- , [ bench ("binomial " ++ show p ++ " " ++ show n) (genFromTable (tableBinomial n p) mwc :: IO Int)
- | (n,p) <- [ (4, 0.5), (10,0.1), (10,0.6), (10, 0.8), (100,0.4)]
- ]
- ]
- , bgroup "CT/table" $ concat
- [ [ bench ("uniform " ++ show i) $ whnf makeTableUniform i
- | i <- [2..30]
- ]
- , [ bench ("poisson " ++ show l) $ whnf tablePoisson l
- | l <- [0.01, 0.2, 0.8, 1.3, 2.4, 8, 12, 100, 1000]
- ]
- , [ bench ("binomial " ++ show p ++ " " ++ show n) $ whnf (tableBinomial n) p
- | (n,p) <- [ (4, 0.5), (10,0.1), (10,0.6), (10, 0.8), (100,0.4)]
- ]
- ]
- ]
- , bgroup "random"
- [
- bench "Double" (R.randomIO >>= evaluate :: IO Double)
- , bench "Int" (R.randomIO >>= evaluate :: IO Int)
- ]
- , bgroup "mersenne"
- [
- bench "Double" (M.random mtg :: IO Double)
- , bench "Int" (M.random mtg :: IO Int)
- ]
- ]
diff --git a/benchmarks/Quickie.hs b/benchmarks/Quickie.hs
deleted file mode 100644
index 2418b8f..0000000
--- a/benchmarks/Quickie.hs
+++ /dev/null
@@ -1,13 +0,0 @@
-{-# LANGUAGE BangPatterns #-}
-import System.Random.MWC (create, uniform)
-import Control.Monad.ST (ST, runST)
-
-u :: ST s Double
-u = do
- let last = 1000000 :: Int
- gen <- create
- let loop !n !i | n == last = return i
- | otherwise = uniform gen >>= loop (n+1)
- loop 0 0
-
-main = print (runST u)
diff --git a/benchmarks/mwc-random-benchmarks.cabal b/benchmarks/mwc-random-benchmarks.cabal
deleted file mode 100644
index c4d0dac..0000000
--- a/benchmarks/mwc-random-benchmarks.cabal
+++ /dev/null
@@ -1,18 +0,0 @@
-name: mwc-random-benchmarks
-version: 0
-synopsis: Benchmarks for the mwc-random package
-description: Benchmarks for the mwc-random package
-license: BSD3
-license-file: ../LICENSE
-build-type: Simple
-cabal-version: >= 1.6
-
-executable bm
- main-is: Benchmark.hs
-
- build-depends:
- base < 5,
- criterion,
- mersenne-random,
- mwc-random,
- random
diff --git a/changelog.md b/changelog.md
index 03a9525..04687fe 100644
--- a/changelog.md
+++ b/changelog.md
@@ -1,8 +1,19 @@
+## Changes in 0.14.0.0
+
+ * Low level functions for acquiring random data for initialization
+ of PRGN state is moved to `System.Random.MWC.SeedSource` module
+
+ * Ensure that carry is always correct when restoring PRNG state from
+ seed. Only affects users who create 258 element seed manually.
+ (#63, #65)
+
+
## Changes in 0.13.6.0
* `tablePoisson` now can handle λ>1923, see #59 for details.
That required intoduction of dependency on math-functions.
+
## Changes in 0.13.5.0
* `logCategorical` added
diff --git a/mwc-random.cabal b/mwc-random.cabal
index 28ba2d1..bab6518 100644
--- a/mwc-random.cabal
+++ b/mwc-random.cabal
@@ -1,5 +1,5 @@
name: mwc-random
-version: 0.13.6.0
+version: 0.14.0.0
synopsis: Fast, high quality pseudo random number generation
description:
This package contains code for generating high quality random
@@ -28,56 +28,20 @@ cabal-version: >= 1.8.0.4
extra-source-files:
changelog.md
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
- System.Random.MWC.CondensedTable
- build-depends:
- base < 5,
- primitive,
- time,
- vector >= 0.7,
- math-functions >= 0.2.1.0
- if impl(ghc >= 6.10)
- build-depends:
- base >= 4
+ exposed-modules: System.Random.MWC
+ System.Random.MWC.Distributions
+ System.Random.MWC.CondensedTable
+ System.Random.MWC.SeedSource
+ build-depends: base >= 4.5 && < 5
+ , primitive >= 0.6
+ , time
+ , vector >= 0.7
+ , math-functions >= 0.2.1.0
+
+ ghc-options: -Wall -funbox-strict-fields -fwarn-tabs
- -- gather extensive profiling data for now
- ghc-prof-options: -auto-all
-
- 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
-
- ghc-options:
- -Wall -threaded -rtsopts
-
- build-depends:
- vector >= 0.7,
- HUnit,
- QuickCheck,
- base,
- mwc-random,
- statistics >= 0.10.1.0,
- test-framework,
- test-framework-hunit,
- test-framework-quickcheck2
source-repository head
type: git
diff --git a/test/KS.hs b/test/KS.hs
deleted file mode 100644
index c91c300..0000000
--- a/test/KS.hs
+++ /dev/null
@@ -1,54 +0,0 @@
--- Kolmogorov-Smirnov tests for distribution
---
--- Note that it's not most powerful test for normality.
-module KS (
- tests
- ) where
-
-import qualified Data.Vector.Unboxed as U
-
-import Statistics.Test.KolmogorovSmirnov
-
-import Statistics.Distribution
-import Statistics.Distribution.Binomial
-import Statistics.Distribution.Exponential
-import Statistics.Distribution.Gamma
-import Statistics.Distribution.Normal
-import Statistics.Distribution.Uniform
-import Statistics.Distribution.Beta
-
-import qualified System.Random.MWC as MWC
-import qualified System.Random.MWC.Distributions as MWC
-
-import Test.HUnit hiding (Test)
-import Test.Framework
-import Test.Framework.Providers.HUnit
-
-
-tests :: MWC.GenIO -> Test
-tests g = testGroup "Kolmogorov-Smirnov"
- [ testCase "standard" $ testKS standard MWC.standard g
- , testCase "normal m=1 s=2" $ testKS (normalDistr 1 2) (MWC.normal 1 2) g
- -- Gamma distribution
- , testCase "gamma k=1 θ=1" $ testKS (gammaDistr 1 1 ) (MWC.gamma 1 1 ) g
- , testCase "gamma k=0.3 θ=0.4" $ testKS (gammaDistr 0.3 0.4) (MWC.gamma 0.3 0.4) g
- , testCase "gamma k=0.3 θ=3" $ testKS (gammaDistr 0.3 3 ) (MWC.gamma 0.3 3 ) g
- , testCase "gamma k=3 θ=0.4" $ testKS (gammaDistr 3 0.4) (MWC.gamma 3 0.4) g
- , testCase "gamma k=3 θ=3" $ testKS (gammaDistr 3 3 ) (MWC.gamma 3 3 ) g
- -- Uniform
- , testCase "uniform -2 .. 3" $ testKS (uniformDistr (-2) 3) (MWC.uniformR (-2,3)) g
- -- Exponential
- , testCase "exponential l=1" $ testKS (exponential 1) (MWC.exponential 1) g
- , testCase "exponential l=3" $ testKS (exponential 3) (MWC.exponential 3) g
- -- Beta
- , testCase "beta a=0.3,b=0.5" $ testKS (betaDistr 0.3 0.5) (MWC.beta 0.3 0.5) g
- , testCase "beta a=0.1,b=0.8" $ testKS (betaDistr 0.3 0.5) (MWC.beta 0.3 0.5) g
- , testCase "beta a=0.8,b=0.1" $ testKS (betaDistr 0.3 0.5) (MWC.beta 0.3 0.5) g
- ]
-
-testKS :: (Distribution d) => d -> (MWC.GenIO -> IO Double) -> MWC.GenIO -> IO ()
-testKS distr generator g = do
- sample <- U.replicateM 1000 (generator g)
- case kolmogorovSmirnovTest distr 0.01 sample of
- Significant -> assertFailure "KS test failed"
- NotSignificant -> return ()
diff --git a/test/QC.hs b/test/QC.hs
deleted file mode 100644
index 39f743a..0000000
--- a/test/QC.hs
+++ /dev/null
@@ -1,56 +0,0 @@
--- QC tests for random number generators
---
--- Require QuickCheck >= 2.2
-module QC (
- tests
- ) where
-
-import Control.Applicative
-
-import Data.Word (Word8,Word16,Word32,Word64,Word)
-import Data.Int (Int8, Int16, Int32, Int64 )
-
-import Test.QuickCheck
-import Test.QuickCheck.Monadic
-import Test.Framework
-import Test.Framework.Providers.QuickCheck2
-
-import System.Random.MWC
-
-
-
-----------------------------------------------------------------
-
-tests :: GenIO -> Test
-tests g = testGroup "Range"
- [ testProperty "Int8" $ (prop_InRange g :: InRange Int8)
- , testProperty "Int16" $ (prop_InRange g :: InRange Int16)
- , testProperty "Int32" $ (prop_InRange g :: InRange Int32)
- , testProperty "Int64" $ (prop_InRange g :: InRange Int64)
- , testProperty "Word8" $ (prop_InRange g :: InRange Word8)
- , testProperty "Word16" $ (prop_InRange g :: InRange Word16)
- , testProperty "Word32" $ (prop_InRange g :: InRange Word32)
- , testProperty "Word64" $ (prop_InRange g :: InRange Word64)
- , testProperty "Int" $ (prop_InRange g :: InRange Int)
- , testProperty "Word64" $ (prop_InRange g :: InRange Word)
- , testProperty "Float" $ (prop_InRange g :: InRange Float)
- , testProperty "Double" $ (prop_InRange g :: InRange Double)
- ]
-
-
-
-----------------------------------------------------------------
-
--- Test that values generated with uniformR never lie outside range.
-prop_InRange :: (Variate a, Ord a,Num a) => GenIO -> OrderedPair a -> Property
-prop_InRange g (OrderedPair (x1,x2)) = monadicIO $ do
- r <- run $ uniformR (x1,x2) g
- assert (x1 <= r && r <= x2)
-
-type InRange a = OrderedPair a -> Property
-
--- Ordered pair (x,y) for which x <= y
-newtype OrderedPair a = OrderedPair (a,a)
- deriving Show
-instance (Ord a, Arbitrary a) => Arbitrary (OrderedPair a) where
- arbitrary = OrderedPair <$> suchThat arbitrary (uncurry (<=))
diff --git a/test/run-dieharder-test.sh b/test/run-dieharder-test.sh
deleted file mode 100644
index fa39968..0000000
--- a/test/run-dieharder-test.sh
+++ /dev/null
@@ -1,26 +0,0 @@
-#!/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
diff --git a/test/tests.hs b/test/tests.hs
deleted file mode 100644
index be84fe2..0000000
--- a/test/tests.hs
+++ /dev/null
@@ -1,16 +0,0 @@
-import Test.Framework (defaultMain)
-import System.Random.MWC (withSystemRandom)
-
-import qualified QC
-import qualified ChiSquare
-import qualified KS
-
-
-main :: IO ()
-main =
- withSystemRandom $ \g ->
- defaultMain
- [ QC.tests g
- , ChiSquare.tests g
- , KS.tests g
- ]
diff --git a/test/visual.R b/test/visual.R
deleted file mode 100644
index 7b5a6a0..0000000
--- a/test/visual.R
+++ /dev/null
@@ -1,105 +0,0 @@
-# Ugly script for displaying distributions alogside with theoretical
-# distribution.
-
-
-view.dumps <- function() {
- # Load random data from dist
- load.d <- function(name) read.table(name)[,1]
- # Plots for continous distribution
- 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)
- }
- # plots for discrete distribution
- plot.ds <- function( name, xs, prob) {
- smp <- load.d( name )
- h <- hist( smp,
- breaks = c( max(xs) + 0.5, xs - 0.5),
- freq=FALSE, main = name
- )
- dh <- sqrt( h$count ) / max( 1, sum( h$count ) )
- arrows( xs, h$density + dh,
- xs, h$density - dh,
- angle=90, code=3, length=0.2 )
- points( xs, prob(xs), pch='0', col='red', type='b')
- }
- ################################################################
- # 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()
- ################################################################
- # Poisson
- plot.ds( "distr/poisson-0.1", 0:6, function(x) dpois(x, lambda=0.1) )
- readline()
- #
- plot.ds( "distr/poisson-1.0", 0:10, function(x) dpois(x, lambda=1.0) )
- readline()
- #
- plot.ds( "distr/poisson-4.5", 0:20, function(x) dpois(x, lambda=4.5) )
- readline()
- #
- plot.ds( "distr/poisson-30", 0:100, function(x) dpois(x, lambda=30) )
- readline()
- #
- ################################################################
- # Binomial
- plot.ds( "distr/binom-4-0.5", 0:4, function(x) dbinom(x, 4, 0.5) )
- readline()
- #
- plot.ds( "distr/binom-10-0.1", 0:10, function(x) dbinom(x, 10, 0.1) )
- readline()
- #
- plot.ds( "distr/binom-10-0.6", 0:10, function(x) dbinom(x, 10, 0.6) )
- readline()
- #
- plot.ds( "distr/binom-10-0.8", 0:10, function(x) dbinom(x, 10, 0.8) )
- readline()
- #
-}
diff --git a/test/visual.hs b/test/visual.hs
deleted file mode 100644
index 41329bc..0000000
--- a/test/visual.hs
+++ /dev/null
@@ -1,44 +0,0 @@
--- 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
-import qualified System.Random.MWC.CondensedTable 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 = 30000
- 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
- -- Poisson
- dumpSample n "poisson-0.1" $ MWC.genFromTable (MWC.tablePoisson 0.1) g
- dumpSample n "poisson-1.0" $ MWC.genFromTable (MWC.tablePoisson 1.0) g
- dumpSample n "poisson-4.5" $ MWC.genFromTable (MWC.tablePoisson 4.5) g
- dumpSample n "poisson-30" $ MWC.genFromTable (MWC.tablePoisson 30) g
- -- Binomial
- dumpSample n "binom-4-0.5" $ MWC.genFromTable (MWC.tableBinomial 4 0.5) g
- dumpSample n "binom-10-0.1" $ MWC.genFromTable (MWC.tableBinomial 10 0.1) g
- dumpSample n "binom-10-0.6" $ MWC.genFromTable (MWC.tableBinomial 10 0.6) g
- dumpSample n "binom-10-0.8" $ MWC.genFromTable (MWC.tableBinomial 10 0.8) g