diff --git a/System/Random/MWC.hs b/System/Random/MWC.hs --- a/System/Random/MWC.hs +++ b/System/Random/MWC.hs @@ -507,14 +507,32 @@ type instance Unsigned Int16 = Word16 type instance Unsigned Int32 = Word32 type instance Unsigned Int64 = Word64 -type instance Unsigned Int = Word type instance Unsigned Word8 = Word8 type instance Unsigned Word16 = Word16 type instance Unsigned Word32 = Word32 type instance Unsigned Word64 = Word64 -type instance Unsigned Word = Word +-- This is workaround for bug #25. +-- +-- GHC-7.6 has a bug (#8072) which results in calculation of wrong +-- number of buckets in function `uniformRange'. Consequently uniformR +-- generates values in wrong range. +-- +-- Bug only affects 32-bit systems and Int/Word data types. Word32 +-- works just fine. So we set Word32 as unsigned counterpart for Int +-- and Word on 32-bit systems. It's done only for GHC-7.6 because +-- other versions are unaffected by the bug and we expect that GHC may +-- optimise code which uses Word better. +#if (WORD_SIZE_IN_BITS < 64) && (__GLASGOW_HASKELL__ == 706) +type instance Unsigned Int = Word32 +type instance Unsigned Word = Word32 +#else +type instance Unsigned Int = Word +type instance Unsigned Word = Word +#endif + + -- Subtract two numbers under assumption that x>=y and store result in -- unsigned data type of same size sub :: (Integral a, Integral (Unsigned a)) => a -> a -> Unsigned a @@ -613,10 +631,11 @@ -- /Communications of the ACM/ 46(5):90–93. -- -- --- * Thomas, D.B.; Leong, P.G.W.; Luk, W.; Villasenor, J.D. --- (2007). Gaussian random number generators. --- /ACM Computing Surveys/ 39(4). --- +-- * Doornik, J.A. (2007) Conversion of high-period random numbers to +-- floating point. +-- /ACM Transactions on Modeling and Computer Simulation/ 17(1). +-- + -- $typehelp -- diff --git a/System/Random/MWC/CondensedTable.hs b/System/Random/MWC/CondensedTable.hs --- a/System/Random/MWC/CondensedTable.hs +++ b/System/Random/MWC/CondensedTable.hs @@ -24,6 +24,8 @@ -- ** Disrete distributions , tablePoisson , tableBinomial + -- * References + -- $references ) where import Control.Arrow (second,(***)) @@ -102,7 +104,7 @@ -- the case, this algorithm will construct a table for some -- distribution that may bear no resemblance to what you intended. tableFromProbabilities - :: (Vector v (a,Word32), Vector v (a,Double), Vector v a, Vector v Word32, Show a) + :: (Vector v (a,Word32), Vector v (a,Double), Vector v a, Vector v Word32) => v (a, Double) -> CondensedTable v a {-# INLINE tableFromProbabilities #-} tableFromProbabilities v @@ -124,7 +126,7 @@ -- probilities. Non-positive weights are discarded, and those -- remaining are normalized to 1. tableFromWeights - :: (Vector v (a,Word32), Vector v (a,Double), Vector v a, Vector v Word32, Show a) + :: (Vector v (a,Word32), Vector v (a,Double), Vector v a, Vector v Word32) => v (a, Double) -> CondensedTable v a {-# INLINE tableFromWeights #-} tableFromWeights = tableFromProbabilities . normalize . G.filter ((> 0) . snd) 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 @@ -10,14 +10,17 @@ -- -- Pseudo-random number generation for non-uniform distributions. -module System.Random.MWC.Distributions +module System.Random.MWC.Distributions ( -- * Variates: non-uniformly distributed values normal , standard , exponential + , truncatedExp , gamma , chiSquare + , geometric0 + , geometric1 -- * References -- $references @@ -106,6 +109,21 @@ return $! - log x / beta +-- | Generate truncated exponentially distributed random variate. +truncatedExp :: PrimMonad m + => Double -- ^ Scale parameter + -> (Double,Double) -- ^ Range to which distribution is + -- truncated. Values may be negative. + -> Gen (PrimState m) -- ^ Generator. + -> m Double +{-# INLINE truncatedExp #-} +truncatedExp beta (a,b) gen = do + -- We shift a to 0 and then generate distribution truncated to [0,b-a] + -- It's easier + let delta = b - a + p <- uniform gen + return $! a - log ( (1 - p) + p*exp(-beta*delta)) / beta + -- | Random variate generator for gamma distribution. gamma :: PrimMonad m => Double -- ^ Shape parameter @@ -150,7 +168,32 @@ | otherwise = do x <- gamma (0.5 * fromIntegral n) 1 gen return $! 2 * x +-- | Random variate generator for the geometric distribution, +-- computing the number of failures before success. Supports [0..]. +geometric0 :: PrimMonad m + => Double -- ^ /p/ success probability lies in (0,1] + -> Gen (PrimState m) -- ^ Generator + -> m Int +{-# INLINE geometric0 #-} +geometric0 p gen + | p == 1 = return 0 + | p > 0 && p < 1 = do q <- uniform gen + -- FIXME: We want to use log1p here but it will + -- introduce dependency on math-functions. + return $! floor $ log q / log (1 - p) + | otherwise = pkgError "geometric0" "probability out of [0,1] range" +-- | Random variate generator for geometric distribution for number of +-- trials. Supports [1..] (i.e. just 'geometric0' shifted by 1). +geometric1 :: PrimMonad m + => Double -- ^ /p/ success probability lies in (0,1] + -> Gen (PrimState m) -- ^ Generator + -> m Int +{-# INLINE geometric1 #-} +geometric1 p gen = do n <- geometric0 p gen + return $! n + 1 + + sqr :: Double -> Double sqr x = x * x {-# INLINE sqr #-} @@ -165,7 +208,7 @@ -- normal random samples. Mimeo, Nuffield College, University of -- Oxford. -- --- * Doornik, J.A. (2007) Conversion of high-period random numbers to --- floating point. --- /ACM Transactions on Modeling and Computer Simulation/ 17(1). --- +-- * Thomas, D.B.; Leong, P.G.W.; Luk, W.; Villasenor, J.D. +-- (2007). Gaussian random number generators. +-- /ACM Computing Surveys/ 39(4). +-- 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.12.0.1 +version: 0.13.0.0 synopsis: Fast, high quality pseudo random number generation description: This package contains code for generating high quality random