random-fu 0.1.0.0 → 0.1.3

raw patch · 26 files changed

+525/−419 lines, 26 filesdep ~basedep ~mtlPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependency ranges changed: base, mtl

API changes (from Hackage documentation)

Files

• random-fu.cabal +22/−14
• src/Data/Random.hs +4/−4
• src/Data/Random/Distribution/Bernoulli.hs +15/−9
• src/Data/Random/Distribution/Beta.hs +12/−7
• src/Data/Random/Distribution/Binomial.hs +24/−15
• src/Data/Random/Distribution/Categorical.hs +10/−5
• src/Data/Random/Distribution/Dirichlet.hs +6/−3
• src/Data/Random/Distribution/Exponential.hs +6/−3
• src/Data/Random/Distribution/Gamma.hs +24/−10
• src/Data/Random/Distribution/Multinomial.hs +5/−2
• src/Data/Random/Distribution/Normal.hs +27/−20
• src/Data/Random/Distribution/Poisson.hs +14/−11
• src/Data/Random/Distribution/Rayleigh.hs +6/−3
• src/Data/Random/Distribution/Triangular.hs +3/−3
• src/Data/Random/Distribution/Uniform.hs +82/−64
• src/Data/Random/Distribution/Ziggurat.hs +25/−23
• src/Data/Random/Internal/Fixed.hs +8/−8
• src/Data/Random/Internal/Primitives.hs +16/−2
• src/Data/Random/Internal/Words.hs +2/−2
• src/Data/Random/RVar.hs +7/−6
• src/Data/Random/Source.hs +58/−99
• src/Data/Random/Source/DevRandom.hs +24/−15
• src/Data/Random/Source/MWC.hs +19/−24
• src/Data/Random/Source/PureMT.hs +85/−42
• src/Data/Random/Source/Std.hs +0/−7
• src/Data/Random/Source/StdGen.hs +21/−18

@@ -1,8 +1,8 @@ name:                   random-fu-version:                0.1.0.0+version:                0.1.3 stability:              provisional -cabal-version:          >= 1.2+cabal-version:          >= 1.6 build-type:             Simple  author:                 James Cook <james.cook@usma.edu>@@ -38,16 +38,23 @@                         eventually so that client code dependencies on it will                          be made explicit).                         +                        Support for "base" packages earlier than version 4+                        (and thus GHC releases earlier than 6.10) has been +                        dropped, as too many of this package's dependencies do+                        not support older versions.+                                                 The "Data.Random" module itself should now have a                         relatively stable interface, but the other modules                         are still subject to change.  Specifically, I am                          considering hiding data constructors for most or all                          of the distributions. -Flag base4 Flag base4_2     Description:        base-4.2 has an incompatible change in Data.Fixed (HasResolution) +Flag mtl2+    Description:        mtl-2 has State, etc., as "type" rather than "newtype"+ Library   ghc-options:          -Wall -funbox-strict-fields -fno-method-sharing   hs-source-dirs:       src@@ -82,27 +89,28 @@                         Data.Random.Source.PureMT                         Data.Random.Source.Std                         Data.Random.Source.StdGen-  if flag(base4)-    build-depends:      syb-    -    if flag(base4_2)-      build-depends:    base >= 4 && <4.2-    else-      cpp-options:      -Dbase_4_2-      build-depends:    base >= 4.2 && <5+  if flag(base4_2)+    build-depends:      base >= 4.2 && <5   else-    build-depends:      base >= 3 && < 4-    +    cpp-options:        -Dold_Fixed+    build-depends:      base >= 4 && <4.2+  +  if flag(mtl2)+    build-depends:      mtl == 2.*+    cpp-options:        -DMTL2+  else+    build-depends:      mtl == 1.*+     build-depends:        array,                         containers,                         mersenne-random-pure64,                         monad-loops >= 0.3.0.1,                         MonadPrompt,                         mwc-random,-                        mtl,                         random,                         random-shuffle,                         stateref >= 0.3 && < 0.4,+                        syb,                         tagged,                         template-haskell,                         vector

@@ -48,10 +48,10 @@       Sampleable(..), sample, sampleState, sampleStateT,              -- * A few very common distributions-      Uniform(..), uniform, -      StdUniform(..), stdUniform,-      Normal(..), normal, stdNormal,-      Gamma(..), gamma,+      Uniform(..), uniform, uniformT,+      StdUniform(..), stdUniform, stdUniformT,+      Normal(..), normal, stdNormal, normalT, stdNormalT,+      Gamma(..), gamma, gammaT,              -- * Entropy Sources       MonadRandom, RandomSource, StdRandom(..),

@@ -25,11 +25,17 @@ bernoulli :: Distribution (Bernoulli b) a => b -> RVar a bernoulli p = rvar (Bernoulli p) +-- |Generate a Bernoulli process with the given probability.  For @Bool@ results,+-- @bernoulli p@ will return True (p*100)% of the time and False otherwise.+-- For numerical types, True is replaced by 1 and False by 0.+bernoulliT :: Distribution (Bernoulli b) a => b -> RVarT m a+bernoulliT p = rvarT (Bernoulli p)+ -- |A random variable whose value is 'True' the given fraction of the time -- and 'False' the rest.-boolBernoulli :: (Fractional a, Ord a, Distribution StdUniform a) => a -> RVar Bool+boolBernoulli :: (Fractional a, Ord a, Distribution StdUniform a) => a -> RVarT m Bool boolBernoulli p = do-    x <- stdUniform+    x <- stdUniformT     return (x <= p)  boolBernoulliCDF :: (Real a) => a -> Bool -> Double@@ -38,9 +44,9 @@  -- | @generalBernoulli t f p@ generates a random variable whose value is @t@ -- with probability @p@ and @f@ with probability @1-p@.-generalBernoulli :: Distribution (Bernoulli b) Bool => a -> a -> b -> RVar a+generalBernoulli :: Distribution (Bernoulli b) Bool => a -> a -> b -> RVarT m a generalBernoulli f t p = do-    x <- bernoulli p+    x <- bernoulliT p     return (if x then t else f)  generalBernoulliCDF :: CDF (Bernoulli b) Bool => (a -> a -> Bool) -> a -> a -> b -> a -> Double@@ -55,7 +61,7 @@ instance (Fractional b, Ord b, Distribution StdUniform b)         => Distribution (Bernoulli b) Bool     where-        rvar (Bernoulli p) = boolBernoulli p+        rvarT (Bernoulli p) = boolBernoulli p instance (Distribution (Bernoulli b) Bool, Real b)        => CDF (Bernoulli b) Bool     where@@ -65,7 +71,7 @@         instance Distribution (Bernoulli b) Bool                => Distribution (Bernoulli b) Int               where-                  rvar (Bernoulli p) = generalBernoulli 0 1 p+                  rvarT (Bernoulli p) = generalBernoulli 0 1 p         instance CDF (Bernoulli b) Bool               => CDF (Bernoulli b) Int               where@@ -76,7 +82,7 @@         instance Distribution (Bernoulli b) Bool                => Distribution (Bernoulli b) Float               where-                  rvar (Bernoulli p) = generalBernoulli 0 1 p+                  rvarT (Bernoulli p) = generalBernoulli 0 1 p         instance CDF (Bernoulli b) Bool               => CDF (Bernoulli b) Float               where@@ -86,7 +92,7 @@ instance (Distribution (Bernoulli b) Bool, Integral a)        => Distribution (Bernoulli b) (Ratio a)           where-           rvar (Bernoulli p) = generalBernoulli 0 1 p+           rvarT (Bernoulli p) = generalBernoulli 0 1 p instance (CDF (Bernoulli b) Bool, Integral a)        => CDF (Bernoulli b) (Ratio a)           where@@ -94,7 +100,7 @@ instance (Distribution (Bernoulli b) Bool, RealFloat a)        => Distribution (Bernoulli b) (Complex a)        where-           rvar (Bernoulli p) = generalBernoulli 0 1 p+           rvarT (Bernoulli p) = generalBernoulli 0 1 p instance (CDF (Bernoulli b) Bool, RealFloat a)        => CDF (Bernoulli b) (Complex a)        where

