random-fu 0.0.0.2 → 0.0.1.1

raw patch · 26 files changed

+879/−479 lines, 26 filesdep +containersdep −bytestringPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: containers

Dependencies removed: bytestring

API changes (from Hackage documentation)

Files

• random-fu.cabal +5/−3
• src/Data/Random.hs +6/−1
• src/Data/Random/Distribution.hs +12/−15
• src/Data/Random/Distribution.hs-boot +0/−10
• src/Data/Random/Distribution/Bernoulli.hs +40/−17
• src/Data/Random/Distribution/Beta.hs +12/−12
• src/Data/Random/Distribution/Binomial.hs +34/−18
• src/Data/Random/Distribution/Discrete.hs +71/−7
• src/Data/Random/Distribution/Exponential.hs +6/−6
• src/Data/Random/Distribution/Gamma.hs +61/−15
• src/Data/Random/Distribution/Normal.hs +30/−11
• src/Data/Random/Distribution/Poisson.hs +21/−17
• src/Data/Random/Distribution/Triangular.hs +6/−5
• src/Data/Random/Distribution/Uniform.hs +75/−31
• src/Data/Random/Internal/Classification.hs +0/−102
• src/Data/Random/Internal/Words.hs +24/−20
• src/Data/Random/Lift.hs +29/−0
• src/Data/Random/List.hs +67/−0
• src/Data/Random/RVar.hs +76/−56
• src/Data/Random/RVar.hs-boot +0/−16
• src/Data/Random/Sample.hs +35/−0
• src/Data/Random/Source.hs +56/−41
• src/Data/Random/Source/DevRandom.hs +26/−9
• src/Data/Random/Source/PureMT.hs +93/−20
• src/Data/Random/Source/Std.hs +3/−2
• src/Data/Random/Source/StdGen.hs +91/−45

@@ -1,5 +1,5 @@ name:                   random-fu-version:                0.0.0.2+version:                0.0.1.1 stability:              experimental  cabal-version:          >= 1.2@@ -34,9 +34,11 @@                         Data.Random.Distribution.Poisson                         Data.Random.Distribution.Triangular                         Data.Random.Distribution.Uniform-                        Data.Random.Internal.Classification                         Data.Random.Internal.Words+                        Data.Random.List+                        Data.Random.Lift                         Data.Random.RVar+                        Data.Random.Sample                         Data.Random.Source                         Data.Random.Source.DevRandom                         Data.Random.Source.StdGen@@ -44,7 +46,7 @@                         Data.Random.Source.Std                            build-depends:        base >= 3,-                        bytestring,+                        containers,                         mersenne-random-pure64,                         monad-loops >= 0.3.0.1,                         mtl,

