diff --git a/System/Random/MWC.hs b/System/Random/MWC.hs --- a/System/Random/MWC.hs +++ b/System/Random/MWC.hs @@ -1,5 +1,11 @@ -{-# LANGUAGE BangPatterns, CPP, DeriveDataTypeable, MagicHash, Rank2Types, - ScopedTypeVariables #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE MagicHash #-} +{-# LANGUAGE Rank2Types #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} -- | -- Module : System.Random.MWC -- Copyright : (c) 2009, 2010 Bryan O'Sullivan @@ -44,31 +50,36 @@ #include "MachDeps.h" #endif -import Control.Exception (IOException, catch) -import Control.Monad (ap, liftM, unless) +import Control.Exception (IOException, catch) +import Control.Monad (ap, liftM, unless) import Control.Monad.Primitive (PrimMonad, PrimState, unsafePrimToIO) -import Control.Monad.ST (ST) -import Data.Bits ((.&.), (.|.), 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, unsafeFreeze) -import Data.Word (Word, Word8, Word16, Word32, Word64) -import Foreign.Marshal.Alloc (allocaBytes) -import Foreign.Marshal.Array (peekArray) -import GHC.Base (Int(I#)) -import GHC.Word (Word64(W64#), uncheckedShiftL64#, uncheckedShiftRL64#) +import Control.Monad.ST (ST) +import Data.Bits ((.&.), (.|.), 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, unsafeFreeze) +import Data.Word (Word, Word8, Word16, Word32, Word64) +import Foreign.Marshal.Alloc (allocaBytes) +import Foreign.Marshal.Array (peekArray) +import GHC.Base (Int(I#)) +import GHC.Word (Word64(W64#), uncheckedShiftL64#, uncheckedShiftRL64#) import Prelude hiding (catch) -import qualified Data.Vector.Generic as G +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 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.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. -- @@ -98,50 +109,83 @@ -- 2**(-33). To do the same with 'Double' variates, subtract -- 2**(-53). uniform :: (PrimMonad m) => Gen (PrimState m) -> m a + -- | Generate single uniformly distributed random variable in a + -- given range. + -- + -- * For integral types inclusive range is used. + -- + -- * For floating point numbers range (a,b] is used if one ignores + -- rounding errors. + uniformR :: (PrimMonad m) => (a,a) -> Gen (PrimState m) -> m a instance Variate Int8 where - uniform = uniform1 fromIntegral - {-# INLINE uniform #-} + uniform = uniform1 fromIntegral + uniformR = uniformRange + {-# INLINE uniform #-} + {-# INLINE uniformR #-} instance Variate Int16 where - uniform = uniform1 fromIntegral - {-# INLINE uniform #-} + uniform = uniform1 fromIntegral + uniformR = uniformRange + {-# INLINE uniform #-} + {-# INLINE uniformR #-} instance Variate Int32 where - uniform = uniform1 fromIntegral - {-# INLINE uniform #-} + uniform = uniform1 fromIntegral + uniformR = uniformRange + {-# INLINE uniform #-} + {-# INLINE uniformR #-} instance Variate Int64 where - uniform = uniform2 wordsTo64Bit - {-# INLINE uniform #-} + uniform = uniform2 wordsTo64Bit + uniformR = uniformRange + {-# INLINE uniform #-} + {-# INLINE uniformR #-} instance Variate Word8 where - uniform = uniform1 fromIntegral - {-# INLINE uniform #-} + uniform = uniform1 fromIntegral + uniformR = uniformRange + {-# INLINE uniform #-} + {-# INLINE uniformR #-} instance Variate Word16 where - uniform = uniform1 fromIntegral - {-# INLINE uniform #-} + uniform = uniform1 fromIntegral + uniformR = uniformRange + {-# INLINE uniform #-} + {-# INLINE uniformR #-} instance Variate Word32 where - uniform = uniformWord32 - {-# INLINE uniform #-} + uniform = uniform1 fromIntegral + uniformR = uniformRange + {-# INLINE uniform #-} + {-# INLINE uniformR #-} instance Variate Word64 where - uniform = uniform2 wordsTo64Bit - {-# INLINE uniform #-} + uniform = uniform2 wordsTo64Bit + uniformR = uniformRange + {-# INLINE uniform #-} + {-# INLINE uniformR #-} instance Variate Bool where uniform = uniform1 wordToBool - {-# INLINE uniform #-} + uniformR (False,True) g = uniform g + uniformR (False,False) _ = return False + uniformR (True,True) _ = return True + uniformR (True,False) g = uniform g + {-# INLINE uniform #-} + {-# INLINE uniformR #-} instance Variate Float where - uniform = uniform1 wordToFloat - {-# INLINE uniform #-} + uniform = uniform1 wordToFloat + uniformR (x1,x2) = uniform1 (\w -> x1 + (x2-x1) * wordToFloat w) + {-# INLINE uniform #-} + {-# INLINE uniformR #-} instance Variate Double where - uniform = uniform2 wordsToDouble - {-# INLINE uniform #-} + uniform = uniform2 wordsToDouble + uniformR (x1,x2) = uniform2 (\w1 w2 -> x1 + (x2-x1) * wordsToDouble w1 w2) + {-# INLINE uniform #-} + {-# INLINE uniformR #-} instance Variate Int where #if WORD_SIZE_IN_BITS < 64 @@ -149,7 +193,9 @@ #else uniform = uniform2 wordsTo64Bit #endif - {-# INLINE uniform #-} + uniformR = uniformRange + {-# INLINE uniform #-} + {-# INLINE uniformR #-} instance Variate Word where #if WORD_SIZE_IN_BITS < 64 @@ -157,7 +203,9 @@ #else uniform = uniform2 wordsTo64Bit #endif - {-# INLINE uniform #-} + uniformR = uniformRange + {-# INLINE uniform #-} + {-# INLINE uniformR #-} {- instance Variate Integer where @@ -169,16 +217,25 @@ instance (Variate a, Variate b) => Variate (a,b) where uniform g = (,) `liftM` uniform g `ap` uniform g - {-# INLINE uniform #-} + uniformR ((x1,y1),(x2,y2)) g = (,) `liftM` uniformR (x1,x2) g `ap` uniformR (y1,y2) g + {-# INLINE uniform #-} + {-# INLINE uniformR #-} instance (Variate a, Variate b, Variate c) => Variate (a,b,c) where uniform g = (,,) `liftM` uniform g `ap` uniform g `ap` uniform g - {-# INLINE uniform #-} + uniformR ((x1,y1,z1),(x2,y2,z2)) g = + (,,) `liftM` uniformR (x1,x2) g `ap` uniformR (y1,y2) g `ap` uniformR (z1,z2) g + {-# INLINE uniform #-} + {-# INLINE uniformR #-} instance (Variate a, Variate b, Variate c, Variate d) => Variate (a,b,c,d) where uniform g = (,,,) `liftM` uniform g `ap` uniform g `ap` uniform g `ap` uniform g - {-# INLINE uniform #-} + uniformR ((x1,y1,z1,t1),(x2,y2,z2,t2)) g = + (,,,) `liftM` uniformR (x1,x2) g `ap` uniformR (y1,y2) g `ap` + uniformR (z1,z2) g `ap` uniformR (t1,t2) g + {-# INLINE uniform #-} + {-# INLINE uniformR #-} wordsTo64Bit :: Integral a => Word32 -> Word32 -> a wordsTo64Bit a b = @@ -261,6 +318,9 @@ deriving (Eq, Show, Typeable) -- 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) @@ -336,8 +396,8 @@ uniformWord32 :: PrimMonad m => Gen (PrimState m) -> m Word32 uniformWord32 (Gen q) = do let a = 809430660 :: Word64 - i <- nextIndex `liftM` M.unsafeRead q ioff - c <- fromIntegral `liftM` M.unsafeRead q coff + 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 t32 = fromIntegral t @@ -356,9 +416,9 @@ uniform2 :: PrimMonad m => (Word32 -> Word32 -> a) -> Gen (PrimState m) -> m a uniform2 f (Gen q) = do let a = 809430660 :: Word64 - i <- nextIndex `liftM` M.unsafeRead q ioff + i <- nextIndex `liftM` M.unsafeRead q ioff let j = nextIndex i - c <- fromIntegral `liftM` M.unsafeRead q coff + 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 @@ -372,6 +432,54 @@ M.unsafeWrite q coff (fromIntegral (u `shiftR` 32)) return $! f t32 u32 {-# INLINE uniform2 #-} + +-- Type family for fixed size integrals. For signed data types it's +-- its unsigned couterpart with same size and for unsigned data types +-- it's same type +type family Unsigned a :: * + +type instance Unsigned Int8 = Word8 +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 + +-- 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 +sub x y = fromIntegral x - fromIntegral y + +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 #-} + +-- Generate unformly distributed value in inclusive range. +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 +{-# INLINE uniformRange #-} -- | Generate a vector of pseudo-random variates. This is not -- necessarily faster than invoking 'uniform' repeatedly in a loop, diff --git a/benchmarks/Benchmark.hs b/benchmarks/Benchmark.hs --- a/benchmarks/Benchmark.hs +++ b/benchmarks/Benchmark.hs @@ -1,10 +1,25 @@ -{-# LANGUAGE Rank2Types #-} import Control.Monad.ST + +import Data.Int +import Data.Word + import Criterion.Main import System.Random.MWC main = do gen <- create - defaultMain [ - bench "mwc" (uniform gen :: IO Double) - ] + defaultMain + [ bench "mwc-Double" (uniform gen :: IO Double) + , bench "mwc-Int" (uniform gen :: IO Int ) + , bench "mwc-Int8" (uniform gen :: IO Int8 ) + , bench "mwc-Int16" (uniform gen :: IO Int16 ) + , bench "mwc-Int32" (uniform gen :: IO Int32 ) + , bench "mwc-Int64" (uniform gen :: IO Int64 ) + , bench "mwc-Word" (uniform gen :: IO Word ) + , bench "mwc-Word8" (uniform gen :: IO Word8 ) + , bench "mwc-Word16" (uniform gen :: IO Word16) + , bench "mwc-Word32" (uniform gen :: IO Word32) + , bench "mwc-Word64" (uniform gen :: IO Word64) + , bench "mwc-Integer" (uniform gen :: IO Word64) + , bench "normal" (normal gen :: IO Double) + ] 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.7.0.1 +version: 0.8.0.0 synopsis: Fast, high quality pseudo random number generation description: This package contains code for generating high quality random @@ -22,7 +22,7 @@ copyright: 2009, 2010 Bryan O'Sullivan category: Math, Statistics build-type: Simple -cabal-version: >= 1.2 +cabal-version: >= 1.6 extra-source-files: README benchmarks/Benchmark.hs @@ -46,3 +46,7 @@ ghc-options: -Wall -funbox-strict-fields if impl(ghc >= 6.8) ghc-options: -fwarn-tabs + +source-repository head + type: mercurial + location: http://bitbucket.org/bos/mwc-random