@@ -17,13 +17,13 @@ import Data.Random.Distribution.Gamma import Data.Random.Distribution.Uniform -{-# SPECIALIZE fractionalBeta :: Float  -> Float  -> RVar Float #-}-{-# SPECIALIZE fractionalBeta :: Double -> Double -> RVar Double #-}-fractionalBeta :: (Fractional a, Distribution Gamma a, Distribution StdUniform a) => a -> a -> RVar a-fractionalBeta 1 1 = stdUniform+{-# SPECIALIZE fractionalBeta :: Float  -> Float  -> RVarT m Float #-}+{-# SPECIALIZE fractionalBeta :: Double -> Double -> RVarT m Double #-}+fractionalBeta :: (Fractional a, Distribution Gamma a, Distribution StdUniform a) => a -> a -> RVarT m a+fractionalBeta 1 1 = stdUniformT fractionalBeta a b = do-    x <- gamma a 1-    y <- gamma b 1+    x <- gammaT a 1+    y <- gammaT b 1     return (x / (x + y))  {-# SPECIALIZE beta :: Float  -> Float  -> RVar Float #-}@@ -31,9 +31,14 @@ beta :: Distribution Beta a => a -> a -> RVar a beta a b = rvar (Beta a b) +{-# SPECIALIZE betaT :: Float  -> Float  -> RVarT m Float #-}+{-# SPECIALIZE betaT :: Double -> Double -> RVarT m Double #-}+betaT :: Distribution Beta a => a -> a -> RVarT m a+betaT a b = rvarT (Beta a b)+ data Beta a = Beta a a  $( replicateInstances ''Float realFloatTypes [d|         instance Distribution Beta Float-              where rvar (Beta a b) = fractionalBeta a b+              where rvarT (Beta a b) = fractionalBeta a b     |])

@@ -22,19 +22,16 @@     -- note that although it's fast enough for large (eg, 2^10000)      -- @Integer@s, it's not accurate enough when using @Double@ as     -- the @b@ parameter.-{-# SPECIALIZE integralBinomial :: Int -> Float  -> RVar Int #-}-{-# SPECIALIZE integralBinomial :: Int -> Double -> RVar Int #-}-{-# SPECIALIZE integralBinomial :: Integer -> Float  -> RVar Integer #-}-{-# SPECIALIZE integralBinomial :: Integer -> Double -> RVar Integer #-}-integralBinomial :: (Integral a, Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => a -> b -> RVar a+integralBinomial :: (Integral a, Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => a -> b -> RVarT m a integralBinomial = bin 0     where+        bin :: (Integral a, Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => a -> a -> b -> RVarT m a         bin !k !t !p             | t > 10    = do                 let a = 1 + t `div` 2                     b = 1 + t - a         -                x <- beta (fromIntegral a) (fromIntegral b)+                x <- betaT (fromIntegral a) (fromIntegral b)                 if x >= p                     then bin  k      (a - 1) (p / x)                     else bin (k + a) (b - 1) ((p - x) / (1 - x))@@ -43,14 +40,14 @@                 where                     count !k'  0    = return k'                     count !k' (n+1) = do-                        x <- stdUniform+                        x <- stdUniformT                         count (if x < p then k' + 1 else k') n                     count _ _ = error "integralBinomial: negative number of trials specified"  -- TODO: improve performance integralBinomialCDF :: (Integral a, Real b) => a -> b -> a -> Double integralBinomialCDF t p x = sum-    [ fromIntegral (t `c` i) * p' ^^ i * (1-p') ^^ (t-i)+    [ fromInteger (toInteger t `c` toInteger i) * p' ^^ i * (1-p') ^^ (t-i)     | i <- [0 .. x]     ]     @@ -70,17 +67,28 @@ floatingBinomialCDF :: (CDF (Binomial b) Integer, RealFrac a) => a -> b -> a -> Double floatingBinomialCDF t p x = cdf (Binomial (truncate t :: Integer) p) (floor x) -{-# SPECIALIZE binomial :: Int -> Float  -> RVar Int #-}-{-# SPECIALIZE binomial :: Int -> Double -> RVar Int #-}+{-# SPECIALIZE binomial :: Int     -> Float  -> RVar Int #-}+{-# SPECIALIZE binomial :: Int     -> Double -> RVar Int #-} {-# SPECIALIZE binomial :: Integer -> Float  -> RVar Integer #-} {-# SPECIALIZE binomial :: Integer -> Double -> RVar Integer #-}-{-# SPECIALIZE binomial :: Float  -> Float  -> RVar Float  #-}-{-# SPECIALIZE binomial :: Float  -> Double -> RVar Float  #-}-{-# SPECIALIZE binomial :: Double -> Float  -> RVar Double #-}-{-# SPECIALIZE binomial :: Double -> Double -> RVar Double #-}+{-# SPECIALIZE binomial :: Float   -> Float  -> RVar Float  #-}+{-# SPECIALIZE binomial :: Float   -> Double -> RVar Float  #-}+{-# SPECIALIZE binomial :: Double  -> Float  -> RVar Double #-}+{-# SPECIALIZE binomial :: Double  -> Double -> RVar Double #-} binomial :: Distribution (Binomial b) a => a -> b -> RVar a binomial t p = rvar (Binomial t p) +{-# SPECIALIZE binomialT :: Int     -> Float  -> RVarT m Int #-}+{-# SPECIALIZE binomialT :: Int     -> Double -> RVarT m Int #-}+{-# SPECIALIZE binomialT :: Integer -> Float  -> RVarT m Integer #-}+{-# SPECIALIZE binomialT :: Integer -> Double -> RVarT m Integer #-}+{-# SPECIALIZE binomialT :: Float   -> Float  -> RVarT m Float  #-}+{-# SPECIALIZE binomialT :: Float   -> Double -> RVarT m Float  #-}+{-# SPECIALIZE binomialT :: Double  -> Float  -> RVarT m Double #-}+{-# SPECIALIZE binomialT :: Double  -> Double -> RVarT m Double #-}+binomialT :: Distribution (Binomial b) a => a -> b -> RVarT m a+binomialT t p = rvarT (Binomial t p)+ data Binomial b a = Binomial a b  $( replicateInstances ''Int integralTypes [d|@@ -88,7 +96,8 @@                  , Distribution Beta b                  , Distribution StdUniform b                  ) => Distribution (Binomial b) Int-            where rvar (Binomial t p) = integralBinomial t p+            where +                rvarT (Binomial t p) = integralBinomial t p         instance ( Real b , Distribution (Binomial b) Int                  ) => CDF (Binomial b) Int             where cdf  (Binomial t p) = integralBinomialCDF t p

@@ -21,11 +21,16 @@ import Data.List import Data.Function --- |Construct a 'Categorical' distribution from a list of probabilities+-- |Construct a 'Categorical' random variable from a list of probabilities -- and categories, where the probabilities all sum to 1. categorical :: Distribution (Categorical p) a => [(p,a)] -> RVar a categorical ps = rvar (Categorical ps) +-- |Construct a 'Categorical' random process from a list of probabilities+-- and categories, where the probabilities all sum to 1.+categoricalT :: Distribution (Categorical p) a => [(p,a)] -> RVarT m a+categoricalT ps = rvarT (Categorical ps)+ -- | Construct a 'Categorical' distribution from a list of weighted categories, -- where the weights do not necessarily sum to 1. {-# INLINE weightedCategorical #-}@@ -51,12 +56,12 @@     deriving (Eq, Show)  instance (Fractional p, Ord p, Distribution StdUniform p) => Distribution (Categorical p) a where-    rvar (Categorical []) = fail "categorical distribution over empty set cannot be sampled"-    rvar (Categorical ds) = do+    rvarT (Categorical []) = fail "categorical distribution over empty set cannot be sampled"+    rvarT (Categorical ds) = do         let (ps, xs) = unzip ds             cs = scanl1 (+) ps         -        u <- stdUniform+        u <- stdUniformT         getEvent u cs xs                  where@@ -68,7 +73,7 @@             getEvent u cs0 xs0 = go 0 cs0 xs0                 where                     go lastC [] _-                        | lastC > 0 = do {newU <- stdUniform; getEvent newU cs0 xs0}+                        | lastC > 0 = do {newU <- stdUniformT; getEvent newU cs0 xs0}                         | otherwise = fail "categorical distribution sampling error: total probablility not greater than zero"                     go lastC (c:cs) (x:xs)                         | c < lastC = fail "categorical distribution sampling error: negative probability for an event!"

@@ -12,11 +12,11 @@  import Data.List -fractionalDirichlet :: (Fractional a, Distribution Gamma a) => [a] -> RVar [a]+fractionalDirichlet :: (Fractional a, Distribution Gamma a) => [a] -> RVarT m [a] fractionalDirichlet []  = return [] fractionalDirichlet [_] = return [1] fractionalDirichlet as = do-    xs <- sequence [gamma a 1 | a <- as]+    xs <- sequence [gammaT a 1 | a <- as]     let total = foldl1' (+) xs          return (map (* recip total) xs)@@ -24,7 +24,10 @@ dirichlet :: Distribution Dirichlet [a] => [a] -> RVar [a] dirichlet as = rvar (Dirichlet as) +dirichletT :: Distribution Dirichlet [a] => [a] -> RVarT m [a]+dirichletT as = rvarT (Dirichlet as)+ newtype Dirichlet a = Dirichlet a deriving (Eq, Show)  instance (Fractional a, Distribution Gamma a) => Distribution Dirichlet [a] where-    rvar (Dirichlet as) = fractionalDirichlet as+    rvarT (Dirichlet as) = fractionalDirichlet as

@@ -15,9 +15,9 @@  data Exponential a = Exp a -floatingExponential :: (Floating a, Distribution StdUniform a) => a -> RVar a+floatingExponential :: (Floating a, Distribution StdUniform a) => a -> RVarT m a floatingExponential lambdaRecip = do-    x <- stdUniform+    x <- stdUniformT     return (negate (log x) * lambdaRecip)  floatingExponentialCDF :: Real a => a -> a -> Double@@ -26,7 +26,10 @@ exponential :: Distribution Exponential a => a -> RVar a exponential = rvar . Exp +exponentialT :: Distribution Exponential a => a -> RVarT m a+exponentialT = rvarT . Exp+ instance (Floating a, Distribution StdUniform a) => Distribution Exponential a where-    rvar (Exp lambdaRecip) = floatingExponential lambdaRecip+    rvarT (Exp lambdaRecip) = floatingExponential lambdaRecip instance (Real a, Distribution Exponential a) => CDF Exponential a where     cdf  (Exp lambdaRecip) = floatingExponentialCDF lambdaRecip

@@ -4,7 +4,15 @@     UndecidableInstances, BangPatterns   #-} -module Data.Random.Distribution.Gamma where+module Data.Random.Distribution.Gamma+    ( Gamma(..)+    , gamma, gammaT+    +    , Erlang(..)+    , erlang, erlangT+    +    , mtGamma+    ) where  import Data.Random.RVar import Data.Random.Distribution@@ -13,18 +21,18 @@  import Data.Ratio --- derived from  Marsaglia & Tang, "A Simple Method for generating gamma+-- |derived from  Marsaglia & Tang, "A Simple Method for generating gamma -- variables", ACM Transactions on Mathematical Software, Vol 26, No 3 (2000), p363-372.-{-# SPECIALIZE mtGamma :: Double -> Double -> RVar Double #-}-{-# SPECIALIZE mtGamma :: Float -> Float -> RVar Float #-}+{-# SPECIALIZE mtGamma :: Double -> Double -> RVarT m Double #-}+{-# SPECIALIZE mtGamma :: Float  -> Float  -> RVarT m Float  #-} mtGamma     :: (Floating a, Ord a,         Distribution StdUniform a,          Distribution Normal a)-    => a -> a -> RVar a+    => a -> a -> RVarT m a mtGamma a b      | a < 1     = do-        u <- stdUniform+        u <- stdUniformT         mtGamma (1+a) $! (b * u ** recip a)     | otherwise = go     where@@ -32,13 +40,13 @@         !c = recip (sqrt (9*d))                  go = do-            x <- stdNormal+            x <- stdNormalT             let !v   = 1 + c*x                          if v <= 0                 then go                 else do-                    u  <- stdUniform+                    u  <- stdUniformT                     let !x_2 = x*x; !x_4 = x_2*x_2                         v3 = v*v*v                         dv = d * v3@@ -52,16 +60,22 @@ gamma :: (Distribution Gamma a) => a -> a -> RVar a gamma a b = rvar (Gamma a b) +gammaT :: (Distribution Gamma a) => a -> a -> RVarT m a+gammaT a b = rvarT (Gamma a b)+ erlang :: (Distribution (Erlang a) b) => a -> RVar b erlang a = rvar (Erlang a) +erlangT :: (Distribution (Erlang a) b) => a -> RVarT m b+erlangT a = rvarT (Erlang a)+ data Gamma a    = Gamma a a data Erlang a b = Erlang a  instance (Floating a, Ord a, Distribution Normal a, Distribution StdUniform a) => Distribution Gamma a where     {-# SPECIALIZE instance Distribution Gamma Double #-}     {-# SPECIALIZE instance Distribution Gamma Float #-}-    rvar (Gamma a b) = mtGamma a b+    rvarT (Gamma a b) = mtGamma a b  instance (Integral a, Floating b, Ord b, Distribution Normal b, Distribution StdUniform b) => Distribution (Erlang a) b where-    rvar (Erlang a) = mtGamma (fromIntegral a) 1+    rvarT (Erlang a) = mtGamma (fromIntegral a) 1

@@ -8,18 +8,21 @@ multinomial :: Distribution (Multinomial p) [a] => [p] -> a -> RVar [a] multinomial ps n = rvar (Multinomial ps n) +multinomialT :: Distribution (Multinomial p) [a] => [p] -> a -> RVarT m [a]+multinomialT ps n = rvarT (Multinomial ps n)+ data Multinomial p a where     Multinomial :: [p] -> a -> Multinomial p [a]  instance (Num a, Fractional p, Distribution (Binomial p) a) => Distribution (Multinomial p) [a] where     -- TODO: implement faster version based on Categorical for small n, large (length ps)-    rvar (Multinomial ps0 t) = go t ps0 (tailSums ps0) id+    rvarT (Multinomial ps0 t) = go t ps0 (tailSums ps0) id         where             go _ []     _            f = return (f [])             go n [_]    _            f = return (f [n])             go 0 (_:ps) (_   :psums) f = go 0 ps psums (f . (0:))             go n (p:ps) (psum:psums) f = do-                x <- binomial n (p / psum)+                x <- binomialT n (p / psum)                 go (n-x) ps psums (f . (x:))                          go _ _ _ _ = error "rvar/Multinomial: programming error! this case should be impossible!"

@@ -1,12 +1,12 @@ {-# LANGUAGE     MultiParamTypeClasses, FlexibleInstances, FlexibleContexts,-    UndecidableInstances, ForeignFunctionInterface, BangPatterns+    UndecidableInstances, ForeignFunctionInterface, BangPatterns, +    RankNTypes   #-}- module Data.Random.Distribution.Normal     ( Normal(..)-    , normal-    , stdNormal+    , normal, normalT+    , stdNormal, stdNormalT          , doubleStdNormal     , floatStdNormal@@ -78,13 +78,13 @@ -- |Draw from the tail of a normal distribution (the region beyond the provided value) {-# INLINE normalTail #-} normalTail :: (Distribution StdUniform a, Floating a, Ord a) =>-              a -> RVar a+              a -> RVarT m a normalTail r = go     where         go = do-            !u <- stdUniform+            !u <- stdUniformT             let !x = log u / r-            !v <- stdUniform+            !v <- stdUniformT             let !y = log v             if x*x + y+y > 0                 then go@@ -94,7 +94,7 @@ -- @logBase 2 c@ and the 'zGetIU' implementation. normalZ ::   (RealFloat a, Erf a, Vector v a, Distribution Uniform a, Integral b) =>-  b -> RVar (Int, a) -> Ziggurat v a+  b -> (forall m. RVarT m (Int, a)) -> Ziggurat v a normalZ p = mkZigguratRec True normalF normalFInv normalFInt normalFVol (2^p)  -- | Ziggurat target function (upper half of a non-normalized gaussian PDF)@@ -130,21 +130,21 @@ -- -- As far as I know, this should be safe to use in any monomorphic -- @Distribution Normal@ instance declaration.-realFloatStdNormal :: (RealFloat a, Erf a, Distribution Uniform a) => RVar a+realFloatStdNormal :: (RealFloat a, Erf a, Distribution Uniform a) => RVarT m a realFloatStdNormal = runZiggurat (normalZ p getIU `asTypeOf` (undefined :: Ziggurat V.Vector a))     where          p :: Int         p = 6         +        getIU :: (Num a, Distribution Uniform a) => RVarT m (Int, a)         getIU = do             i <- getRandomPrim PrimWord8-            u <- uniform (-1) 1+            u <- uniformT (-1) 1             return (fromIntegral i .&. (2^p-1), u)  -- |A random variable sampling from the standard normal distribution -- over the 'Double' type.-{-# NOINLINE doubleStdNormal #-}-doubleStdNormal :: RVar Double+doubleStdNormal :: RVarT m Double doubleStdNormal = runZiggurat doubleStdNormalZ  -- doubleStdNormalC must not be over 2^12 if using wordToDoubleWithExcess@@ -162,6 +162,7 @@         getIU         (normalTail doubleStdNormalR)     where +        getIU :: RVarT m (Int, Double)         getIU = do             !w <- getRandomPrim PrimWord64             let (u,i) = wordToDoubleWithExcess w@@ -169,8 +170,7 @@  -- |A random variable sampling from the standard normal distribution -- over the 'Float' type.-{-# NOINLINE floatStdNormal #-}-floatStdNormal :: RVar Float+floatStdNormal :: RVarT m Float floatStdNormal = runZiggurat floatStdNormalZ  -- floatStdNormalC must not be over 2^9 if using word32ToFloatWithExcess@@ -188,6 +188,7 @@         getIU         (normalTail floatStdNormalR)     where+        getIU :: RVarT m (Int, Float)         getIU = do             !w <- getRandomPrim PrimWord32             let (u,i) = word32ToFloatWithExcess w@@ -204,16 +205,14 @@     | Normal a a -- mean, sd  instance Distribution Normal Double where-    {-# SPECIALIZE instance Distribution Normal Double #-}-    rvar StdNormal = doubleStdNormal-    rvar (Normal m s) = do+    rvarT StdNormal = doubleStdNormal+    rvarT (Normal m s) = do         x <- doubleStdNormal         return (x * s + m)  instance Distribution Normal Float where-    {-# SPECIALIZE instance Distribution Normal Float #-}-    rvar StdNormal = floatStdNormal-    rvar (Normal m s) = do+    rvarT StdNormal = floatStdNormal+    rvarT (Normal m s) = do         x <- floatStdNormal         return (x * s + m) @@ -227,6 +226,14 @@ stdNormal :: Distribution Normal a => RVar a stdNormal = rvar StdNormal +-- |'stdNormalT' is a normal process with distribution 'StdNormal'.+stdNormalT :: Distribution Normal a => RVarT m a+stdNormalT = rvarT StdNormal+ -- |@normal m s@ is a random variable with distribution @'Normal' m s@. normal :: Distribution Normal a => a -> a -> RVar a normal m s = rvar (Normal m s)++-- |@normalT m s@ is a random process with distribution @'Normal' m s@.+normalT :: Distribution Normal a => a -> a -> RVarT m a+normalT m s = rvarT (Normal m s)

@@ -20,18 +20,18 @@ import Control.Monad  -- from Knuth, with interpretation help from gsl sources-integralPoisson :: (Integral a, RealFloat b, Distribution StdUniform b, Distribution (Erlang a) b, Distribution (Binomial b) a) => b -> RVar a+integralPoisson :: (Integral a, RealFloat b, Distribution StdUniform b, Distribution (Erlang a) b, Distribution (Binomial b) a) => b -> RVarT m a integralPoisson = psn 0     where-        psn :: (Integral a, RealFloat b, Distribution StdUniform b, Distribution (Erlang a) b, Distribution (Binomial b) a) => a -> b -> RVar a+        psn :: (Integral a, RealFloat b, Distribution StdUniform b, Distribution (Erlang a) b, Distribution (Binomial b) a) => a -> b -> RVarT m a         psn j mu             | mu > 10   = do                 let m = floor (mu * (7/8))             -                x <- erlang m+                x <- erlangT m                 if x >= mu                     then do-                        b <- binomial (m - 1) (mu / x)+                        b <- binomialT (m - 1) (mu / x)                         return (j + b)                     else psn (j + m) (mu - x)             @@ -40,21 +40,21 @@                     emu = exp (-mu)                                      prod p k = do-                        u <- stdUniform+                        u <- stdUniformT                         if p * u > emu                             then prod (p * u) (k + 1)                             else return k  integralPoissonCDF :: (Integral a, Real b) => b -> a -> Double integralPoissonCDF mu k = exp (negate lambda) * sum-    [ lambda ^^ i / i_fac-    | (i, i_fac) <- zip [0..k] (scanl (*) 1 [1..])+    [ exp (fromIntegral i * log lambda - i_fac_ln)+    | (i, i_fac_ln) <- zip [0..k] (scanl (+) 0 (map log [1..]))     ]          where lambda = realToFrac mu -fractionalPoisson :: (Num a, Distribution (Poisson b) Integer) => b -> RVar a-fractionalPoisson mu = liftM fromInteger (poisson mu)+fractionalPoisson :: (Num a, Distribution (Poisson b) Integer) => b -> RVarT m a+fractionalPoisson mu = liftM fromInteger (poissonT mu)  fractionalPoissonCDF :: (CDF (Poisson b) Integer, RealFrac a) => b -> a -> Double fractionalPoissonCDF mu k = cdf (Poisson mu) (floor k :: Integer)@@ -62,6 +62,9 @@ poisson :: (Distribution (Poisson b) a) => b -> RVar a poisson mu = rvar (Poisson mu) +poissonT :: (Distribution (Poisson b) a) => b -> RVarT m a+poissonT mu = rvarT (Poisson mu)+ data Poisson b a = Poisson b  $( replicateInstances ''Int integralTypes [d|@@ -70,14 +73,14 @@                  , Distribution (Erlang Int) b                  , Distribution (Binomial b) Int                  ) => Distribution (Poisson b) Int where-            rvar (Poisson mu) = integralPoisson mu+            rvarT (Poisson mu) = integralPoisson mu         instance (Real b, Distribution (Poisson b) Int) => CDF (Poisson b) Int where             cdf  (Poisson mu) = integralPoissonCDF mu     |] )  $( replicateInstances ''Float realFloatTypes [d|         instance (Distribution (Poisson b) Integer) => Distribution (Poisson b) Float where-            rvar (Poisson mu) = fractionalPoisson mu+            rvarT (Poisson mu) = fractionalPoisson mu         instance (CDF (Poisson b) Integer) => CDF (Poisson b) Float where             cdf  (Poisson mu) = fractionalPoissonCDF mu     |])

@@ -13,9 +13,9 @@ import Data.Random.Distribution import Data.Random.Distribution.Uniform -floatingRayleigh :: (Floating a, Distribution StdUniform a) => a -> RVar a+floatingRayleigh :: (Floating a, Distribution StdUniform a) => a -> RVarT m a floatingRayleigh s = do-    u <- stdUniformPos+    u <- stdUniformPosT     return (s * sqrt (-2 * log u))  -- |The rayleigh distribution with a specified mode (\"sigma\") parameter.@@ -27,11 +27,14 @@ rayleigh :: Distribution Rayleigh a => a -> RVar a rayleigh = rvar . Rayleigh +rayleighT :: Distribution Rayleigh a => a -> RVarT m a+rayleighT = rvarT . Rayleigh+ rayleighCDF :: Real a => a -> a -> Double rayleighCDF s x = 1 - exp ((-0.5)* realToFrac (x*x) / realToFrac (s*s))  instance (RealFloat a, Distribution StdUniform a) => Distribution Rayleigh a where-    rvar (Rayleigh s) = floatingRayleigh s+    rvarT (Rayleigh s) = floatingRayleigh s  instance (Real a, Distribution Rayleigh a) => CDF Rayleigh a where     cdf  (Rayleigh s) x = rayleighCDF s x

@@ -29,13 +29,13 @@     deriving (Eq, Show)  -- |Compute a triangular distribution for a 'Floating' type.-floatingTriangular :: (Floating a, Ord a, Distribution StdUniform a) => a -> a -> a -> RVar a+floatingTriangular :: (Floating a, Ord a, Distribution StdUniform a) => a -> a -> a -> RVarT m a floatingTriangular a b c     | a > b     = floatingTriangular b a c     | b > c     = floatingTriangular a c b     | otherwise = do         let p = (c-b)/(c-a)-        u <- stdUniform+        u <- stdUniformT         let d   | u >= p    = a                 | otherwise = c             x   | u >= p    = (u - p) / (1 - p)@@ -57,6 +57,6 @@     = 1      instance (RealFloat a, Ord a, Distribution StdUniform a) => Distribution Triangular a where-    rvar (Triangular a b c) = floatingTriangular a b c+    rvarT (Triangular a b c) = floatingTriangular a b c instance (RealFrac a, Distribution Triangular a) => CDF Triangular a where     cdf  (Triangular a b c) = triangularCDF a b c

@@ -12,10 +12,13 @@ module Data.Random.Distribution.Uniform     ( Uniform(..) 	, uniform+	, uniformT 	     , StdUniform(..)     , stdUniform+    , stdUniformT     , stdUniformPos+    , stdUniformPosT          , integralUniform     , realFloatUniform@@ -51,21 +54,21 @@  -- |Compute a random 'Integral' value between the 2 values provided (inclusive). {-# INLINE integralUniform #-}-integralUniform :: (Integral a) => a -> a -> RVar a+integralUniform :: (Integral a) => a -> a -> RVarT m a integralUniform !x !y = if x < y then integralUniform' x y else integralUniform' y x -{-# SPECIALIZE integralUniform' :: Int   -> Int   -> RVar Int   #-}-{-# SPECIALIZE integralUniform' :: Int8  -> Int8  -> RVar Int8  #-}-{-# SPECIALIZE integralUniform' :: Int16 -> Int16 -> RVar Int16 #-}-{-# SPECIALIZE integralUniform' :: Int32 -> Int32 -> RVar Int32 #-}-{-# SPECIALIZE integralUniform' :: Int64 -> Int64 -> RVar Int64 #-}-{-# SPECIALIZE integralUniform' :: Word   -> Word   -> RVar Word   #-}-{-# SPECIALIZE integralUniform' :: Word8  -> Word8  -> RVar Word8  #-}-{-# SPECIALIZE integralUniform' :: Word16 -> Word16 -> RVar Word16 #-}-{-# SPECIALIZE integralUniform' :: Word32 -> Word32 -> RVar Word32 #-}-{-# SPECIALIZE integralUniform' :: Word64 -> Word64 -> RVar Word64 #-}-{-# SPECIALIZE integralUniform' :: Integer -> Integer -> RVar Integer #-}-integralUniform' :: (Integral a) => a -> a -> RVar a+{-# SPECIALIZE integralUniform' :: Int     -> Int     -> RVarT m Int   #-}+{-# SPECIALIZE integralUniform' :: Int8    -> Int8    -> RVarT m Int8  #-}+{-# SPECIALIZE integralUniform' :: Int16   -> Int16   -> RVarT m Int16 #-}+{-# SPECIALIZE integralUniform' :: Int32   -> Int32   -> RVarT m Int32 #-}+{-# SPECIALIZE integralUniform' :: Int64   -> Int64   -> RVarT m Int64 #-}+{-# SPECIALIZE integralUniform' :: Word    -> Word    -> RVarT m Word   #-}+{-# SPECIALIZE integralUniform' :: Word8   -> Word8   -> RVarT m Word8  #-}+{-# SPECIALIZE integralUniform' :: Word16  -> Word16  -> RVarT m Word16 #-}+{-# SPECIALIZE integralUniform' :: Word32  -> Word32  -> RVarT m Word32 #-}+{-# SPECIALIZE integralUniform' :: Word64  -> Word64  -> RVarT m Word64 #-}+{-# SPECIALIZE integralUniform' :: Integer -> Integer -> RVarT m Integer #-}+integralUniform' :: (Integral a) => a -> a -> RVarT m a integralUniform' !l !u     | nReject == 0  = fmap shift prim     | otherwise     = fmap shift loop@@ -109,29 +112,29 @@  -- |Compute a random value for a 'Bounded' 'Enum' type, between 'minBound' and -- 'maxBound' (inclusive)-boundedEnumStdUniform :: (Enum a, Bounded a) => RVar a+boundedEnumStdUniform :: (Enum a, Bounded a) => RVarT m a boundedEnumStdUniform = enumUniform minBound maxBound  boundedEnumStdUniformCDF :: (Enum a, Bounded a, Ord a) => a -> Double boundedEnumStdUniformCDF = enumUniformCDF minBound maxBound  -- |Compute a uniform random 'Float' value in the range [0,1)-floatStdUniform :: RVar Float+floatStdUniform :: RVarT m Float floatStdUniform = do     x <- getRandomPrim PrimWord32     return (word32ToFloat x)  -- |Compute a uniform random 'Double' value in the range [0,1) {-# INLINE doubleStdUniform #-}-doubleStdUniform :: RVar Double+doubleStdUniform :: RVarT m Double doubleStdUniform = getRandomPrim PrimDouble  -- |Compute a uniform random value in the range [0,1) for any 'RealFloat' type -realFloatStdUniform :: RealFloat a => RVar a+realFloatStdUniform :: RealFloat a => RVarT m a realFloatStdUniform = do     let (b, e) = decodeFloat one     -    x <- uniform 0 (b-1)+    x <- uniformT 0 (b-1)     if x == 0         then return (0 `asTypeOf` one)         else return (encodeFloat x e)@@ -140,12 +143,12 @@  -- |Compute a uniform random 'Fixed' value in the range [0,1), with any -- desired precision.-fixedStdUniform :: HasResolution r => RVar (Fixed r)+fixedStdUniform :: HasResolution r => RVarT m (Fixed r) fixedStdUniform = x     where         res = resolutionOf2 x         x = do-            u <- uniform 0 (res)+            u <- uniformT 0 (res)             return (mkFixed u)  -- |The CDF of the random variable 'realFloatStdUniform'.@@ -156,7 +159,7 @@     | otherwise = realToFrac x  -- |@floatUniform a b@ computes a uniform random 'Float' value in the range [a,b)-floatUniform :: Float -> Float -> RVar Float+floatUniform :: Float -> Float -> RVarT m Float floatUniform 0 1 = floatStdUniform floatUniform a b = do     x <- floatStdUniform@@ -164,7 +167,7 @@  -- |@doubleUniform a b@ computes a uniform random 'Double' value in the range [a,b) {-# INLINE doubleUniform #-}-doubleUniform :: Double -> Double -> RVar Double+doubleUniform :: Double -> Double -> RVarT m Double doubleUniform 0 1 = doubleStdUniform doubleUniform a b = do     x <- doubleStdUniform@@ -172,7 +175,7 @@  -- |@realFloatUniform a b@ computes a uniform random value in the range [a,b) for -- any 'RealFloat' type-realFloatUniform :: RealFloat a => a -> a -> RVar a+realFloatUniform :: RealFloat a => a -> a -> RVarT m a realFloatUniform 0 1 = realFloatStdUniform realFloatUniform a b = do     x <- realFloatStdUniform@@ -180,7 +183,7 @@  -- |@fixedUniform a b@ computes a uniform random 'Fixed' value in the range  -- [a,b), with any desired precision.-fixedUniform :: HasResolution r => Fixed r -> Fixed r -> RVar (Fixed r)+fixedUniform :: HasResolution r => Fixed r -> Fixed r -> RVarT m (Fixed r) fixedUniform a b = do     u <- integralUniform (unMkFixed a) (unMkFixed b)     return (mkFixed u)@@ -195,7 +198,7 @@  -- |@realFloatUniform a b@ computes a uniform random value in the range [a,b) for -- any 'Enum' type-enumUniform :: Enum a => a -> a -> RVar a+enumUniform :: Enum a => a -> a -> RVarT m a enumUniform a b = do     x <- integralUniform (fromEnum a) (fromEnum b)     return (toEnum x)@@ -209,13 +212,19 @@          where e2f = fromIntegral . fromEnum --- @uniform a b@ computes a uniformly distributed random value in the range+-- @uniform a b@ is a uniformly distributed random variable in the range -- [a,b] for 'Integral' or 'Enum' types and in the range [a,b) for 'Fractional' -- types.  Requires a @Distribution Uniform@ instance for the type. uniform :: Distribution Uniform a => a -> a -> RVar a uniform a b = rvar (Uniform a b) --- |Get a \"standard\" uniformly distributed value.+-- @uniformT a b@ is a uniformly distributed random process in the range+-- [a,b] for 'Integral' or 'Enum' types and in the range [a,b) for 'Fractional'+-- types.  Requires a @Distribution Uniform@ instance for the type.+uniformT :: Distribution Uniform a => a -> a -> RVarT m a+uniformT a b = rvarT (Uniform a b)++-- |Get a \"standard\" uniformly distributed variable. -- For integral types, this means uniformly distributed over the full range -- of the type (there is no support for 'Integer').  For fractional -- types, this means uniformly distributed on the interval [0,1).@@ -224,16 +233,29 @@ stdUniform :: (Distribution StdUniform a) => RVar a stdUniform = rvar StdUniform +-- |Get a \"standard\" uniformly distributed process.+-- For integral types, this means uniformly distributed over the full range+-- of the type (there is no support for 'Integer').  For fractional+-- types, this means uniformly distributed on the interval [0,1).+{-# SPECIALIZE stdUniformT :: RVarT m Double #-}+{-# SPECIALIZE stdUniformT :: RVarT m Float #-}+stdUniformT :: (Distribution StdUniform a) => RVarT m a+stdUniformT = rvarT StdUniform+ -- |Like 'stdUniform', but returns only positive or zero values.  Not  -- exported because it is not truly uniform: nonzero values are twice -- as likely as zero on signed types.-stdUniformNonneg :: (Distribution StdUniform a, Num a) => RVar a-stdUniformNonneg = fmap abs stdUniform+stdUniformNonneg :: (Distribution StdUniform a, Num a) => RVarT m a+stdUniformNonneg = fmap abs stdUniformT  -- |Like 'stdUniform' but only returns positive values. stdUniformPos :: (Distribution StdUniform a, Num a) => RVar a-stdUniformPos = iterateUntil (/= 0) stdUniformNonneg+stdUniformPos = stdUniformPosT +-- |Like 'stdUniform' but only returns positive values.+stdUniformPosT :: (Distribution StdUniform a, Num a) => RVarT m a+stdUniformPosT = iterateUntil (/= 0) stdUniformNonneg+ -- |A definition of a uniform distribution over the type @t@.  See also 'uniform'. data Uniform t =      -- |A uniform distribution defined by a lower and upper range bound.@@ -252,8 +274,8 @@ data StdUniform t = StdUniform  $( replicateInstances ''Int integralTypes [d|-        instance Distribution Uniform Int   where rvar (Uniform a b) = integralUniform a b-        instance CDF Uniform Int            where cdf  (Uniform a b) = integralUniformCDF a b+        instance Distribution Uniform Int   where rvarT (Uniform a b) = integralUniform a b+        instance CDF Uniform Int            where cdf   (Uniform a b) = integralUniformCDF a b     |])  instance Distribution StdUniform Word8      where rvarT ~StdUniform = getRandomPrim PrimWord8@@ -267,20 +289,16 @@ instance Distribution StdUniform Int64      where rvarT ~StdUniform = fromIntegral `fmap` getRandomPrim PrimWord64  instance Distribution StdUniform Int where-    rvar-        | toInteger (maxBound :: Int) > toInteger (maxBound :: Int32)-        = const (fromIntegral `fmap` getRandomPrim PrimWord64)-        -        | otherwise-        = const (fromIntegral `fmap` getRandomPrim PrimWord32)+    rvar ~StdUniform =+        $(if toInteger (maxBound :: Int) > toInteger (maxBound :: Int32)+            then [|fromIntegral `fmap` getRandomPrim PrimWord64|]+            else [|fromIntegral `fmap` getRandomPrim PrimWord32|])  instance Distribution StdUniform Word where-    rvar-        | toInteger (maxBound :: Word) > toInteger (maxBound :: Word32)-        = const (fromIntegral `fmap` getRandomPrim PrimWord64)-        -        | otherwise-        = const (fromIntegral `fmap` getRandomPrim PrimWord32)+    rvar ~StdUniform =+        $(if toInteger (maxBound :: Word) > toInteger (maxBound :: Word32)+            then [|fromIntegral `fmap` getRandomPrim PrimWord64|]+            else [|fromIntegral `fmap` getRandomPrim PrimWord32|])  -- Integer has no StdUniform... @@ -289,30 +307,30 @@     |])  -instance Distribution Uniform Float         where rvar (Uniform a b) = floatUniform  a b-instance Distribution Uniform Double        where rvar (Uniform a b) = doubleUniform a b-instance CDF Uniform Float                  where cdf  (Uniform a b) = realUniformCDF a b-instance CDF Uniform Double                 where cdf  (Uniform a b) = realUniformCDF a b+instance Distribution Uniform Float         where rvarT (Uniform a b) = floatUniform  a b+instance Distribution Uniform Double        where rvarT (Uniform a b) = doubleUniform a b+instance CDF Uniform Float                  where cdf   (Uniform a b) = realUniformCDF a b+instance CDF Uniform Double                 where cdf   (Uniform a b) = realUniformCDF a b -instance Distribution StdUniform Float      where rvar ~StdUniform = floatStdUniform-instance Distribution StdUniform Double     where rvar ~StdUniform = getRandomPrim PrimDouble; rvarT ~StdUniform = getRandomPrim PrimDouble-instance CDF StdUniform Float               where cdf  ~StdUniform = realStdUniformCDF-instance CDF StdUniform Double              where cdf  ~StdUniform = realStdUniformCDF+instance Distribution StdUniform Float      where rvarT ~StdUniform = floatStdUniform+instance Distribution StdUniform Double     where rvarT ~StdUniform = getRandomPrim PrimDouble; rvarT ~StdUniform = getRandomPrim PrimDouble+instance CDF StdUniform Float               where cdf   ~StdUniform = realStdUniformCDF+instance CDF StdUniform Double              where cdf   ~StdUniform = realStdUniformCDF  instance HasResolution r => -         Distribution Uniform (Fixed r)     where rvar (Uniform a b) = fixedUniform  a b+         Distribution Uniform (Fixed r)     where rvarT (Uniform a b) = fixedUniform  a b instance HasResolution r => -         CDF Uniform (Fixed r)              where cdf  (Uniform a b) = realUniformCDF a b+         CDF Uniform (Fixed r)              where cdf   (Uniform a b) = realUniformCDF a b instance HasResolution r =>-         Distribution StdUniform (Fixed r)  where rvar ~StdUniform = fixedStdUniform+         Distribution StdUniform (Fixed r)  where rvarT ~StdUniform = fixedStdUniform instance HasResolution r => -         CDF StdUniform (Fixed r)           where cdf  ~StdUniform = realStdUniformCDF+         CDF StdUniform (Fixed r)           where cdf   ~StdUniform = realStdUniformCDF -instance Distribution Uniform ()            where rvar (Uniform _ _) = return ()-instance CDF Uniform ()                     where cdf  (Uniform _ _) = return 1+instance Distribution Uniform ()            where rvarT (Uniform _ _) = return ()+instance CDF Uniform ()                     where cdf   (Uniform _ _) = return 1 $( replicateInstances ''Char [''Char, ''Bool, ''Ordering] [d|-        instance Distribution Uniform Char  where rvar (Uniform a b) = enumUniform a b-        instance CDF Uniform Char           where cdf  (Uniform a b) = enumUniformCDF a b+        instance Distribution Uniform Char  where rvarT (Uniform a b) = enumUniform a b+        instance CDF Uniform Char           where cdf   (Uniform a b) = enumUniformCDF a b      |]) @@ -321,8 +339,8 @@ instance Distribution StdUniform Bool       where rvarT ~StdUniform = fmap even (getRandomPrim PrimWord8) instance CDF StdUniform Bool                where cdf   ~StdUniform = boundedEnumStdUniformCDF -instance Distribution StdUniform Char       where rvar ~StdUniform = boundedEnumStdUniform-instance CDF StdUniform Char                where cdf  ~StdUniform = boundedEnumStdUniformCDF-instance Distribution StdUniform Ordering   where rvar ~StdUniform = boundedEnumStdUniform-instance CDF StdUniform Ordering            where cdf  ~StdUniform = boundedEnumStdUniformCDF+instance Distribution StdUniform Char       where rvarT ~StdUniform = boundedEnumStdUniform+instance CDF StdUniform Char                where cdf   ~StdUniform = boundedEnumStdUniformCDF+instance Distribution StdUniform Ordering   where rvarT ~StdUniform = boundedEnumStdUniform+instance CDF StdUniform Ordering            where cdf   ~StdUniform = boundedEnumStdUniformCDF 

@@ -1,5 +1,6 @@ {-# LANGUAGE         MultiParamTypeClasses,+        RankNTypes,         FlexibleInstances, FlexibleContexts,         RecordWildCards, BangPatterns   #-}@@ -34,7 +35,6 @@ import Data.Vector.Generic as Vec import qualified Data.Vector as V import qualified Data.Vector.Unboxed as UV-import Data.Function (fix)  -- |A data structure containing all the data that is needed -- to implement Marsaglia & Tang's \"ziggurat\" algorithm for@@ -73,18 +73,18 @@         -- single random word (64 bits) can be efficiently converted to         -- a double (using 52 bits) and a bin number (using up to 12 bits),         -- for example.-        zGetIU            :: !(RVar (Int, t)),+        zGetIU            :: !(forall m. RVarT m (Int, t)),                  -- |The distribution for the final \"virtual\" bin         -- (the ziggurat algorithm does not handle distributions         -- that wander off to infinity, so another distribution is needed         -- to handle the last \"bin\" that stretches to infinity)-        zTailDist         :: (RVar t),+        zTailDist         :: (forall m. RVarT m t),                  -- |A copy of the uniform RVar generator for the base type,         -- so that @Distribution Uniform t@ is not needed when sampling         -- from a Ziggurat (makes it a bit more self-contained).-        zUniform          :: !(t -> t -> RVar t),+        zUniform          :: !(forall m. t -> t -> RVarT m t),                  -- |The (one-sided antitone) PDF, not necessarily normalized         zFunc             :: !(t -> t),@@ -99,12 +99,12 @@  -- |Sample from the distribution encoded in a 'Ziggurat' data structure. {-# INLINE runZiggurat #-}-{-# SPECIALIZE runZiggurat :: Ziggurat UV.Vector Float  -> RVar Float #-}-{-# SPECIALIZE runZiggurat :: Ziggurat UV.Vector Double -> RVar Double #-}-{-# SPECIALIZE runZiggurat :: Ziggurat  V.Vector Float  -> RVar Float #-}-{-# SPECIALIZE runZiggurat :: Ziggurat  V.Vector Double -> RVar Double #-}+{-# SPECIALIZE runZiggurat :: Ziggurat UV.Vector Float  -> RVarT m Float #-}+{-# SPECIALIZE runZiggurat :: Ziggurat UV.Vector Double -> RVarT m Double #-}+{-# SPECIALIZE runZiggurat :: Ziggurat  V.Vector Float  -> RVarT m Float #-}+{-# SPECIALIZE runZiggurat :: Ziggurat  V.Vector Double -> RVarT m Double #-} runZiggurat :: (Num a, Ord a, Vector v a) =>-               Ziggurat v a -> RVar a+               Ziggurat v a -> RVarT m a runZiggurat !Ziggurat{..} = go     where         {-# NOINLINE go #-}@@ -166,10 +166,10 @@ --  * an RVar sampling from the tail (the region where x > R) --  {-# INLINE mkZiggurat_ #-}-{-# SPECIALIZE mkZiggurat_ :: Bool -> (Float  ->  Float) -> (Float  ->  Float) -> Int -> Float  -> Float  -> RVar (Int,  Float) -> RVar Float  -> Ziggurat UV.Vector Float #-}-{-# SPECIALIZE mkZiggurat_ :: Bool -> (Double -> Double) -> (Double -> Double) -> Int -> Double -> Double -> RVar (Int, Double) -> RVar Double -> Ziggurat UV.Vector Double #-}-{-# SPECIALIZE mkZiggurat_ :: Bool -> (Float  ->  Float) -> (Float  ->  Float) -> Int -> Float  -> Float  -> RVar (Int,  Float) -> RVar Float  -> Ziggurat V.Vector Float #-}-{-# SPECIALIZE mkZiggurat_ :: Bool -> (Double -> Double) -> (Double -> Double) -> Int -> Double -> Double -> RVar (Int, Double) -> RVar Double -> Ziggurat V.Vector Double #-}+{-# SPECIALIZE mkZiggurat_ :: Bool -> (Float  ->  Float) -> (Float  ->  Float) -> Int -> Float  -> Float  -> (forall m. RVarT m (Int,  Float)) -> (forall m. RVarT m Float ) -> Ziggurat UV.Vector Float #-}+{-# SPECIALIZE mkZiggurat_ :: Bool -> (Double -> Double) -> (Double -> Double) -> Int -> Double -> Double -> (forall m. RVarT m (Int, Double)) -> (forall m. RVarT m Double) -> Ziggurat UV.Vector Double #-}+{-# SPECIALIZE mkZiggurat_ :: Bool -> (Float  ->  Float) -> (Float  ->  Float) -> Int -> Float  -> Float  -> (forall m. RVarT m (Int,  Float)) -> (forall m. RVarT m Float ) -> Ziggurat V.Vector Float #-}+{-# SPECIALIZE mkZiggurat_ :: Bool -> (Double -> Double) -> (Double -> Double) -> Int -> Double -> Double -> (forall m. RVarT m (Int, Double)) -> (forall m. RVarT m Double) -> Ziggurat V.Vector Double #-} mkZiggurat_ :: (RealFloat t, Vector v t,                Distribution Uniform t) =>               Bool@@ -178,15 +178,15 @@               -> Int               -> t               -> t-              -> RVar (Int, t)-              -> RVar t+              -> (forall m. RVarT m (Int, t))+              -> (forall m. RVarT m t)               -> Ziggurat v t mkZiggurat_ m f fInv c r v getIU tailDist = Ziggurat     { zTable_xs         = xs     , zTable_y_ratios   = precomputeRatios xs     , zTable_ys         = Vec.map f xs     , zGetIU            = getIU-    , zUniform          = uniform+    , zUniform          = uniformT     , zFunc             = f     , zTailDist         = tailDist     , zMirror           = m@@ -209,8 +209,8 @@               -> (t -> t)               -> t               -> Int-              -> RVar (Int, t)-              -> (t -> RVar t)+              -> (forall m. RVarT m (Int, t))+              -> (forall m. t -> RVarT m t)               -> Ziggurat v t mkZiggurat m f fInv fInt fVol c getIU tailDist =     mkZiggurat_ m f fInv c r v getIU (tailDist r) @@ -246,10 +246,12 @@   -> (t -> t)   -> t   -> Int-  -> RVar (Int, t)+  -> (forall m. RVarT m (Int, t))   -> Ziggurat v t mkZigguratRec m f fInv fInt fVol c getIU = z         where+            fix :: ((forall m. a -> RVarT m a) -> (forall m. a -> RVarT m a)) -> (forall m. a -> RVarT m a)+            fix f = f (fix f)             z = mkZiggurat m f fInv fInt fVol c getIU (fix (mkTail m f fInv fInt fVol c getIU z))  mkTail :: @@ -258,12 +260,12 @@     -> (a -> a) -> (a -> a) -> (a -> a)     -> a     -> Int-    -> RVar (Int, a)+    -> (forall m. RVarT m (Int, a))     -> Ziggurat v a-    -> (a -> RVar a)-    -> (a -> RVar a)+    -> (forall m. a -> RVarT m a)+    -> (forall m. a -> RVarT m a) mkTail m f fInv fInt fVol c getIU typeRep nextTail r = do-     x <- rvar (mkZiggurat m f' fInv' fInt' fVol' c getIU nextTail `asTypeOf` typeRep)+     x <- rvarT (mkZiggurat m f' fInv' fInt' fVol' c getIU nextTail `asTypeOf` typeRep)      return (x + r * signum x)         where             fIntR = fInt r

@@ -4,30 +4,30 @@ import Data.Fixed import Unsafe.Coerce -#ifdef base_4_2+#ifdef old_Fixed -- So much for backward compatibility through base (>=5) ...  resolutionOf :: HasResolution r => f r -> Integer-resolutionOf = resolution+resolutionOf x = resolution (res x)+    where+        res :: HasResolution r => f r -> r+        res = undefined  resolutionOf2 :: HasResolution r => f (g r) -> Integer resolutionOf2 x = resolution (res x)     where-        res :: HasResolution r => f (g r) -> g r+        res :: HasResolution r => f (g r) -> r         res = undefined  #else  resolutionOf :: HasResolution r => f r -> Integer-resolutionOf x = resolution (res x)-    where-        res :: HasResolution r => f r -> r-        res = undefined+resolutionOf = resolution  resolutionOf2 :: HasResolution r => f (g r) -> Integer resolutionOf2 x = resolution (res x)     where-        res :: HasResolution r => f (g r) -> r+        res :: HasResolution r => f (g r) -> g r         res = undefined  #endif

@@ -17,7 +17,7 @@ -- to make the flexibility this system provides worth the overhead.  I hope -- this is not the case, but if it turns out to be a major problem, this -- system may disappear or be modified in significant ways.-module Data.Random.Internal.Primitives (Prim(..), decomposePrimWhere) where+module Data.Random.Internal.Primitives (Prim(..), getPrimWhere, decomposePrimWhere) where  import Data.Random.Internal.Words import Data.Word@@ -65,6 +65,20 @@     showsPrec _p PrimDouble              = showString "PrimDouble"     showsPrec  p (PrimNByteInteger n)    = showParen (p > 10) (showString "PrimNByteInteger " . showsPrec 11 n) +-- |This function wraps up the most common calling convention for 'decomposePrimWhere'.+-- Given a predicate identifying \"supported\" 'Prim's, and a (possibly partial) +-- function that maps those 'Prim's to implementations, derives a total function+-- mapping all 'Prim's to implementations.+{-# INLINE getPrimWhere #-}+{-# SPECIALIZE getPrimWhere :: Monad m => (forall t. Prim t -> Bool) -> (forall t. Prim t -> m t) -> Prim Word8   -> m Word8   #-}+{-# SPECIALIZE getPrimWhere :: Monad m => (forall t. Prim t -> Bool) -> (forall t. Prim t -> m t) -> Prim Word16  -> m Word16  #-}+{-# SPECIALIZE getPrimWhere :: Monad m => (forall t. Prim t -> Bool) -> (forall t. Prim t -> m t) -> Prim Word32  -> m Word32  #-}+{-# SPECIALIZE getPrimWhere :: Monad m => (forall t. Prim t -> Bool) -> (forall t. Prim t -> m t) -> Prim Word64  -> m Word64  #-}+{-# SPECIALIZE getPrimWhere :: Monad m => (forall t. Prim t -> Bool) -> (forall t. Prim t -> m t) -> Prim Double  -> m Double  #-}+{-# SPECIALIZE getPrimWhere :: Monad m => (forall t. Prim t -> Bool) -> (forall t. Prim t -> m t) -> Prim Integer -> m Integer #-}+getPrimWhere :: Monad m => (forall t. Prim t -> Bool) -> (forall t. Prim t -> m t) -> Prim a -> m a+getPrimWhere supported getPrim prim = runPromptM getPrim (decomposePrimWhere supported prim)+ -- |This is essentially a suite of interrelated default implementations, -- each definition making use of only \"supported\" primitives.  It _really_ -- ought to be inlined to the point where the @supported@ predicate@@ -74,7 +88,7 @@ -- static set of "best" definitions for each required primitive in terms of  -- only supported primitives. -- --- Hopefully, when not inlined, it does not impose too much overhead.+-- Hopefully it does not impose too much overhead when not inlined. {-# INLINE decomposePrimWhere #-} {-# SPECIALIZE decomposePrimWhere :: (forall t. Prim t -> Bool) -> Prim Word8   -> Prompt Prim Word8   #-} {-# SPECIALIZE decomposePrimWhere :: (forall t. Prim t -> Bool) -> Prim Word16  -> Prompt Prim Word16  #-}

@@ -12,7 +12,7 @@ -- anything extra at runtime  {-# INLINE buildWord16 #-}--- |Build a word out of 8 bytes.  No promises are made regarding the order+-- |Build a word out of 2 bytes.  No promises are made regarding the order -- in which the bytes are stuffed.  Note that this means that a 'RandomSource' -- or 'MonadRandom' making use of the default definition of 'getRandomWord', etc., -- may return different random values on different platforms when started @@ -25,7 +25,7 @@         peek (castPtr p)  {-# INLINE buildWord32 #-}--- |Build a word out of 8 bytes.  No promises are made regarding the order+-- |Build a word out of 4 bytes.  No promises are made regarding the order -- in which the bytes are stuffed.  Note that this means that a 'RandomSource' -- or 'MonadRandom' making use of the default definition of 'getRandomWord', etc., -- may return different random values on different platforms when started 

@@ -49,19 +49,24 @@ --  -- > logNormal = exp <$> stdNormal ----- Once defined (in any style), there are a couple ways to sample 'RVar's:+-- Once defined (in any style), there are several ways to sample 'RVar's: --  -- * In a monad, using a 'RandomSource': --  -- > sampleFrom DevRandom (uniform 1 100) :: IO Int -- +-- * In a monad, using a 'MonadRandom' instance:+--+-- > sample (uniform 1 100) :: State PureMT Int+--  -- * As a pure function transforming a functional RNG: --  -- > sampleState (uniform 1 100) :: StdGen -> (Int, StdGen) type RVar = RVarT Identity  -- |\"Run\" an 'RVar' - samples the random variable from the provided--- source of entropy.+-- source of entropy.  Typically 'sample', 'sampleFrom' or 'sampleState' will+-- be more convenient to use. runRVar :: RandomSource m s => RVar a -> s -> m a runRVar = runRVarT @@ -190,10 +195,6 @@     liftIO = T.lift . T.liftIO  instance MonadRandom (RVarT n) where-    supportedPrims _ _ = True-    {-# INLINE getSupportedRandomPrim #-}-    getSupportedRandomPrim p    = RVarT (prompt p)-    {-# INLINE getRandomPrim #-}     getRandomPrim p = RVarT (prompt p)  -- I would really like to be able to do this, but I can't because of the

@@ -12,8 +12,6 @@     ) where  import Data.Word-import Control.Monad.Prompt-import Data.Tagged  import Data.Random.Internal.Primitives @@ -23,119 +21,80 @@ -- when directly requesting entropy for a random variable these functions -- are used. -- --- The minimal definition is 'supportedPrims' and 'getSupportedRandomPrim'--- with cases for those primitives where 'supportedPrims' returns 'True'.------ It is recommended (despite the warnings it generates) that, even when--- all primitives are supported, a final wildcard case of 'supportedPrims' is--- specified, as:+-- Occasionally one might want a 'RandomSource' specifying the 'MonadRandom'+-- instance (for example, when using 'runRVar').  For those cases, +-- "Data.Random.Source.Std".'StdRandom' provides a 'RandomSource' that+-- maps to the 'MonadRandom' instance. -- --- > supportedPrims _ _ = False------ The overlapping pattern warnings can be suppressed (without suppressing --- other, genuine, overlapping-pattern warnings) by the GHC flag--- @-fno-warn-simple-patterns@.  This is not actually the documented behavior--- of that flag as far as I can find in 3 google-minutes, but it works with--- GHC 6.12.1 anyway, and that's good enough for me.------ Note that it is very important that at least 'supportedPrims' (and preferably--- 'getSupportedRandomPrim' as well) gets inlined into the default implementation--- of 'getRandomPrim'.  If your 'supportedPrims' is more than about 2 or 3--- cases, add an INLINE pragma so that it can be optimized out of 'getRandomPrim'.+-- For example, @State StdGen@ has a 'MonadRandom' instance, so to run an+-- 'RVar' (called @x@ in this example) in this monad one could write+-- @runRVar x StdRandom@ (or more concisely with the 'sample' function: @sample x@).+--  class Monad m => MonadRandom m where-    -- |Predicate indicating whether a given primitive is supported by the-    -- instance.  The first parameter is a phantom used to select the instance.-    supportedPrims :: m () -> Prim t -> Bool-    -    -- |Generate a random value corresponding to the specified primitive.  Will-    -- not be called unless supportedPrims returns true for that primitive.-    getSupportedRandomPrim :: Prim t -> m t--    -- This could just be a function, but placing it in a dictionary gives-    -- GHC a place to optimize it separately for each instance, which is -    -- kinda the whole point of the 'Prim' machinery:-    -- -    -- |Generate a random value corresponding to the specified primitive.  The-    -- default implementation makes use of 'supportedPrims' and 'getSupportedRandomPrim'-    -- to construct any required Prim out of the supported ones.-    {-# NOINLINE getRandomPrim #-}+    -- |Generate a random value corresponding to the specified primitive.+    -- The 'Prim' type has many variants, and is also somewhat unstable.+    -- 'getPrimWhere' is a useful function for abstracting over the type,+    -- semi-automatically extending a partial implementation to the full+    -- 'Prim' type.     getRandomPrim :: Prim t -> m t-    getRandomPrim prim = val-        where-            val = runPromptM getSupportedRandomPrim (decomposePrimWhere (supportedPrims mPhantom) prim)-            mPhantom = error "supportedPrims tried to evaluate a phantom parameter" `asTypeOf` (val >> return ())  -- |A source of entropy which can be used in the given monad.------ The minimal definition is 'supportedPrimsFrom' and 'getSupportedRandomPrimFrom'--- with cases for those primitives where 'supportedPrimsFrom' returns 'True'. -- --- Note that it is very important that at least 'supportedPrimsFrom' (and preferably--- 'getSupportedRandomPrimFrom' as well) gets inlined into the default implementation--- of 'getRandomPrimFrom'.  If your 'supportedPrimsFrom' is more than about 2 or 3--- cases, add an INLINE pragma so that it can be optimized out of 'getRandomPrimFrom'.---  -- See also 'MonadRandom'. class Monad m => RandomSource m s where-    -- |Predicate indicating whether a given primitive is supported by the-    -- instance.  The tag on the first parameter is a phantom used only to-    -- select the instance, but the value itself may be inspected.-    supportedPrimsFrom :: Tagged (m ()) s -> Prim t -> Bool-    -    -- |Generate a random value corresponding to the specified primitive-    getSupportedRandomPrimFrom :: s -> Prim t -> m t-    -    -    -- This could just be a function, but placing it in a dictionary gives-    -- GHC a place to optimize it separately for each instance, which is -    -- kinda the whole point of the 'Prim' machinery:-    -- -    -- |Generate a random value corresponding to the specified primitive.  The-    -- default implementation makes use of 'supportedPrimsFrom' and-    -- 'getSupportedRandomPrimFrom' to construct any required Prim out of -    -- the supported ones.-    {-# NOINLINE getRandomPrimFrom #-}+    -- |Generate a random value corresponding to the specified primitive.+    -- The 'Prim' type has many variants, and is also somewhat unstable.+    -- 'getPrimWhere' is a useful function for abstracting over the type,+    -- semi-automatically extending a partial implementation to the full+    -- 'Prim' type.     getRandomPrimFrom :: s -> Prim t -> m t-    getRandomPrimFrom src prim = val-        where-            val = runPromptM (getSupportedRandomPrimFrom src) (decomposePrimWhere supported prim)-            supported :: Prim t -> Bool-            supported = supportedPrimsFrom (tagIt (val >> return ()) src)-            -            tagIt :: a -> b -> Tagged a b-            tagIt _ it = Tagged it  instance Monad m => RandomSource m (m Word8) where-    supportedPrimsFrom _ PrimWord8 = True-    supportedPrimsFrom _ _ = False-    -    getSupportedRandomPrimFrom f PrimWord8 = f-    getSupportedRandomPrimFrom _ p = error ("getSupportedRandomPrimFrom/RandomSource m (m Word8): unsupported prim requested: " ++ show p)+    getRandomPrimFrom f = getPrimWhere supported (getPrim f)+        where+            supported :: Prim a -> Bool+            supported PrimWord8 = True+            supported _ = False+            +            getPrim :: m Word8 -> Prim a -> m a+            getPrim f PrimWord8 = f  instance Monad m => RandomSource m (m Word16) where-    supportedPrimsFrom _ PrimWord16 = True-    supportedPrimsFrom _ _ = False-    -    getSupportedRandomPrimFrom f PrimWord16 = f-    getSupportedRandomPrimFrom _ p = error ("getSupportedRandomPrimFrom/RandomSource m (m Word16): unsupported prim requested: " ++ show p)+    getRandomPrimFrom f = getPrimWhere supported (getPrim f)+        where+            supported :: Prim a -> Bool+            supported PrimWord16 = True+            supported _ = False+            +            getPrim :: m Word16 -> Prim a -> m a+            getPrim f PrimWord16 = f  instance Monad m => RandomSource m (m Word32) where-    supportedPrimsFrom _ PrimWord32 = True-    supportedPrimsFrom _ _ = False-    -    getSupportedRandomPrimFrom f PrimWord32 = f-    getSupportedRandomPrimFrom _ p = error ("getSupportedRandomPrimFrom/RandomSource m (m Word32): unsupported prim requested: " ++ show p)+    getRandomPrimFrom f = getPrimWhere supported (getPrim f)+        where+            supported :: Prim a -> Bool+            supported PrimWord32 = True+            supported _ = False+            +            getPrim :: m Word32 -> Prim a -> m a+            getPrim f PrimWord32 = f  instance Monad m => RandomSource m (m Word64) where-    supportedPrimsFrom _ PrimWord64 = True-    supportedPrimsFrom _ _ = False-    -    getSupportedRandomPrimFrom f PrimWord64 = f-    getSupportedRandomPrimFrom _ p = error ("getSupportedRandomPrimFrom/RandomSource m (m Word64): unsupported prim requested: " ++ show p)+    getRandomPrimFrom f = getPrimWhere supported (getPrim f)+        where+            supported :: Prim a -> Bool+            supported PrimWord64 = True+            supported _ = False+            +            getPrim :: m Word64 -> Prim a -> m a+            getPrim f PrimWord64 = f  instance Monad m => RandomSource m (m Double) where-    supportedPrimsFrom _ PrimDouble = True-    supportedPrimsFrom _ _ = False-    -    getSupportedRandomPrimFrom f PrimDouble = f-    getSupportedRandomPrimFrom _ p = error ("getSupportedRandomPrimFrom/RandomSource m (m Double): unsupported prim requested: " ++ show p)+    getRandomPrimFrom f = getPrimWhere supported (getPrim f)+        where+            supported :: Prim a -> Bool+            supported PrimDouble = True+            supported _ = False+            +            getPrim :: m Double -> Prim a -> m a+            getPrim f PrimDouble = f

@@ -10,6 +10,7 @@     ) where  import Data.Random.Source+import Data.Random.Internal.Primitives  import System.IO (openBinaryFile, hGetBuf, Handle, IOMode(..)) import Foreign@@ -36,18 +37,26 @@ dev DevURandom = devURandom  instance RandomSource IO DevRandom where-    supportedPrimsFrom _ PrimWord8          = True-    supportedPrimsFrom _ PrimWord32         = True-    supportedPrimsFrom _ PrimWord64         = True-    supportedPrimsFrom _ _ = False-    -    getSupportedRandomPrimFrom src PrimWord8  = allocaBytes 1 $ \buf -> do-        1 <- hGetBuf (dev src) buf  1-        peek buf-    getSupportedRandomPrimFrom src PrimWord32  = allocaBytes 1 $ \buf -> do-        4 <- hGetBuf (dev src) buf  4-        peek (castPtr buf)-    getSupportedRandomPrimFrom src PrimWord64  = allocaBytes 8 $ \buf -> do-        8 <- hGetBuf (dev src) buf  8-        peek (castPtr buf)-    getSupportedRandomPrimFrom src prim = error ("getSupportedRandomPrimFrom/" ++ show src ++ ": unsupported prim requested: " ++ show prim)+    getRandomPrimFrom src = getPrimWhere supported getPrim+        where+            supported :: Prim a -> Bool+            supported PrimWord8          = True+            supported PrimWord16         = True+            supported PrimWord32         = True+            supported PrimWord64         = True+            supported _ = False+            +            getPrim :: Prim a -> IO a+            getPrim PrimWord8  = allocaBytes 1 $ \buf -> do+                1 <- hGetBuf (dev src) buf  1+                peek buf+            getPrim PrimWord16 = allocaBytes 2 $ \buf -> do+                2 <- hGetBuf (dev src) buf  2+                peek (castPtr buf)+            getPrim PrimWord32  = allocaBytes 4 $ \buf -> do+                4 <- hGetBuf (dev src) buf  4+                peek (castPtr buf)+            getPrim PrimWord64  = allocaBytes 8 $ \buf -> do+                8 <- hGetBuf (dev src) buf  8+                peek (castPtr buf)+            getPrim prim = error ("getRandomPrimFrom/" ++ show src ++ ": unsupported prim requested: " ++ show prim)

@@ -17,30 +17,25 @@ import Control.Monad.ST  instance RandomSource (ST s) (Gen s) where-    {-# INLINE supportedPrimsFrom #-}-    supportedPrimsFrom _ PrimWord8  = True-    supportedPrimsFrom _ PrimWord16 = True-    supportedPrimsFrom _ PrimWord32 = True-    supportedPrimsFrom _ PrimWord64 = True-    supportedPrimsFrom _ PrimDouble = True-    supportedPrimsFrom _ _ = False+    getRandomPrimFrom src = getPrimWhere supported (getPrim src)+        where+            {-# INLINE supported #-}+            supported :: Prim a -> Bool+            supported PrimWord8  = True+            supported PrimWord16 = True+            supported PrimWord32 = True+            supported PrimWord64 = True+            supported PrimDouble = True+            supported _ = False     -    {-# INLINE getSupportedRandomPrimFrom #-}-    getSupportedRandomPrimFrom gen PrimWord8    = uniform gen-    getSupportedRandomPrimFrom gen PrimWord16   = uniform gen-    getSupportedRandomPrimFrom gen PrimWord32   = uniform gen-    getSupportedRandomPrimFrom gen PrimWord64   = uniform gen-    getSupportedRandomPrimFrom gen PrimDouble   = fmap wordToDouble (uniform gen)-    getSupportedRandomPrimFrom _ p = error ("getSupportedRandomPrimFrom/Gen s: unsupported prim requested: " ++ show p)+            {-# INLINE getPrim #-}+            getPrim :: Gen s -> Prim a -> ST s a+            getPrim gen PrimWord8    = uniform gen+            getPrim gen PrimWord16   = uniform gen+            getPrim gen PrimWord32   = uniform gen+            getPrim gen PrimWord64   = uniform gen+            getPrim gen PrimDouble   = fmap wordToDouble (uniform gen)+            getPrim gen p            = error ("getSupportedRandomPrimFrom/Gen s: unsupported prim requested: " ++ show p)  instance RandomSource IO (Gen RealWorld) where-    {-# INLINE supportedPrimsFrom #-}-    supportedPrimsFrom _ PrimWord8  = True-    supportedPrimsFrom _ PrimWord16 = True-    supportedPrimsFrom _ PrimWord32 = True-    supportedPrimsFrom _ PrimWord64 = True-    supportedPrimsFrom _ PrimDouble = True-    supportedPrimsFrom _ _ = False-    -    {-# INLINE getSupportedRandomPrimFrom #-}-    getSupportedRandomPrimFrom gen prim = stToIO (getSupportedRandomPrimFrom gen prim)+    getRandomPrimFrom src = stToIO . getRandomPrimFrom src

@@ -1,9 +1,10 @@ {-# LANGUAGE+    CPP,     BangPatterns,     MultiParamTypeClasses,     FlexibleContexts, FlexibleInstances,     UndecidableInstances,-    GADTs+    GADTs, RankNTypes   #-} {-# OPTIONS_GHC -fno-warn-orphans #-} @@ -12,7 +13,18 @@ -- values (the pure pseudorandom generator provided by the -- mersenne-random-pure64 package), as well as instances for some common -- cases.-module Data.Random.Source.PureMT where+-- +-- A 'PureMT' generator is immutable, so 'PureMT' by itself cannot be a +-- 'RandomSource' (if it were, it would always give the same \"random\"+-- values).  Some form of mutable state must be used, such as an 'IORef',+-- 'State' monad, etc..  A few default instances are provided by this module+-- along with more-general functions ('getRandomPrimFromMTRef' and+-- 'getRandomPrimFromMTState') usable as implementations for new cases+-- users might need.+module Data.Random.Source.PureMT +    ( PureMT, newPureMT, pureMT+    , module Data.Random.Source.PureMT +    ) where  import Data.Random.Internal.Primitives import Data.Random.Source@@ -20,29 +32,60 @@  import Data.StateRef -import Control.Monad.Prompt import Control.Monad.State import qualified Control.Monad.ST.Strict as S import qualified Control.Monad.State.Strict as S --- |Given a mutable reference to a 'PureMT' generator, we can implement--- 'RandomSource' for in any monad in which the reference can be modified.-getRandomPrimFromMTRef :: (Monad m, ModifyRef sr m PureMT) => sr -> Prim a -> m a-getRandomPrimFromMTRef ref prim-    | supported prim = getThing (genPrim prim)-    | otherwise = runPromptM (getRandomPrimFromMTRef ref) (decomposePrimWhere supported prim)+-- |Given a function for applying a 'PureMT' transformation to some hidden +-- state, this function derives a function able to generate all 'Prim's+-- in the given monad.  This is then suitable for either a 'MonadRandom' or+-- 'RandomSource' instance, where the 'supportedPrims' or+-- 'supportedPrimsFrom' function (respectively) is @const True@.+{-# INLINE getRandomPrimBy #-}+getRandomPrimBy :: Monad m => (forall t. (PureMT -> (t, PureMT)) -> m t) -> Prim a -> m a+getRandomPrimBy getThing = getPrimWhere supported (\prim -> getThing (genPrim prim))     where +        {-# INLINE supported #-}         supported :: Prim a -> Bool         supported PrimWord64 = True         supported PrimDouble = True         supported _          = False         +        {-# INLINE genPrim #-}         genPrim :: Prim a -> (PureMT -> (a, PureMT))         genPrim PrimWord64 = randomWord64         genPrim PrimDouble = randomDouble-        genPrim p = error ("getRandomPrimFromMTRef: genPrim called for unsupported prim " ++ show p)-        -        getThing thing = atomicModifyReference ref $ \(!oldMT) -> case thing oldMT of (!w, !newMT) -> (newMT, w)+        genPrim p = error ("getRandomPrimBy: genPrim called for unsupported prim " ++ show p)++-- |Given a mutable reference to a 'PureMT' generator, we can implement+-- 'RandomSource' for in any monad in which the reference can be modified.+-- +-- Typically this would be used to define a new 'RandomSource' instance for+-- some new reference type or new monad in which an existing reference type+-- can be modified atomically.  As an example, the following instance could+-- be used to describe how 'IORef' 'PureMT' can be a 'RandomSource' in the+-- 'IO' monad:+-- +-- > instance RandomSource IO (IORef PureMT) where+-- >     supportedPrimsFrom _ _ = True+-- >     getSupportedRandomPrimFrom = getRandomPrimFromMTRef+-- +-- (note that there is actually a more general instance declared already+-- covering this as a a special case, so there's no need to repeat this+-- declaration anywhere)+-- +-- Example usage:+-- +-- > main = do+-- >     src <- newIORef (pureMT 1234)          -- OR: newPureMT >>= newIORef+-- >     x <- sampleFrom src (uniform 0 100)    -- OR: runRVar (uniform 0 100) src+-- >     print x+getRandomPrimFromMTRef :: (Monad m, ModifyRef sr m PureMT) => sr -> Prim a -> m a+getRandomPrimFromMTRef ref = getRandomPrimBy getThing+    where+        {-# INLINE getThing #-}+        getThing thing = atomicModifyReference ref $ \(!oldMT) -> +            case thing oldMT of (!w, !newMT) -> (newMT, w)               -- |Similarly, @getRandomPrimFromMTState x@ can be used in any \"state\"@@ -53,57 +96,58 @@ -- @runState . sample :: Distribution d t => d t -> PureMT -> (t, PureMT)@. -- 'PureMT' in the type there can be replaced by 'StdGen' or anything else  -- satisfying @MonadRandom (State s) => s@).+-- +-- For example, this module includes the following declaration:+-- +-- > instance MonadRandom (State PureMT) where+-- >     supportedPrims _ _ = True+-- >     getSupportedRandomPrim = getRandomPrimFromMTState+-- +-- This describes a \"standard\" way of getting random values in 'State'+-- 'PureMT', which can then be used in various ways, for example (assuming +-- some 'RVar' @foo@ and some 'Word64' @seed@):+-- +-- > runState (runRVar foo StdRandom) (pureMT seed)+-- > runState (sampleFrom StdRandom foo) (pureMT seed)+-- > runState (sample foo) (pureMT seed)+-- +-- Of course, the initial 'PureMT' state could also be obtained by any other+-- convenient means, such as 'newPureMT' if you don't care what seed is used. getRandomPrimFromMTState :: MonadState PureMT m => Prim a -> m a-getRandomPrimFromMTState prim-    | supported prim = getThing (genPrim prim)-    | otherwise = runPromptM getRandomPrimFromMTState (decomposePrimWhere supported prim)+getRandomPrimFromMTState = getRandomPrimBy getThing     where-        supported :: Prim a -> Bool-        supported PrimWord64 = True-        supported PrimDouble = True-        supported _          = False-        -        genPrim :: Prim a -> (PureMT -> (a, PureMT))-        genPrim PrimWord64 = randomWord64-        genPrim PrimDouble = randomDouble-        genPrim p = error ("getRandomPrimFromMTRef: genPrim called for unsupported prim " ++ show p)-        +        {-# INLINE getThing #-}         getThing thing = do             !mt <- get             let (!ws, !newMt) = thing mt             put newMt             return ws +#ifndef MTL2 instance MonadRandom (State PureMT) where-    supportedPrims _ _ = True-    getSupportedRandomPrim = getRandomPrimFromMTState+    getRandomPrim = getRandomPrimFromMTState  instance MonadRandom (S.State PureMT) where-    supportedPrims _ _ = True-    getSupportedRandomPrim = getRandomPrimFromMTState+    getRandomPrim = getRandomPrimFromMTState+#endif  instance (Monad m1, ModifyRef (Ref m2 PureMT) m1 PureMT) => RandomSource m1 (Ref m2 PureMT) where-    supportedPrimsFrom _ _ = True-    getSupportedRandomPrimFrom = getRandomPrimFromMTRef+    getRandomPrimFrom = getRandomPrimFromMTRef      instance Monad m => MonadRandom (StateT PureMT m) where-    supportedPrims _ _ = True-    getSupportedRandomPrim = getRandomPrimFromMTState+    getRandomPrim = getRandomPrimFromMTState  instance Monad m => MonadRandom (S.StateT PureMT m) where-    supportedPrims _ _ = True-    getSupportedRandomPrim = getRandomPrimFromMTState+    getRandomPrim = getRandomPrimFromMTState  instance (Monad m, ModifyRef (IORef PureMT) m PureMT) => RandomSource m (IORef PureMT) where-    {-# SPECIALIZE instance RandomSource IO (IORef PureMT)#-}-    supportedPrimsFrom _ _ = True-    getSupportedRandomPrimFrom = getRandomPrimFromMTRef+    {-# SPECIALIZE instance RandomSource IO (IORef PureMT) #-}+    getRandomPrimFrom = getRandomPrimFromMTRef      instance (Monad m, ModifyRef (STRef s PureMT) m PureMT) => RandomSource m (STRef s PureMT) where     {-# SPECIALIZE instance RandomSource (ST s) (STRef s PureMT) #-}     {-# SPECIALIZE instance RandomSource (S.ST s) (STRef s PureMT) #-}-    supportedPrimsFrom _ _ = True-    getSupportedRandomPrimFrom = getRandomPrimFromMTRef+    getRandomPrimFrom = getRandomPrimFromMTRef  -- Note that this instance is probably a Bad Idea.  STM allows random variables -- to interact in spooky quantum-esque ways - One transaction can 'retry' until@@ -112,6 +156,5 @@ -- instance (Monad m, ModifyRef (TVar PureMT) m PureMT) => RandomSource m (TVar PureMT) where --     {-# SPECIALIZE instance RandomSource IO  (TVar PureMT) #-} --     {-# SPECIALIZE instance RandomSource STM (TVar PureMT) #-}---     supportedPrimsFrom _ _ = True---     getSupportedRandomPrimFrom = getRandomPrimFromMTRef+--     getRandomPrimFrom = getRandomPrimFromMTRef     

@@ -8,7 +8,6 @@ module Data.Random.Source.Std where  import Data.Random.Source-import Data.Tagged  -- |A token representing the \"standard\" entropy source in a 'MonadRandom' -- monad.  Its sole purpose is to make the following true (when the types check):@@ -18,10 +17,4 @@  instance MonadRandom m => RandomSource m StdRandom where     {-SPECIALIZE instance MonadRandom m => RandomSource m StdRandom -}-    supportedPrimsFrom w = supportedPrims (mkWit w)-        where-            mkWit :: Tagged a b -> a-            mkWit = error "supportedPrims tried to evaluate its phantom parameter"-    getSupportedRandomPrimFrom   StdRandom = getSupportedRandomPrim-         getRandomPrimFrom StdRandom = getRandomPrim

@@ -1,4 +1,5 @@ {-# LANGUAGE+    CPP,     MultiParamTypeClasses, FlexibleInstances, UndecidableInstances, GADTs,     BangPatterns, RankNTypes   #-}@@ -24,13 +25,11 @@   instance (Monad m1, ModifyRef (Ref m2 StdGen) m1 StdGen) => RandomSource m1 (Ref m2 StdGen) where-    supportedPrimsFrom _ _ = True-    getSupportedRandomPrimFrom = getRandomPrimFromRandomGenRef+    getRandomPrimFrom = getRandomPrimFromRandomGenRef  instance (Monad m, ModifyRef (IORef   StdGen) m StdGen) => RandomSource m (IORef   StdGen) where     {-# SPECIALIZE instance RandomSource IO (IORef StdGen) #-}-    supportedPrimsFrom _ _ = True-    getSupportedRandomPrimFrom = getRandomPrimFromRandomGenRef+    getRandomPrimFrom = getRandomPrimFromRandomGenRef  -- Note that this instance is probably a Bad Idea.  STM allows random variables -- to interact in spooky quantum-esque ways - One transaction can 'retry' until@@ -45,8 +44,7 @@ instance (Monad m, ModifyRef (STRef s StdGen) m StdGen) => RandomSource m (STRef s StdGen) where     {-# SPECIALIZE instance RandomSource (ST s) (STRef s StdGen) #-}     {-# SPECIALIZE instance RandomSource (S.ST s) (STRef s StdGen) #-}-    supportedPrimsFrom _ _ = True-    getSupportedRandomPrimFrom = getRandomPrimFromRandomGenRef+    getRandomPrimFrom = getRandomPrimFromRandomGenRef  getRandomPrimFromStdGenIO :: Prim a -> IO a getRandomPrimFromStdGenIO prim@@ -81,6 +79,10 @@  -- |Given a mutable reference to a 'RandomGen' generator, we can make a -- 'RandomSource' usable in any monad in which the reference can be modified.+-- +-- See "Data.Random.Source.PureMT".'getRandomPrimFromMTRef' for more detailed+-- usage hints - this function serves exactly the same purpose except for a+-- 'StdGen' generator instead of a 'PureMT' generator. getRandomPrimFromRandomGenRef :: (Monad m, ModifyRef sr m g, RandomGen g) =>                                   sr -> Prim a -> m a getRandomPrimFromRandomGenRef ref prim@@ -116,12 +118,15 @@ -- Additionally, the standard mtl state monads have 'MonadRandom' instances -- which do precisely that, allowing an easy conversion of 'RVar's and -- other 'Distribution' instances to \"pure\" random variables.+-- +-- Again, see "Data.Random.Source.PureMT".'getRandomPrimFromMTState' for more+-- detailed usage hints - this function serves exactly the same purpose except +-- for a 'StdGen' generator instead of a 'PureMT' generator. {-# SPECIALIZE getRandomPrimFromRandomGenState :: Prim a -> State StdGen a #-} {-# SPECIALIZE getRandomPrimFromRandomGenState :: Monad m => Prim a -> StateT StdGen m a #-} getRandomPrimFromRandomGenState :: (RandomGen g, MonadState g m) => Prim a -> m a getRandomPrimFromRandomGenState prim-    | supported prim = genSupported prim-    | otherwise = runPromptM genSupported (decomposePrimWhere supported prim)+    = runPromptM genSupported (decomposePrimWhere supported prim)     where          {-# INLINE genSupported #-}         genSupported prim = genPrim prim getThing@@ -158,19 +163,17 @@                     put newGen                     return (f $! i) +#ifndef MTL2 instance MonadRandom (State StdGen) where-    supportedPrims _ _ = True-    getSupportedRandomPrim = getRandomPrimFromRandomGenState--instance Monad m => MonadRandom (StateT StdGen m) where-    supportedPrims _ _ = True-    getSupportedRandomPrim = getRandomPrimFromRandomGenState+    getRandomPrim = getRandomPrimFromRandomGenState  instance MonadRandom (S.State StdGen) where-    supportedPrims _ _ = True-    getSupportedRandomPrim = getRandomPrimFromRandomGenState+    getRandomPrim = getRandomPrimFromRandomGenState+#endif +instance Monad m => MonadRandom (StateT StdGen m) where+    getRandomPrim = getRandomPrimFromRandomGenState+ instance Monad m => MonadRandom (S.StateT StdGen m) where-    supportedPrims _ _ = True-    getSupportedRandomPrim = getRandomPrimFromRandomGenState+    getRandomPrim = getRandomPrimFromRandomGenState