@@ -19,7 +19,8 @@ -- a couple handy 'RVar's.  module Data.Random-    ( module Data.Random.Source+    ( module Data.Random.Sample+    , module Data.Random.Source     , module Data.Random.Source.DevRandom     , module Data.Random.Source.StdGen     , module Data.Random.Source.PureMT@@ -35,9 +36,11 @@     , module Data.Random.Distribution.Poisson     , module Data.Random.Distribution.Triangular     , module Data.Random.Distribution.Uniform+    , module Data.Random.List     , module Data.Random.RVar     ) where +import Data.Random.Sample import Data.Random.Source import Data.Random.Source.DevRandom import Data.Random.Source.StdGen@@ -54,5 +57,7 @@ import Data.Random.Distribution.Poisson import Data.Random.Distribution.Triangular import Data.Random.Distribution.Uniform+import Data.Random.Lift ()+import Data.Random.List import Data.Random.RVar 

@@ -7,26 +7,23 @@  module Data.Random.Distribution where -import {-# SOURCE #-} Data.Random.RVar+import Data.Random.Lift+import Data.Random.RVar import Data.Random.Source import Data.Random.Source.Std import Data.Word  -- |A definition of a random variable's distribution.  From the distribution--- an 'RVar' can be created, or the distribution can be directly sampled.--- 'RVar' in particular is an instance of 'Distribution', and so can be 'sample'd.------ Minimum instance definition: either 'rvar' or 'sampleFrom'.+-- an 'RVar' can be created, or the distribution can be directly sampled using +-- 'sampleFrom' or 'sample'. class Distribution d t where     -- |Return a random variable with this distribution.-    rvar :: d t -> RVar t-    rvar = sampleFrom StdRandom-    -    -- |Directly sample from the distribution, given a source of entropy.-    sampleFrom :: RandomSource m s => s -> d t -> m t-    sampleFrom src dist = sampleFrom src (rvar dist)+    rvar :: (Distribution d t) => d t -> RVar t+    rvar = rvarT --- |Sample a distribution using the default source of entropy for the--- monad in which the sampling occurs.-sample :: (Distribution d t, MonadRandom m) => d t -> m t-sample = sampleFrom StdRandom+-- |Return a random variable with the given distribution, pre-lifted to an arbitrary 'RVarT'.+-- Any arbitrary 'RVar' can also be converted to an 'RVarT m' for an arbitrary 'm', using+-- either 'lift' or 'sample'.+rvarT :: Distribution d t => d t -> RVarT n t+rvarT d = lift (rvar d)+

@@ -1,10 +0,0 @@-{-- -      ``Data/Random/Distribution''- -}-{-# LANGUAGE-    MultiParamTypeClasses, KindSignatures-  #-}--module Data.Random.Distribution where--class Distribution (d :: * -> *) t

@@ -9,8 +9,6 @@  module Data.Random.Distribution.Bernoulli where -import Data.Random.Internal.Classification- import Data.Random.Source import Data.Random.Distribution import Data.Random.RVar@@ -19,29 +17,54 @@  import Data.Int import Data.Word+import Data.Ratio+import Data.Complex -bernoulli :: (Distribution (Bernoulli b) a) => b -> RVar a-bernoulli p = sample (Bernoulli p)+-- |Generate a Bernoulli variate 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.+bernoulli :: Distribution (Bernoulli b) a => b -> RVar a+bernoulli p = rvar (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 p = do-    x <- realFloatUniform 0 1+    x <- stdUniform     return (x <= p) +-- | @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 t f p = do-    x <- boolBernoulli p+    x <- bernoulli p     return (if x then t else f) -class (Classification NumericType t c) => BernoulliByClassification c t where-    bernoulliByClassification :: RealFloat a => a -> RVar t+data Bernoulli b a = Bernoulli b -instance (Classification NumericType t IntegralType, Num t) => BernoulliByClassification IntegralType t-    where bernoulliByClassification = generalBernoulli 0 1-instance (Classification NumericType t FractionalType, Num t) => BernoulliByClassification FractionalType t-    where bernoulliByClassification = generalBernoulli 0 1-instance (Classification NumericType t EnumType, Enum t) => BernoulliByClassification EnumType t-    where bernoulliByClassification = generalBernoulli (toEnum 0) (toEnum 1)+instance (Fractional b, Ord b, Distribution StdUniform b) +       => Distribution (Bernoulli b) Bool+    where+        rvar (Bernoulli p) = boolBernoulli p -data Bernoulli b a = Bernoulli b+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Int         where rvar (Bernoulli p) = generalBernoulli 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Int8        where rvar (Bernoulli p) = generalBernoulli 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Int16       where rvar (Bernoulli p) = generalBernoulli 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Int32       where rvar (Bernoulli p) = generalBernoulli 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Int64       where rvar (Bernoulli p) = generalBernoulli 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Word8       where rvar (Bernoulli p) = generalBernoulli 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Word16      where rvar (Bernoulli p) = generalBernoulli 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Word32      where rvar (Bernoulli p) = generalBernoulli 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Word64      where rvar (Bernoulli p) = generalBernoulli 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Integer     where rvar (Bernoulli p) = generalBernoulli 0 1 p -instance (BernoulliByClassification c t, RealFloat b) => Distribution (Bernoulli b) t where-    rvar (Bernoulli p) = bernoulliByClassification p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Float       where rvar (Bernoulli p) = generalBernoulli 0 1 p+instance Distribution (Bernoulli b) Bool => Distribution (Bernoulli b) Double      where rvar (Bernoulli p) = generalBernoulli 0 1 p+instance (Distribution (Bernoulli b) Bool, Integral a)+    => Distribution (Bernoulli b) (Ratio a)   +    where rvar (Bernoulli p) = generalBernoulli 0 1 p+instance (Distribution (Bernoulli b) Bool, RealFloat a)+    => Distribution (Bernoulli b) (Complex a)+    where rvar (Bernoulli p) = generalBernoulli 0 1 p++

@@ -17,23 +17,23 @@  import Control.Monad -realFloatBeta :: RealFloat a => a -> a -> RVar a-realFloatBeta 1 1 = realFloatStdUniform-realFloatBeta a b = do-    x <- realFloatGamma a 1-    y <- realFloatGamma b 1+fractionalBeta :: (Fractional a, Distribution Gamma a, Distribution StdUniform a) => a -> a -> RVar a+fractionalBeta 1 1 = stdUniform+fractionalBeta a b = do+    x <- gamma a 1+    y <- gamma b 1     return (x / (x + y)) -realFloatBetaFromIntegral :: (Integral a, Integral b, RealFloat c) => a -> b -> RVar c-realFloatBetaFromIntegral a b =  do-    x <- realFloatErlang a-    y <- realFloatErlang b+fractionalBetaFromIntegral :: (Fractional c, Distribution (Erlang a) c, Distribution (Erlang b) c) => a -> b -> RVar c+fractionalBetaFromIntegral a b =  do+    x <- erlang a+    y <- erlang b     return (x / (x + y))  beta :: Distribution Beta a => a -> a -> RVar a-beta a b = sample (Beta a b)+beta a b = rvar (Beta a b)  data Beta a = Beta a a -instance (RealFloat a) => Distribution Beta a where-    rvar (Beta a b) = realFloatBeta a b+instance (Fractional a, Distribution Gamma a, Distribution StdUniform a) => Distribution Beta a where+    rvar (Beta a b) = fractionalBeta a b

@@ -9,8 +9,6 @@  module Data.Random.Distribution.Binomial where -import Data.Random.Internal.Classification- import Data.Random.Source import Data.Random.Distribution import Data.Random.RVar@@ -24,15 +22,28 @@      -- algorithm from Knuth's TAOCP, 3rd ed., p 136     -- specific choice of cutoff size taken from gsl source-integralBinomial :: (Integral a, RealFloat b) => a -> b -> RVar a+    -- 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.+integralBinomial :: (Integral a, Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => a -> b -> RVar a integralBinomial t p = bin 0 t p     where+        -- GHC likes to discharge Beta to the Beta instance's context, which+        -- @integralBinomial@'s context doesn't (directly) satisfy.+        -- Seems like GHC could do better.  GHC could discharge it in @integralBinomial@'s+        -- context when attempting to satisfy bin, since the Beta instance covers all @b@,+        -- whereupon it would find that the contexts do, in fact, match.+        -- Of course, it's a pretty obscure case, so maybe it's better to just +        -- leave it to the coder who's doing weird stuff to tell the compiler what+        -- he/she wants.+        -- Anyway, this type signature makes GHC happy.+        bin :: (Integral a, Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => a -> a -> b -> RVar a         bin k t p             | t > 10    = do                 let a = 1 + t `div` 2                     b = 1 + t - a         -                x <- realFloatBetaFromIntegral a b+                x <- beta (fromIntegral a) (fromIntegral b)                 if x >= p                     then bin  k      (a - 1) (p / x)                     else bin (k + a) (b - 1) ((p - x) / (1 - x))@@ -41,24 +52,29 @@                 where                     count k  0    = return k                     count k (n+1) = do-                        x <- realFloatStdUniform+                        x <- stdUniform                         (count $! (if x < p then k + 1 else k)) n -+-- would it be valid to repeat the above computation using fractional @t@?+-- obviously something different would have to be done with @count@ as well...+floatingBinomial :: (RealFrac a, Distribution (Binomial b) Integer) => a -> b -> RVar a+floatingBinomial t p = fmap fromInteger (rvar (Binomial (truncate t) p))  binomial :: Distribution (Binomial b) a => a -> b -> RVar a-binomial t p = sample (Binomial t p)--class (Classification NumericType t c) => BinomialByClassification c t where-    binomialByClassification :: RealFloat a => t -> a -> RVar t--instance (Classification NumericType t IntegralType, Integral t) => BinomialByClassification IntegralType t-    where binomialByClassification = integralBinomial-instance (Classification NumericType t FractionalType, RealFrac t) => BinomialByClassification FractionalType t-    where binomialByClassification t p = liftM fromInteger (integralBinomial (truncate t) p)--instance (BinomialByClassification c t, RealFloat b) => Distribution (Binomial b) t where-    rvar (Binomial t p) = binomialByClassification t p+binomial t p = rvar (Binomial t p)  data Binomial b a = Binomial a b +instance (Ord b, Floating b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Int        where rvar (Binomial t p) = integralBinomial t p+instance (Ord b, Floating b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Int8       where rvar (Binomial t p) = integralBinomial t p+instance (Ord b, Floating b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Int16      where rvar (Binomial t p) = integralBinomial t p+instance (Ord b, Floating b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Int32      where rvar (Binomial t p) = integralBinomial t p+instance (Ord b, Floating b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Int64      where rvar (Binomial t p) = integralBinomial t p+instance (Ord b, Floating b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Word8      where rvar (Binomial t p) = integralBinomial t p+instance (Ord b, Floating b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Word16     where rvar (Binomial t p) = integralBinomial t p+instance (Ord b, Floating b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Word32     where rvar (Binomial t p) = integralBinomial t p+instance (Ord b, Floating b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Word64     where rvar (Binomial t p) = integralBinomial t p+instance (Ord b, Floating b, Distribution Beta b, Distribution StdUniform b) => Distribution (Binomial b) Integer    where rvar (Binomial t p) = integralBinomial t p++instance Distribution (Binomial b) Integer => Distribution (Binomial b) Float  where rvar (Binomial t p) = floatingBinomial t p+instance Distribution (Binomial b) Integer => Distribution (Binomial b) Double where rvar (Binomial t p) = floatingBinomial t p

@@ -11,13 +11,19 @@ import Data.Random.RVar import Data.Random.Distribution import Data.Random.Distribution.Uniform+import Data.Random.List (randomElement)  import Control.Monad+import Control.Applicative +import Data.List+import Data.Function+ discrete :: Distribution (Discrete p) a => [(p,a)] -> RVar a discrete ps = rvar (Discrete ps) -data Discrete p a = Discrete [(p, a)]+newtype Discrete p a = Discrete [(p, a)]+    deriving (Eq, Show)  instance (Num p, Ord p, Distribution Uniform p) => Distribution (Discrete p) a where     rvar (Discrete []) = fail "discrete distribution over empty set cannot be sampled"@@ -27,9 +33,67 @@                  when (any (<0) ps) $ fail "negative probability in discrete distribution"         -        u <- uniform 0 (last cs)-        return $ head-            [ x-            | (c,x) <- zip cs xs-            , c >= u-            ]+        let totalProb = last cs+        if totalProb <= 0+            then randomElement xs   -- this probably makes the monad instance incorrect for discarding zero-probability events...+            else do+                u <- uniform 0 totalProb+                return $ head+                    [ x+                    | (c,x) <- zip cs xs+                    , c >= u+                    ]++instance Functor (Discrete p) where+    fmap f (Discrete ds) = Discrete [(p, f x) | (p, x) <- ds]++-- TODO - check out whether this is valid when not requiring normalization...+-- We want each subset of cases in fx derived from a given case +-- in x to have the same total probability as the set in x from whence they came.+--+-- thus, w(f x) == w (x) is sufficient (although not necessary), where w() is+-- the weight.+instance (Fractional p, Ord p) => Monad (Discrete p) where+    return x = Discrete [(1, x)]+    (Discrete x) >>= f = Discrete $ do+        (p, x) <- x+        +        let Discrete fx = f x+        let qx = [ (q, x)+                 | (q, x) <- fx+                 , q > 0    -- should this be done?  Consider case where all results have 0 weight...+                 ]+            qs = map fst qx     -- either (qx == []) or (sum qs > 0)+            scale = recip (sum qs)+        +        (q, x) <- qx+        -- now (qx /= []), because ((q,x) `elem` qx)+        -- therefore sum qs > 0, therefore 0 < scale < ∞.+        +        return (p * q * scale, x)++instance (Fractional p, Ord p) => Applicative (Discrete p) where+    pure = return+    (<*>) = ap++collectDiscreteEvents :: (Ord e, Num p, Ord p) => Discrete p e -> Discrete p e+collectDiscreteEvents (Discrete ds) = +    Discrete . concatMap (uncurry combine . unzip) . groupEvents . sortEvents $ ds+    +    where+        groupEvents = groupBy ((==) `on` snd)+        sortEvents  = sortBy (compare `on` snd)+        +        -- don't combine negative weights with positive ones.+        -- don't error out ether - just leave them alone, it'll+        -- all barf when the distribution is sampled.+        combine ps (x:_) = case partition (> 0) (filter (/= 0) ps) of+            ([], []) ->                             []+            ([], ns) ->               (sum ns, x) : []+            (ps, []) -> (sum ps, x) :               []+            (ps, ns) -> (sum ps, x) : (sum ns, x) : []+        +        weight (p,x)+            | p < 0     = error "negative probability in discrete distribution"+            | otherwise = p+        event ((p,x):_) = x

@@ -16,13 +16,13 @@  data Exponential a = Exp a -realFloatExponential :: RealFloat a => a -> RVar a-realFloatExponential lambdaRecip = do-    x <- realFloatStdUniform+floatingExponential :: (Floating a, Distribution StdUniform a) => a -> RVar a+floatingExponential lambdaRecip = do+    x <- stdUniform     return (negate (log x) * lambdaRecip)  exponential :: Distribution Exponential a => a -> RVar a-exponential = sample . Exp+exponential = rvar . Exp -instance (RealFloat a) => Distribution Exponential a where-    rvar (Exp lambdaRecip) = realFloatExponential lambdaRecip+instance (Floating a, Distribution StdUniform a) => Distribution Exponential a where+    rvar (Exp lambdaRecip) = floatingExponential lambdaRecip

@@ -28,21 +28,32 @@     -- originally comes from Marsaglia & Tang, "A Simple Method for     -- generating gamma variables", ACM Transactions on Mathematical     -- Software, Vol 26, No 3 (2000), p363-372.-realFloatGamma :: RealFloat a => a -> a -> RVar a+realFloatGamma :: (Floating a, Ord a, Distribution NormalPair (a,a), Distribution StdUniform a) => a -> a -> RVar a realFloatGamma a b     | a < 1      = do-        u <- realFloatStdUniform+        u <- stdUniform         x <- realFloatGamma (1 + a) b         return (x * u ** recip a)     | otherwise-    = go+    = goNothing         where             d = a - (1 / 3)             c = recip (3 * sqrt d) -- (1 / 3) / sqrt d             -            go = do-                x <- realFloatStdNormal+            -- manually unrolled StateT (Maybe Double)+            -- since the current stdNormal uses a normal pair and+            -- discards the 2nd, this implementation uses+            -- the normal pair and stashes one so that every other+            -- time through the loop it can recycle what+            -- would be discarded anyway.+            goNothing = do+                (x, stashed) <- normalPair+                step x (goJust stashed)+                +            goJust x = step x goNothing+            +            step x next = do                 let cx = c * x                     v = (1 + cx) ^ 3                     @@ -50,25 +61,60 @@                     x_4 = x_2 * x_2                                  if cx <= (-1)-                    then go+                    then next                     else do-                        u <- realFloatStdUniform+                        u <- stdUniform                                                  if         u < 1 - 0.0331 * x_4                             || log u < 0.5 * x_2  + d * (1 - v + log v)                             then return (b * d * v)-                            else go+                            else next -realFloatErlang :: (Integral a, RealFloat b) => a -> RVar b-realFloatErlang a = realFloatGamma (fromIntegral a) 1 +realFloatErlang :: (Integral a, Floating b, Ord b, Distribution NormalPair (b,b), Distribution StdUniform b) => a -> RVar b+realFloatErlang a+    | a < 1 +    = fail "realFloatErlang: a < 1"+    | otherwise+    = goNothing+        where+            d = fromIntegral a - (1 / 3)+            c = recip (3 * sqrt d) -- (1 / 3) / sqrt d+            +            goNothing = do+                (x, stashed) <- normalPair+                step x (goJust stashed)+                +            goJust x = step x goNothing+            +            step x next = do+                let cx = c * x+                    v = (1 + cx) ^ 3+                    +                    x_2 = x * x+                    x_4 = x_2 * x_2+                +                if cx <= (-1)+                    then next+                    else do+                        u <- stdUniform+                        +                        if         u < 1 - 0.0331 * x_4+                            || log u < 0.5 * x_2  + d * (1 - v + log v)+                            then return (d * v)+                            else next+ gamma :: (Distribution Gamma a) => a -> a -> RVar a-gamma a b = sample (Gamma a b)+gamma a b = rvar (Gamma a b) -erlang :: (Distribution Gamma a, Integral b, Num a) => b -> a -> RVar a-erlang a b = sample (Gamma (fromIntegral a) b)+erlang :: (Distribution (Erlang a) b) => a -> RVar b+erlang a = rvar (Erlang a) -data Gamma a = Gamma a a+data Gamma a    = Gamma a a+data Erlang a b = Erlang a -instance RealFloat a => Distribution Gamma a where+instance (Floating a, Ord a, Distribution NormalPair (a,a), Distribution StdUniform a) => Distribution Gamma a where     rvar (Gamma a b) = realFloatGamma a b++instance (Integral a, Floating b, Ord b, Distribution NormalPair (b,b), Distribution StdUniform b) => Distribution (Erlang a) b where+    rvar (Erlang a) = realFloatErlang a

@@ -18,18 +18,22 @@  import Control.Monad --- Box-Muller method-normalPair :: (Floating a, Distribution Uniform a) => RVar (a,a)-normalPair = do-    u <- uniform 0 1-    t <- uniform 0 (2 * pi)+normalPair :: Distribution NormalPair (a, a) => RVar (a,a)+normalPair = rvar NormalPair++{-# INLINE boxMullerNormalPair #-}+boxMullerNormalPair :: (Floating a, Distribution StdUniform a) => RVar (a,a)+boxMullerNormalPair = do+    u <- stdUniform+    t <- stdUniform     let r = sqrt (-2 * log u)+        theta = (2 * pi) * t         -        x = r * cos t-        y = r * sin t+        x = r * cos theta+        y = r * sin theta     return (x,y) --- slightly slower+{-# INLINE knuthPolarNormalPair #-} knuthPolarNormalPair :: (Floating a, Ord a, Distribution Uniform a) => RVar (a,a) knuthPolarNormalPair = do     v1 <- uniform (-1) 1@@ -53,15 +57,30 @@     return x      +doubleStdNormal :: RVar Double+doubleStdNormal = do+    u <- doubleStdUniform+    t <- doubleStdUniform+    let r = sqrt (-2 * log u)+        +        x = r * cos (t * 2 * pi)+    return x+    + data Normal a     = StdNormal     | Normal a a -- mean, sd -instance (Floating a, Distribution Uniform a) => Distribution Normal a where-    rvar StdNormal = liftM fst normalPair+data NormalPair a = NormalPair++instance (Floating a, Ord a, Distribution StdUniform a) => Distribution Normal a where+    rvar StdNormal = liftM fst boxMullerNormalPair     rvar (Normal m s) = do-        x <- liftM fst normalPair+        x <- liftM fst boxMullerNormalPair         return (x * s + m)++instance (Floating a, Distribution StdUniform a) => Distribution NormalPair (a, a) where+    rvar _ = boxMullerNormalPair  stdNormal :: Distribution Normal a => RVar a stdNormal = rvar StdNormal

@@ -22,17 +22,18 @@ import Control.Monad  -- from Knuth, with interpretation help from gsl sources-integralPoisson :: (Integral a, RealFloat b) => b -> RVar a+integralPoisson :: (Integral a, RealFloat b, Distribution StdUniform b, Distribution (Erlang a) b, Distribution (Binomial b) a) => b -> RVar a integralPoisson mu = psn 0 mu     where+        psn :: (Integral a, RealFloat b, Distribution StdUniform b, Distribution (Erlang a) b, Distribution (Binomial b) a) => a -> b -> RVar a         psn k mu             | mu > 10   = do                 let m = floor (mu * (7/8))             -                x <- realFloatErlang m+                x <- erlang m                 if x >= mu                     then do-                        b <- integralBinomial (m - 1) (mu / x)+                        b <- binomial (m - 1) (mu / x)                         return (k + b)                     else psn (k + m) (mu - x)             @@ -41,27 +42,30 @@                     emu = exp (-mu)                                      prod p k = do-                        u <- realFloatStdUniform+                        u <- stdUniform                         if p * u > emu                             then prod (p * u) (k + 1)                             else return k +fractionalPoisson :: (Num a, Distribution (Poisson b) Integer) => b -> RVar a+fractionalPoisson mu = liftM fromInteger (poisson mu)  poisson :: (Distribution (Poisson b) a) => b -> RVar a-poisson mu = sample (Poisson mu)+poisson mu = rvar (Poisson mu)  data Poisson b a = Poisson b -instance RealFloat b => Distribution (Poisson b) Int        where rvar (Poisson mu) = integralPoisson mu-instance RealFloat b => Distribution (Poisson b) Int8       where rvar (Poisson mu) = integralPoisson mu-instance RealFloat b => Distribution (Poisson b) Int16      where rvar (Poisson mu) = integralPoisson mu-instance RealFloat b => Distribution (Poisson b) Int32      where rvar (Poisson mu) = integralPoisson mu-instance RealFloat b => Distribution (Poisson b) Int64      where rvar (Poisson mu) = integralPoisson mu-instance RealFloat b => Distribution (Poisson b) Word8      where rvar (Poisson mu) = integralPoisson mu-instance RealFloat b => Distribution (Poisson b) Word16     where rvar (Poisson mu) = integralPoisson mu-instance RealFloat b => Distribution (Poisson b) Word32     where rvar (Poisson mu) = integralPoisson mu-instance RealFloat b => Distribution (Poisson b) Word64     where rvar (Poisson mu) = integralPoisson mu-instance RealFloat b => Distribution (Poisson b) Integer    where rvar (Poisson mu) = integralPoisson mu+-- so ugly...+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Int     ) b, Distribution (Binomial b) Int     ) => Distribution (Poisson b) Int        where rvar (Poisson mu) = integralPoisson mu+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Int8    ) b, Distribution (Binomial b) Int8    ) => Distribution (Poisson b) Int8       where rvar (Poisson mu) = integralPoisson mu+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Int16   ) b, Distribution (Binomial b) Int16   ) => Distribution (Poisson b) Int16      where rvar (Poisson mu) = integralPoisson mu+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Int32   ) b, Distribution (Binomial b) Int32   ) => Distribution (Poisson b) Int32      where rvar (Poisson mu) = integralPoisson mu+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Int64   ) b, Distribution (Binomial b) Int64   ) => Distribution (Poisson b) Int64      where rvar (Poisson mu) = integralPoisson mu+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Word8   ) b, Distribution (Binomial b) Word8   ) => Distribution (Poisson b) Word8      where rvar (Poisson mu) = integralPoisson mu+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Word16  ) b, Distribution (Binomial b) Word16  ) => Distribution (Poisson b) Word16     where rvar (Poisson mu) = integralPoisson mu+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Word32  ) b, Distribution (Binomial b) Word32  ) => Distribution (Poisson b) Word32     where rvar (Poisson mu) = integralPoisson mu+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Word64  ) b, Distribution (Binomial b) Word64  ) => Distribution (Poisson b) Word64     where rvar (Poisson mu) = integralPoisson mu+instance (RealFloat b, Distribution StdUniform b, Distribution (Erlang Integer ) b, Distribution (Binomial b) Integer ) => Distribution (Poisson b) Integer    where rvar (Poisson mu) = integralPoisson mu -instance RealFloat b => Distribution (Poisson b) Float      where rvar (Poisson mu) = liftM fromIntegral (integralPoisson mu)-instance RealFloat b => Distribution (Poisson b) Double     where rvar (Poisson mu) = liftM fromIntegral (integralPoisson mu)+instance (Distribution (Poisson b) Integer) => Distribution (Poisson b) Float            where rvar (Poisson mu) = fractionalPoisson mu+instance (Distribution (Poisson b) Integer) => Distribution (Poisson b) Double           where rvar (Poisson mu) = fractionalPoisson mu

@@ -3,7 +3,8 @@  -} {-# LANGUAGE     MultiParamTypeClasses,-    FlexibleInstances+    FlexibleInstances, FlexibleContexts,+    UndecidableInstances   #-}  module Data.Random.Distribution.Triangular where@@ -18,19 +19,19 @@     , triUpper  :: a     } deriving (Eq, Show) -realFloatTriangular :: (RealFloat a) => a -> a -> a -> RVar a+realFloatTriangular :: (Floating a, Ord a, Distribution StdUniform a) => a -> a -> a -> RVar a realFloatTriangular a b c     | a <= b && b <= c     = do         let p = (c-b)/(c-a)-        u <- realFloatStdUniform+        u <- stdUniform         let d   | u >= p    = a                 | otherwise = c             x   | u >= p    = (u - p) / (1 - p)                 | otherwise = u / p -- may prefer this: reusing u costs resolution, especially if p or 1-p is small and c-a is large.---        x <- realFloatStdUniform+--        x <- stdUniform         return (b - ((1 - sqrt x) * (b-d))) -instance RealFloat a => Distribution Triangular a where+instance (Floating a, Ord a, Distribution StdUniform a) => Distribution Triangular a where     rvar (Triangular a b c) = realFloatTriangular a b c

@@ -4,32 +4,34 @@ {-# LANGUAGE     MultiParamTypeClasses, FunctionalDependencies,     FlexibleContexts, FlexibleInstances, -    UndecidableInstances+    UndecidableInstances, EmptyDataDecls   #-}  module Data.Random.Distribution.Uniform     ( Uniform(..)-	, UniformByClassification(..) 	, uniform 	     , StdUniform(..)-    , StdUniformByClassification(..)     , stdUniform          , integralUniform     , realFloatUniform+    , floatUniform+    , doubleUniform          , boundedStdUniform     , boundedEnumStdUniform     , realFloatStdUniform+    , floatStdUniform+    , doubleStdUniform     ) where -import Data.Random.Internal.Classification-+import Data.Random.Internal.Words import Data.Random.Source import Data.Random.Distribution import Data.Random.RVar +import Data.Ratio import Data.Word import Data.Int import Data.Bits@@ -37,6 +39,7 @@  import Control.Monad.Loops +integralUniform :: (Integral a) => a -> a -> RVar a integralUniform a b     | a > b     = compute b a     | otherwise = compute a b@@ -61,20 +64,38 @@ boundedEnumStdUniform :: (Enum a, Bounded a) => RVar a boundedEnumStdUniform = enumUniform minBound maxBound --- (0,1]+floatStdUniform :: RVar Float+floatStdUniform = do+    x <- getRandomWord+    return (wordToFloat x)++doubleStdUniform :: RVar Double+doubleStdUniform = getRandomDouble+ realFloatStdUniform :: RealFloat a => RVar a-realFloatStdUniform | False     = return one-                    | otherwise = do+realFloatStdUniform = do     let bitsNeeded  = floatDigits one         (_, e) = decodeFloat one          x <- nBitInteger bitsNeeded     if x == 0-        then return 1+        then return one         else return (encodeFloat x (e-1))          where one = 1 +floatUniform :: Float -> Float -> RVar Float+floatUniform 0 1 = floatStdUniform+floatUniform a b = do+    x <- floatStdUniform+    return (a + x * (b - a))++doubleUniform :: Double -> Double -> RVar Double+doubleUniform 0 1 = doubleStdUniform+doubleUniform a b = do+    x <- doubleStdUniform+    return (a + x * (b - a))+ realFloatUniform :: RealFloat a => a -> a -> RVar a realFloatUniform 0 1 = realFloatStdUniform realFloatUniform a b = do@@ -89,34 +110,57 @@ uniform :: Distribution Uniform a => a -> a -> RVar a uniform a b = rvar (Uniform a b) -stdUniform :: Distribution StdUniform a => RVar a+-- |Get a \"standard\" uniformly distributed value.+-- For integral types, this means uniformly distributed over the full range+-- of the type (and hence there is no support for Integer).  For fractional+-- types, this means uniformly distributed on the interval [0,1).+stdUniform :: (Distribution StdUniform a) => RVar a stdUniform = rvar StdUniform -class (Classification NumericType t c) => UniformByClassification c t where-    uniformByClassification :: t -> t -> RVar t--class (Classification NumericType t c) => StdUniformByClassification c t where-    stdUniformByClassification :: RVar t+stdUniformPos :: (Distribution StdUniform a, Ord a, Num a) => RVar a+stdUniformPos = do+    x <- stdUniform+    if x > 0+        then return x+        else stdUniformPos  data Uniform t = Uniform !t !t data StdUniform t = StdUniform -instance UniformByClassification c t => Distribution Uniform t-    where rvar (Uniform a b) = uniformByClassification a b+instance Distribution Uniform Int           where rvar (Uniform a b) = integralUniform a b+instance Distribution Uniform Int8          where rvar (Uniform a b) = integralUniform a b+instance Distribution Uniform Int16         where rvar (Uniform a b) = integralUniform a b+instance Distribution Uniform Int32         where rvar (Uniform a b) = integralUniform a b+instance Distribution Uniform Int64         where rvar (Uniform a b) = integralUniform a b+instance Distribution Uniform Word8         where rvar (Uniform a b) = integralUniform a b+instance Distribution Uniform Word16        where rvar (Uniform a b) = integralUniform a b+instance Distribution Uniform Word32        where rvar (Uniform a b) = integralUniform a b+instance Distribution Uniform Word64        where rvar (Uniform a b) = integralUniform a b+instance Distribution Uniform Integer       where rvar (Uniform a b) = integralUniform a b -instance StdUniformByClassification c t => Distribution StdUniform t-    where rvar _ = stdUniformByClassification+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 (Classification NumericType t IntegralType, Integral t) => UniformByClassification IntegralType t-    where uniformByClassification = integralUniform-instance (Classification NumericType t FractionalType, RealFloat t) => UniformByClassification FractionalType t-    where uniformByClassification = realFloatUniform-instance (Classification NumericType t EnumType, Enum t) => UniformByClassification EnumType t-    where uniformByClassification = enumUniform+instance Distribution Uniform Char          where rvar (Uniform a b) = enumUniform a b+instance Distribution Uniform Bool          where rvar (Uniform a b) = enumUniform a b+instance Distribution Uniform ()            where rvar (Uniform a b) = enumUniform a b+instance Distribution Uniform Ordering      where rvar (Uniform a b) = enumUniform a b -instance (Classification NumericType t IntegralType, Integral t, Bounded t) => StdUniformByClassification IntegralType t-    where stdUniformByClassification = boundedStdUniform-instance (Classification NumericType t FractionalType, RealFloat t) => StdUniformByClassification FractionalType t-    where stdUniformByClassification = realFloatStdUniform-instance (Classification NumericType t EnumType, Enum t, Bounded t) => StdUniformByClassification EnumType t-    where stdUniformByClassification = boundedStdUniform+instance Distribution StdUniform Int        where rvar StdUniform = fmap fromIntegral getRandomWord+instance Distribution StdUniform Int8       where rvar StdUniform = fmap fromIntegral getRandomByte+instance Distribution StdUniform Int16      where rvar StdUniform = fmap fromIntegral getRandomWord+instance Distribution StdUniform Int32      where rvar StdUniform = fmap fromIntegral getRandomWord+instance Distribution StdUniform Int64      where rvar StdUniform = fmap fromIntegral getRandomWord+instance Distribution StdUniform Word8      where rvar StdUniform = getRandomByte+instance Distribution StdUniform Word16     where rvar StdUniform = fmap fromIntegral getRandomWord+instance Distribution StdUniform Word32     where rvar StdUniform = fmap fromIntegral getRandomWord+instance Distribution StdUniform Word64     where rvar StdUniform = fmap fromIntegral getRandomWord++instance Distribution StdUniform Float      where rvar StdUniform = floatStdUniform+instance Distribution StdUniform Double     where rvar StdUniform = doubleStdUniform++instance Distribution StdUniform Char       where rvar StdUniform = boundedEnumStdUniform+instance Distribution StdUniform Bool       where rvar StdUniform = fmap even getRandomByte+instance Distribution StdUniform ()         where rvar StdUniform = return ()+instance Distribution StdUniform Ordering   where rvar StdUniform = boundedEnumStdUniform+

@@ -1,102 +0,0 @@-{-# LANGUAGE-    MultiParamTypeClasses, FunctionalDependencies,-    EmptyDataDecls-  #-}-{-- -      ``Data/Random/Internal/Classification''- -}      --- | \"Classification systems\" - for a motivating example, see---  the implementation of the Uniform distribution.  Basically,---  I would like to make instances like:---  ---  > instance RealFloat a => Distribution Uniform a where ...---  > instance Integral a =>  Distribution Uniform a where ...---  ---  and so on.  However, this is not sound - what happens if someone---  comes along and makes a type that's an instance of both Integral and---  RealFloat?---  ---  So, we introduce a classification system based on phantom types, so---  that each type can be unambiguously declared to be \"intensionally\"---  Integral, Floating, or whatever.---  ---  Now, obviously it'd be nice not to clutter the Distribution typeclass---  with extra phantom types that the end user shouldn't care about.  Hence---  the pattern of introducing typeclasses such as "UniformByClassification"---  ---  Now, if a new type comes along that is Integral, a single declaration---  of the following form suffices to attach it to all such Distribution---  instances:---  ---  > instance Classification NumericType t IntegralType---  ---  Not quite as automagic as the @Integral a => Distribution foo@  case,---  but a bit closer.  Not only that, it leaves open the possibility that---  a user may bring in a type that is \"mostly\" integral, and has an Integral---  instance, but should be handled differently for purposes of uniform---  random number generation.  In such a case, the user may introduce a new---  classification of their own and provide the required instances for that---  classification.---  ---  All in all, although it is not yet well-tested, it has the \"feel\" of ---  a good compromise.-module Data.Random.Internal.Classification where--import Data.Int-import Data.Word-import Data.Ratio---- |classificiation system, experimental--- ---      * c (a phantom type) is the classification system---      ---      * t is the type to be classified---      ---      * tc (a phantom type) is the classification of t according to c------ The functional dependency, aside from being important because the relation--- is functional, allows the classification system to be \"discharged\" in--- cases such as the following:------ > class Classification SomeCS t c => FooByClassification t c where ...--- > instance FooByClassification t c => Foo t where ...--- --- Thus the class of interest to the end user need not display anything--- at all about the classification system, except in the superclasses of --- the classes in the contexts of some of its instances.-class Classification c t tc | c t -> tc---- |A simple classification system covering the cases we care--- about when sampling distributions.  Loosely, these are the reasons we care:--- ---   * distributions over Fractional types are handled as if the type were continuous.--- ---   * distributions over Integral types are handled discretely.--- ---   * distributions over Enum types (which are not Num instances) are handled ---     like Integral types, but require use of 'fromEnum' and/or 'toEnum' to work with them.-data NumericType--data IntegralType-data FractionalType-data EnumType--instance Classification NumericType Int            IntegralType-instance Classification NumericType Int8           IntegralType-instance Classification NumericType Int16          IntegralType-instance Classification NumericType Int32          IntegralType-instance Classification NumericType Int64          IntegralType-instance Classification NumericType Word8          IntegralType-instance Classification NumericType Word16         IntegralType-instance Classification NumericType Word32         IntegralType-instance Classification NumericType Word64         IntegralType-instance Classification NumericType Integer        IntegralType--instance Classification NumericType Float          FractionalType-instance Classification NumericType Double         FractionalType-instance Classification NumericType (Ratio a)      FractionalType--instance Classification NumericType Char           EnumType-instance Classification NumericType Bool           EnumType-instance Classification NumericType ()             EnumType-instance Classification NumericType Ordering       EnumType

@@ -3,34 +3,38 @@  -}  -- |A few little functions I found myself writing inline over and over again.------ Note that these need to be checked to ensure proper behavior on big-endian --- systems.  They are probably not right at the moment. module Data.Random.Internal.Words where  import Foreign import GHC.IOBase +import Data.Bits import Data.Word import Control.Monad -wordsToBytes :: [Word64] -> [Word8]-wordsToBytes = concatMap wordToBytes--wordToBytes :: Word64 -> [Word8]-wordToBytes x = unsafePerformIO . allocaBytes 8 $ \p -> do-    poke (castPtr p) x-    mapM (peekElemOff p) [0..7]+{-# INLINE buildWord #-}+buildWord :: Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word64+buildWord b0 b1 b2 b3 b4 b5 b6 b7+    = unsafePerformIO . allocaBytes 8 $ \p -> do+        pokeByteOff p 0 b0+        pokeByteOff p 1 b1+        pokeByteOff p 2 b2+        pokeByteOff p 3 b3+        pokeByteOff p 4 b4+        pokeByteOff p 5 b5+        pokeByteOff p 6 b6+        pokeByteOff p 7 b7+        peek (castPtr p) -bytesToWords :: [Word8] -> [Word64]-bytesToWords = map bytesToWord . chunk 8-    where-        chunk n [] = []-        chunk n xs = case splitAt n xs of-            (ys, zs) -> ys : chunk n zs+{-# INLINE wordToFloat #-}+-- |Pack 23 unspecified bits from a 'Word64' into a 'Float' in the range [0,1).+-- Used to convert a 'stdUniform' 'Word64' to a 'stdUniform' 'Double'.+wordToFloat :: Word64 -> Float+wordToFloat x = (encodeFloat $! toInteger (x `shiftR` ( 41 {- 64-23 -}))) $ (-23) +{-# INLINE wordToDouble #-}+-- |Pack 52 unspecified bits from a 'Word64' into a 'Double' in the range [0,1).+-- Used to convert a 'stdUniform' 'Word64' to a 'stdUniform' 'Double'.+wordToDouble :: Word64 -> Double+wordToDouble x = (encodeFloat $! toInteger (x `shiftR` ( 12 {- 64-52 -}))) $ (-52) -bytesToWord :: [Word8] -> Word64-bytesToWord bs = unsafePerformIO . allocaBytes 8 $ \p -> do-    zipWithM (pokeElemOff p) [0..7] (bs ++ repeat 0)-    peek (castPtr p)

@@ -0,0 +1,29 @@+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, IncoherentInstances #-}++module Data.Random.Lift where++import Control.Monad.Identity+import qualified Control.Monad.Trans as T++-- | A class for \"liftable\" data structures. Conceptually+-- an extension of 'T.MonadTrans' to allow deep lifting,+-- but lifting need not be done between monads only. Eg lifting+-- between 'Applicative's is allowed.+--+-- For instances where 'm' and 'n' have 'return'/'pure' defined,+-- these instances must satisfy+-- @lift (return x) == return x@.+class Lift m n where+    lift :: m a -> n a++instance (Monad m, T.MonadTrans t) => Lift m (t m) where+    lift = T.lift++instance Lift m m where+    lift = id++-- | This instance is incoherent with the other two. However,+-- by the law @lift (return x) == return x@, the results+-- must always be the same.+instance Monad m => Lift Identity m where+    lift = return . runIdentity

@@ -0,0 +1,67 @@+{-+ -      ``Data/Random/List''+ -}+{-# LANGUAGE +    FlexibleContexts+  #-}++module Data.Random.List where++import Data.Random.RVar+import Data.Random.Source+import Data.Random.Distribution+import Data.Random.Distribution.Uniform+import GHC.IOBase++import qualified Data.Sequence as S++randomElement :: [a] -> RVar a+randomElement [] = error "randomElement: empty list!"+randomElement xs = do+    n <- uniform 0 (length xs - 1)+    return (xs !! n)++randomSeqElement :: S.Seq a -> RVar a+randomSeqElement s+    | S.null s  = error "randomSeqElement: empty list!"+    | otherwise = do+        n <- uniform 0 (S.length s - 1)+        return (s `S.index` n)++shuffle :: [a] -> RVar [a]+shuffle = shuffleSeq . S.fromList++shuffleSeq :: S.Seq a -> RVar [a]+shuffleSeq s = shuffle (S.length s) s+    where+        shuffle 0 _ = return []+        shuffle (n+1) s = do+            i <- uniform 0 n+            let (x, xs) = extract i s+            ys <- shuffle n xs+            return (x:ys)+        +        extract n s = case S.splitAt n s of+            (l,r) -> case S.viewl r of+                x S.:< r  -> (x, l S.>< r)++-- |Shuffle a list using interleaved IO when extracting elements.+lazyShuffleFrom :: (RandomSource IO s) => s -> [a] -> IO [a]+lazyShuffleFrom src = lazyShuffleSeqFrom src . S.fromList++-- |Shuffle a 'S.Seq' using interleaved IO when extracting elements.+lazyShuffleSeqFrom :: (RandomSource IO s) => s -> S.Seq a -> IO [a]+lazyShuffleSeqFrom src s = shuffle (S.length s) s+    where+        shuffle 0     _  = return []+        shuffle (n+1) s +            | S.null s = return []+            | otherwise = do+                i <- runRVar (uniform 0 n) src+                let (x, xs) = extract i s+                ys <- unsafeInterleaveIO (shuffle n xs)+                return (x:ys)+        +        extract n s = case S.splitAt n s of+            (l,r) -> case S.viewl r of+                x S.:< r  -> (x, l S.>< r)

@@ -4,7 +4,8 @@ {-# LANGUAGE     RankNTypes,     MultiParamTypeClasses,-    FlexibleInstances+    FlexibleInstances, +    GADTs   #-}  -- |Random variables.  An 'RVar' is a sampleable random variable.  Because@@ -14,87 +15,106 @@ -- 'RVar's. module Data.Random.RVar     ( RVar-    +    , runRVar+    , RVarT+    , runRVarT     , nByteInteger     , nBitInteger     ) where -import Data.Random.Distribution++import Data.Random.Internal.Words import Data.Random.Source+import Data.Random.Lift as L  import Data.Word import Data.Bits +import qualified Control.Monad.Trans as T import Control.Applicative import Control.Monad+import Control.Monad.Reader+import Control.Monad.Identity  -- |An opaque type containing a \"random variable\" - a value  -- which depends on the outcome of some random process.-newtype RVar a = RVar { runDistM :: forall m s. RandomSource m s => s -> m a }+type RVar = RVarT Identity -instance Functor RVar where+-- | single combined container allowing all the relevant +-- dictionaries (plus the RandomSource item itself) to be passed+-- with one pointer.+data RVarDict n m where+    RVarDict :: (Lift n m, Monad m, RandomSource m s) => s -> RVarDict n m++runRVar :: RandomSource m s => RVar a -> s -> m a+runRVar = runRVarT++-- |A random variable with access to operations in an underlying monad.  Useful+-- examples include any form of state for implementing random processes with hysteresis,+-- or writer monads for implementing tracing of complicated algorithms.+newtype RVarT n a = RVarT { unRVarT :: forall m r. (a -> m r) -> RVarDict n m -> m r }++-- | \"Runs\" the monad.+runRVarT :: (Lift n m, RandomSource m s) => RVarT n a -> s -> m a+runRVarT (RVarT m) (src) = m return (RVarDict src)++instance Functor (RVarT n) where     fmap = liftM -instance Monad RVar where-    return x = RVar (\_ -> return x)-    fail s   = RVar (\_ -> fail s)-    (RVar x) >>= f = RVar (\s -> do-            x <- x s-            case f x of-                RVar y -> y s-        )+instance Monad (RVarT n) where+    return x = RVarT $ \k _ -> k x+    fail s   = RVarT $ \_ (RVarDict _) -> fail s+    (RVarT m) >>= k = RVarT $ \c s -> m (\a -> unRVarT (k a) c s) s -instance Applicative RVar where+instance Applicative (RVarT n) where     pure  = return     (<*>) = ap -instance Distribution RVar a where-    rvar = id-    sampleFrom src x = runDistM x src+instance T.MonadTrans RVarT where+    lift m = RVarT $ \k r@(RVarDict _) -> L.lift m >>= \a -> k a -instance MonadRandom RVar where-    getRandomBytes n = RVar (\s -> getRandomBytesFrom s n)-    getRandomWords n = RVar (\s -> getRandomWordsFrom s n)+instance Lift (RVarT Identity) (RVarT m) where+    lift (RVarT m) = RVarT $ \k (RVarDict src) -> m k (RVarDict src) --- some 'fundamental' RVars+instance MonadIO m => MonadIO (RVarT m) where+    liftIO = T.lift . liftIO++instance MonadRandom (RVarT n) where+    getRandomByte   = RVarT $ \k (RVarDict s) -> getRandomByteFrom   s >>= \a -> k a+    getRandomWord   = RVarT $ \k (RVarDict s) -> getRandomWordFrom   s >>= \a -> k a+    getRandomDouble = RVarT $ \k (RVarDict s) -> getRandomDoubleFrom s >>= \a -> k a++-- some 'fundamental' RVarTs -- this maybe ought to even be a part of the RandomSource class...+{-# INLINE nByteInteger #-} -- |A random variable evenly distributed over all unsigned integers from -- 0 to 2^(8*n)-1, inclusive.-nByteInteger :: Int -> RVar Integer-nByteInteger n-    | n .&. 7 == 0-    = do-        xs <- getRandomWords (n `shiftR` 3)-        return $! concatWords xs-    | n > 8-    = do-        let nWords = n `shiftR` 3-            nBytes = n .&. 7-        ws <- getRandomWords nWords-        bs <- getRandomBytes nBytes-        return $! ((concatWords ws `shiftL` (nBytes `shiftL` 3)) .|. concatBytes bs)-    | otherwise-    = do-        xs <- getRandomBytes n-        return $! concatBytes xs+nByteInteger :: Int -> RVarT m Integer+nByteInteger 1 = do+    x <- getRandomByte+    return $! toInteger x+nByteInteger 8 = do+    x <- getRandomWord+    return $! toInteger x+nByteInteger (n+8) = do+    x <- getRandomWord+    y <- nByteInteger n+    return $! ((toInteger x `shiftL` n) .|. y)+nByteInteger n = do+    x <- getRandomWord+    return $! toInteger (x `shiftR` ((8-n) `shiftL` 3)) +{-# INLINE nBitInteger #-} -- |A random variable evenly distributed over all unsigned integers from -- 0 to 2^n-1, inclusive.-nBitInteger :: Int -> RVar Integer-nBitInteger n-    | n .&. 7 == 0-    = nByteInteger (n `shiftR` 3)-    | otherwise-    = do-        x <- nByteInteger ((n `shiftR` 3) + 1)-        return $! (x .&. (bit n - 1))--concatBytes :: (Bits a, Num a) => [Word8] -> a-concatBytes = concatBits fromIntegral--concatWords :: (Bits a, Num a) => [Word64] -> a-concatWords = concatBits fromIntegral--concatBits :: (Bits a, Bits b, Num b) => (a -> b) -> [a] -> b-concatBits f [] = 0-concatBits f (x:xs) = f x .|. (concatBits f xs `shiftL` bitSize x)+nBitInteger :: Int -> RVarT m Integer+nBitInteger 8  = do+    x <- getRandomByte+    return $! toInteger x+nBitInteger (n+64) = do+    x <- getRandomWord+    y <- nBitInteger n+    return $! (toInteger x `shiftL` n) .|. y+nBitInteger n = do+        x <- getRandomWord+        return $! toInteger (x `shiftR` (64-n))

@@ -1,16 +0,0 @@-{-- -      ``Data/Random/RVar''- -}-{-# LANGUAGE-    MultiParamTypeClasses,-    FlexibleInstances-  #-}--module Data.Random.RVar where--import Data.Random.Source-import {-# SOURCE #-} Data.Random.Distribution--data RVar a-instance MonadRandom RVar-instance Distribution RVar a

@@ -0,0 +1,35 @@+{-+ -      ``Data/Random/Sample''+ -}+{-# LANGUAGE+        MultiParamTypeClasses,+        FlexibleInstances, FlexibleContexts, +        IncoherentInstances+  #-}++module Data.Random.Sample where++import Data.Random.Distribution+import Data.Random.Lift+import Data.Random.RVar+import Data.Random.Source+import Data.Random.Source.Std++-- |A typeclass allowing 'Distribution's and 'RVar's to be sampled.  Both may+-- also be sampled via 'runRVar' or 'runRVarT', but I find it psychologically+-- pleasing to be able to sample both using this function.+class Sampleable d m t where+    -- |Directly sample from a distribution or random variable, using the given source of entropy.+    sampleFrom :: RandomSource m s => s -> d t -> m t++instance Distribution d t => Sampleable d m t where+    sampleFrom src d = runRVarT (rvar d) src++-- This instance conflicts with the other, but because RVarT is not a Distribution there is no conflict.+instance Lift m n => Sampleable (RVarT m) n t where+    sampleFrom src x = runRVarT x src++-- |Sample a distribution using the default source of entropy for the+-- monad in which the sampling occurs.+sample :: (Sampleable d m t, MonadRandom m) => d t -> m t+sample = sampleFrom StdRandom

@@ -13,6 +13,7 @@ import Data.Word import Data.Bits import Data.List+import Control.Monad  import Data.Random.Internal.Words @@ -22,54 +23,68 @@ -- when directly requesting entropy for a random variable these functions -- are used. -- --- The minimal definition is either 'getRandomBytes' or 'getRandomWords'.+-- The minimal definition is either 'getRandomByte' or 'getRandomWord'.+-- 'getRandomDouble' is defaulted in terms of 'getRandomWord'. class Monad m => MonadRandom m where-    -- |get the specified number of random (uniformly distributed) bytes-    getRandomBytes :: Int -> m [Word8]-    getRandomBytes n-        | n .&. 7 == 0-        = do-            let wc = n `shiftR` 3-            ws <- getRandomWords wc-            return (concatMap wordToBytes ws)-        | otherwise-        = do-            let wc = (n `shiftR` 3) + 1-            ws <- getRandomWords wc-            return . take n . concatMap wordToBytes $ ws+    -- |Get a random uniformly-distributed byte.+    getRandomByte :: m Word8+    getRandomByte = do+        word <- getRandomWord+        return (fromIntegral word)+    +    -- |Get a random 'Word64' uniformly-distributed over the full range of the type.+    getRandomWord :: m Word64+    getRandomWord = do+        b0 <- getRandomByte+        b1 <- getRandomByte+        b2 <- getRandomByte+        b3 <- getRandomByte+        b4 <- getRandomByte+        b5 <- getRandomByte+        b6 <- getRandomByte+        b7 <- getRandomByte         -    -- |alternate basis function, providing access to larger chunks-    getRandomWords :: Int -> m [Word64]-    getRandomWords n = do-        bs <- getRandomBytes (n `shiftL` 3)-        return (bytesToWords bs)+        return (buildWord b0 b1 b2 b3 b4 b5 b6 b7)+    +    -- |Get a random 'Double' uniformly-distributed over the interval [0,1)+    getRandomDouble :: m Double+    getRandomDouble = do+        word <- getRandomWord+        return (wordToDouble word)  -- |A source of entropy which can be used in the given monad. ----- The minimal definition is either 'getRandomBytesFrom' or 'getRandomWordsFrom'+-- The minimal definition is either 'getRandomByteFrom' or 'getRandomWordFrom'.+-- 'getRandomDoubleFrom' is defaulted in terms of 'getRandomWordFrom' class Monad m => RandomSource m s where-    getRandomBytesFrom :: s -> Int -> m [Word8]-    getRandomBytesFrom src n-        | n .&. 7 == 0-        = do-            let wc = n `shiftR` 3-            ws <- getRandomWordsFrom src wc-            return (concatMap wordToBytes ws)-        | otherwise-        = do-            let wc = (n `shiftR` 3) + 1-            ws <- getRandomWordsFrom src wc-            return . take n . concatMap wordToBytes $ ws+    -- |Get a random uniformly-distributed byte.+    getRandomByteFrom :: s -> m Word8+    getRandomByteFrom src = do+        word <- getRandomWordFrom src+        return (fromIntegral word)+    +    -- |Get a random 'Word64' uniformly-distributed over the full range of the type.+    getRandomWordFrom :: s -> m Word64+    getRandomWordFrom src = do+        b0 <- getRandomByteFrom src+        b1 <- getRandomByteFrom src+        b2 <- getRandomByteFrom src+        b3 <- getRandomByteFrom src+        b4 <- getRandomByteFrom src+        b5 <- getRandomByteFrom src+        b6 <- getRandomByteFrom src+        b7 <- getRandomByteFrom src         +        return (buildWord b0 b1 b2 b3 b4 b5 b6 b7)     -    getRandomWordsFrom :: s -> Int -> m [Word64]-    getRandomWordsFrom src n = do-        bs <- getRandomBytesFrom src (n `shiftL` 3)-        return (bytesToWords bs)--instance Monad m => RandomSource m (Int -> m [Word8]) where-    getRandomBytesFrom = id+    -- |Get a random 'Double' uniformly-distributed over the interval [0,1)+    getRandomDoubleFrom :: s -> m Double+    getRandomDoubleFrom src = do+        word <- getRandomWordFrom src+        return (wordToDouble word) -instance Monad m => RandomSource m (Int -> m [Word64]) where-    getRandomWordsFrom = id+instance Monad m => RandomSource m (m Word8) where+    getRandomByteFrom = id +instance Monad m => RandomSource m (m Word64) where+    getRandomWordFrom = id

@@ -5,20 +5,37 @@     MultiParamTypeClasses   #-} -module Data.Random.Source.DevRandom where+module Data.Random.Source.DevRandom +    ( DevRandom(..)+    ) where  import Data.Random.Source  import GHC.IOBase (unsafePerformIO)-import Data.ByteString (hGet, unpack)-import System.IO (openBinaryFile, IOMode(..))+import System.IO (openBinaryFile, hGetBuf, IOMode(..))+import Foreign  -- |On systems that have it, \/dev\/random is a handy-dandy ready-to-use source--- of nonsense.-data DevRandom = DevRandom-{-# NOINLINE devRandom #-}-devRandom = unsafePerformIO (openBinaryFile "/dev/random" ReadMode)+-- of nonsense.  Keep in mind that on some systems, Linux included, \/dev\/random+-- collects \"real\" entropy, and if you don't have a good source of it, such as+-- special hardware for the purpose or a *lot* of network traffic, it's pretty easy+-- to suck the entropy pool dry with entropy-intensive applications.  For many+-- purposes other than cryptography, \/dev\/urandom is preferable because when it+-- runs out of real entropy it'll still churn out pseudorandom data.+data DevRandom = DevRandom | DevURandom -instance RandomSource IO DevRandom where-    getRandomBytesFrom DevRandom n = fmap unpack (hGet devRandom n)+{-# NOINLINE devRandom  #-}+devRandom  = unsafePerformIO (openBinaryFile "/dev/random"  ReadMode)+{-# NOINLINE devURandom #-}+devURandom = unsafePerformIO (openBinaryFile "/dev/urandom" ReadMode) +dev DevRandom  = devRandom+dev DevURandom = devURandom++instance RandomSource IO DevRandom where+    getRandomByteFrom src  = allocaBytes 1 $ \buf -> do+        1 <- hGetBuf (dev src) buf  1+        peek buf+    getRandomWordFrom src  = allocaBytes 8 $ \buf -> do+        8 <- hGetBuf (dev src) buf  8+        peek (castPtr buf)

@@ -3,11 +3,13 @@  -} {-# LANGUAGE     MultiParamTypeClasses,-    FlexibleContexts, FlexibleInstances+    FlexibleContexts, FlexibleInstances,+    UndecidableInstances   #-}  module Data.Random.Source.PureMT where +import Data.Random.Internal.Words import Data.Random.Source import System.Random.Mersenne.Pure64 @@ -15,41 +17,112 @@ import Data.Word  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 make a -- 'RandomSource' usable in any monad in which the reference can be modified. ----- For example, if @x :: TVar PureMT@, @getRandomWordsFromMTRef x@ can be+-- For example, if @x :: TVar PureMT@, @getRandomWordFromMTRef x@ can be -- used as a 'RandomSource' in 'IO', 'STM', or any monad which is an instance--- of 'MonadIO'.-getRandomWordsFromMTRef :: ModifyRef sr m PureMT => sr -> Int -> m [Word64]-getRandomWordsFromMTRef ref n = do-    atomicModifyRef ref (randomWords n [])+-- of 'MonadIO'.  These functions can also be used to implement additional+-- 'RandomSource' instances for mutable references to 'PureMT' states.+getRandomWordFromMTRef :: ModifyRef sr m PureMT => sr -> m Word64+getRandomWordFromMTRef ref = do+    atomicModifyRef ref (swap . randomWord64)          where         swap (a,b) = (b,a)-        randomWords    0  ws mt = (mt, ws)-        randomWords (n+1) ws mt = case randomWord64 mt of-            (w, mt) -> randomWords n (w:ws) mt --- |Similarly, @getRandomWordsFromMTState x@ can be used in any \"state\"+getRandomByteFromMTRef :: ModifyRef sr m PureMT => sr -> m Word8+getRandomByteFromMTRef ref = do+    x <- atomicModifyRef ref (swap . randomInt)+    return (fromIntegral x)+    +    where+        swap (a,b) = (b,a)++-- for whatever reason, my simple wordToDouble is faster than whatever+-- the mersenne random library is using, at least in the version I have.+-- if this changes, switch to the commented version.+-- Same thing below, in getRandomDoubleFromMTState.+getRandomDoubleFromMTRef :: ModifyRef sr m PureMT => sr -> m Double+getRandomDoubleFromMTRef src = liftM wordToDouble (getRandomWordFromMTRef src)+-- getRandomDoubleFromMTRef ref = do+--     atomicModifyRef ref (swap . randomDouble)+--     +--     where+--         swap (a,b) = (b,a)++-- |Similarly, @getRandomWordFromMTState x@ can be used in any \"state\" -- monad in the mtl sense whose state is a 'PureMT' generator. -- 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.-getRandomWordsFromMTState :: MonadState PureMT m => Int -> m [Word64]-getRandomWordsFromMTState n = do+-- other 'Distribution' instances to \"pure\" random variables (e.g., by+-- @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@).+getRandomWordFromMTState :: MonadState PureMT m => m Word64+getRandomWordFromMTState = do     mt <- get-    let randomWords    0  ws mt = (mt, ws)-        randomWords (n+1) ws mt = case randomWord64 mt of-            (w, mt) -> randomWords n (w:ws) mt-        -        (newMt, ws) = randomWords n [] mt+    let (ws, newMt) = randomWord64 mt     put newMt     return ws +getRandomByteFromMTState :: MonadState PureMT m => m Word8+getRandomByteFromMTState = do+    mt <- get+    let (ws, newMt) = randomInt mt+    put newMt+    return (fromIntegral ws)+++getRandomDoubleFromMTState :: MonadState PureMT m => m Double+getRandomDoubleFromMTState = liftM wordToDouble getRandomWordFromMTState+-- getRandomDoubleFromMTState = do+--     mt <- get+--     let (x, newMt) = randomDouble mt+--     put newMt+--     return x+ instance MonadRandom (State PureMT) where-    getRandomWords = getRandomWordsFromMTState+    getRandomByte   = getRandomByteFromMTState+    getRandomWord   = getRandomWordFromMTState+    getRandomDouble = getRandomDoubleFromMTState +instance MonadRandom (S.State PureMT) where+    getRandomByte   = getRandomByteFromMTState+    getRandomWord   = getRandomWordFromMTState+    getRandomDouble = getRandomDoubleFromMTState+ instance Monad m => MonadRandom (StateT PureMT m) where-    getRandomWords = getRandomWordsFromMTState+    getRandomByte   = getRandomByteFromMTState+    getRandomWord   = getRandomWordFromMTState+    getRandomDouble = getRandomDoubleFromMTState++instance Monad m => MonadRandom (S.StateT PureMT m) where+    getRandomByte   = getRandomByteFromMTState+    getRandomWord   = getRandomWordFromMTState+    getRandomDouble = getRandomDoubleFromMTState+++instance (ModifyRef (IORef PureMT) m PureMT) => RandomSource m (IORef PureMT) where+    {-# SPECIALIZE instance RandomSource IO (IORef PureMT)#-}+    getRandomByteFrom   = getRandomByteFromMTRef+    getRandomWordFrom   = getRandomWordFromMTRef+    getRandomDoubleFrom = getRandomDoubleFromMTRef++instance (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) #-}+    getRandomByteFrom   = getRandomByteFromMTRef+    getRandomWordFrom   = getRandomWordFromMTRef+    getRandomDoubleFrom = getRandomDoubleFromMTRef++instance (ModifyRef (TVar PureMT) m PureMT) => RandomSource m (TVar PureMT) where+    {-# SPECIALIZE instance RandomSource IO  (TVar PureMT) #-}+    {-# SPECIALIZE instance RandomSource STM (TVar PureMT) #-}+    getRandomByteFrom   = getRandomByteFromMTRef+    getRandomWordFrom   = getRandomWordFromMTRef+    getRandomDoubleFrom = getRandomDoubleFromMTRef+

@@ -16,5 +16,6 @@ data StdRandom = StdRandom  instance MonadRandom m => RandomSource m StdRandom where-    getRandomBytesFrom StdRandom = getRandomBytes-    getRandomWordsFrom StdRandom = getRandomWords+    getRandomByteFrom   StdRandom = getRandomByte+    getRandomWordFrom   StdRandom = getRandomWord+    getRandomDoubleFrom StdRandom = getRandomDouble

@@ -7,82 +7,128 @@  module Data.Random.Source.StdGen where +import Data.Random.Internal.Words import Data.Random.Source import System.Random import Control.Monad import Control.Monad.State+import qualified Control.Monad.ST.Strict as S+import qualified Control.Monad.State.Strict as S import Data.StateRef import Data.Word  instance (ModifyRef (IORef   StdGen) m StdGen) => RandomSource m (IORef   StdGen) where-    getRandomBytesFrom = getRandomBytesFromRandomGenRef-    getRandomWordsFrom = getRandomWordsFromRandomGenRef+    {-# SPECIALIZE instance RandomSource IO (IORef StdGen) #-}+    getRandomByteFrom   = getRandomByteFromRandomGenRef+    getRandomWordFrom   = getRandomWordFromRandomGenRef+    getRandomDoubleFrom = getRandomDoubleFromRandomGenRef instance (ModifyRef (TVar    StdGen) m StdGen) => RandomSource m (TVar    StdGen) where-    getRandomBytesFrom = getRandomBytesFromRandomGenRef-    getRandomWordsFrom = getRandomWordsFromRandomGenRef+    {-# SPECIALIZE instance RandomSource IO  (TVar StdGen) #-}+    {-# SPECIALIZE instance RandomSource STM (TVar StdGen) #-}+    getRandomByteFrom   = getRandomByteFromRandomGenRef+    getRandomWordFrom   = getRandomWordFromRandomGenRef+    getRandomDoubleFrom = getRandomDoubleFromRandomGenRef instance (ModifyRef (STRef s StdGen) m StdGen) => RandomSource m (STRef s StdGen) where-    getRandomBytesFrom = getRandomBytesFromRandomGenRef-    getRandomWordsFrom = getRandomWordsFromRandomGenRef+    {-# SPECIALIZE instance RandomSource (ST s) (STRef s StdGen) #-}+    {-# SPECIALIZE instance RandomSource (S.ST s) (STRef s StdGen) #-}+    getRandomByteFrom   = getRandomByteFromRandomGenRef+    getRandomWordFrom   = getRandomWordFromRandomGenRef+    getRandomDoubleFrom = getRandomDoubleFromRandomGenRef -getRandomBytesFromStdGenIO :: Int -> IO [Word8]-getRandomBytesFromStdGenIO n = do-    ints <- replicateM n (randomRIO (0, 255))-    let bytes = map fromIntegral (ints :: [Int])-    return bytes+getRandomByteFromStdGenIO :: IO Word8+getRandomByteFromStdGenIO = do+    int <- randomRIO (0, 255) :: IO Int+    return (fromIntegral int) +getRandomWordFromStdGenIO :: IO Word64+getRandomWordFromStdGenIO = do+    int <- randomRIO (0, 0xffffffffffffffff)+    return (fromInteger int)++-- based on reading the source of the "random" library's implementation, I do+-- not believe that the randomRIO (0,1) implementation for Double is capable of producing+-- the value 0.  Therefore, I'm not using it.  If this is an incorrect reading on+-- my part, or if this changes, then feel free to use the commented version.+-- Same goes for the other getRandomDouble... functions here.+getRandomDoubleFromStdGenIO :: IO Double+getRandomDoubleFromStdGenIO = liftM wordToDouble getRandomWordFromStdGenIO+-- getRandomDoubleFromStdGenIO = randomRIO (0, 1)+ -- |Given a mutable reference to a 'RandomGen' generator, we can make a -- 'RandomSource' usable in any monad in which the reference can be modified. ----- For example, if @x :: TVar StdGen@, @getRandomBytesFromRandomGenRef x@ can be+-- For example, if @x :: TVar StdGen@, @getRandomByteFromRandomGenRef x@ can be -- used as a 'RandomSource' in 'IO', 'STM', or any monad which is an instance -- of 'MonadIO'.  It's generally probably better to use--- 'getRandomWordsFromRandomGenRef' though, as this one is likely to throw--- away a lot of perfectly good entropy.-getRandomBytesFromRandomGenRef :: (ModifyRef sr m g, RandomGen g) =>-                                  sr -> Int -> m [Word8]-getRandomBytesFromRandomGenRef g n = do-    let swap (a,b) = (b,a)-    ints <- replicateM n (atomicModifyRef g (swap . randomR (0, 255)))-    let bytes = map fromIntegral (ints :: [Int])-    return bytes-    --- |Similarly, @getRandomWordsFromRandomGenState x@ can be used in any \"state\"+-- 'getRandomWordFromRandomGenRef' though, as this one is likely to throw+-- away a lot of perfectly good entropy.  Better still is to use these 3 functions+-- together to create a 'RandomSource' instance for the reference you're using,+-- if one does not already exist.+getRandomByteFromRandomGenRef :: (ModifyRef sr m g, RandomGen g) =>+                                  sr -> m Word8+getRandomByteFromRandomGenRef g = atomicModifyRef g (swap . randomR (0,255))+    where +        swap :: (Int, a) -> (a, Word8)+        swap (a,b) = (b,fromIntegral a)++getRandomWordFromRandomGenRef :: (ModifyRef sr m g, RandomGen g) =>+                                  sr -> m Word64+getRandomWordFromRandomGenRef g = atomicModifyRef g (swap . randomR (0,0xffffffffffffffff))+    where swap (a,b) = (b,fromInteger a)++getRandomDoubleFromRandomGenRef :: (ModifyRef sr m g, RandomGen g) =>+                                  sr -> m Double+getRandomDoubleFromRandomGenRef g = liftM wordToDouble (getRandomWordFromRandomGenRef g)+-- getRandomDoubleFromRandomGenRef g = atomicModifyRef g (swap . randomR (0,1))+--     where swap (a,b) = (b,a)++-- |Similarly, @getRandomWordFromRandomGenState x@ can be used in any \"state\" -- monad in the mtl sense whose state is a 'RandomGen' generator. -- 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.-getRandomBytesFromRandomGenState :: (RandomGen g, MonadState g m) =>-                                  Int -> m [Word8]-getRandomBytesFromRandomGenState n = replicateM n $ do+getRandomByteFromRandomGenState :: (RandomGen g, MonadState g m) => m Word8+getRandomByteFromRandomGenState = do     g <- get-    case randomR (0,255 :: Int) g of+    case randomR (0, 255 :: Int) g of         (i,g) -> do             put g             return (fromIntegral i) --- |See 'getRandomBytesFromRandomGenRef'-getRandomWordsFromRandomGenRef :: (ModifyRef sr m g, RandomGen g) =>-                                  sr -> Int -> m [Word64]-getRandomWordsFromRandomGenRef g n = do-    let swap (a,b) = (b,a)-    ints <- replicateM n (atomicModifyRef g (swap . randomR (0, 2^64-1)))-    let bytes = map fromInteger ints-    return bytes-    --- |See 'getRandomBytesFromRandomGenState'-getRandomWordsFromRandomGenState :: (RandomGen g, MonadState g m) =>-                                  Int -> m [Word64]-getRandomWordsFromRandomGenState n = replicateM n $ do+getRandomWordFromRandomGenState :: (RandomGen g, MonadState g m) => m Word64+getRandomWordFromRandomGenState = do     g <- get-    case randomR (0,2^64-1) g of+    case randomR (0, 0xffffffffffffffff) g of         (i,g) -> do             put g             return (fromInteger i) +getRandomDoubleFromRandomGenState :: (RandomGen g, MonadState g m) => m Double+getRandomDoubleFromRandomGenState = liftM wordToDouble getRandomWordFromRandomGenState+-- getRandomDoubleFromRandomGenState = do+--     g <- get+--     case randomR (0, 1) g of+--         (x,g) -> do+--             put g+--             return x++ instance MonadRandom (State StdGen) where-    getRandomBytes = getRandomBytesFromRandomGenState-    getRandomWords = getRandomWordsFromRandomGenState+    getRandomByte   = getRandomByteFromRandomGenState+    getRandomWord   = getRandomWordFromRandomGenState+    getRandomDouble = getRandomDoubleFromRandomGenState  instance Monad m => MonadRandom (StateT StdGen m) where-    getRandomBytes = getRandomBytesFromRandomGenState-    getRandomWords = getRandomWordsFromRandomGenState+    getRandomByte   = getRandomByteFromRandomGenState+    getRandomWord   = getRandomWordFromRandomGenState+    getRandomDouble = getRandomDoubleFromRandomGenState++instance MonadRandom (S.State StdGen) where+    getRandomByte   = getRandomByteFromRandomGenState+    getRandomWord   = getRandomWordFromRandomGenState+    getRandomDouble = getRandomDoubleFromRandomGenState++instance Monad m => MonadRandom (S.StateT StdGen m) where+    getRandomByte   = getRandomByteFromRandomGenState+    getRandomWord   = getRandomWordFromRandomGenState+    getRandomDouble = getRandomDoubleFromRandomGenState