diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,5 @@
+#!/usr/bin/env runhaskell
+
+> import Distribution.Simple
+> main = defaultMain
+
diff --git a/random-fu.cabal b/random-fu.cabal
new file mode 100644
--- /dev/null
+++ b/random-fu.cabal
@@ -0,0 +1,52 @@
+name:                   random-fu
+version:                0.0.0.2
+stability:              experimental
+
+cabal-version:          >= 1.2
+build-type:             Simple
+
+author:                 James Cook <james.cook@usma.edu>
+maintainer:             James Cook <james.cook@usma.edu>
+license:                PublicDomain
+homepage:               http://code.haskell.org/~mokus/random-fu
+
+category:               Math
+synopsis:               Random number generation
+description:            Random number generation based on orthogonal typeclasses
+                        for entropy sources and random variable distributions, all
+                        served up on a monadic platter.  Aspires to be useful
+                        in an idiomatic way in both \"pure\" and \"impure\" styles,
+                        as well as reasonably fast.  May not yet meet the latter
+                        goal, but I think the former is starting to shape up
+                        nicely.
+
+Library
+  hs-source-dirs:       src
+  exposed-modules:      Data.Random
+                        Data.Random.Distribution
+                        Data.Random.Distribution.Bernoulli
+                        Data.Random.Distribution.Beta
+                        Data.Random.Distribution.Binomial
+                        Data.Random.Distribution.Exponential
+                        Data.Random.Distribution.Discrete
+                        Data.Random.Distribution.Gamma
+                        Data.Random.Distribution.Normal
+                        Data.Random.Distribution.Poisson
+                        Data.Random.Distribution.Triangular
+                        Data.Random.Distribution.Uniform
+                        Data.Random.Internal.Classification
+                        Data.Random.Internal.Words
+                        Data.Random.RVar
+                        Data.Random.Source
+                        Data.Random.Source.DevRandom
+                        Data.Random.Source.StdGen
+                        Data.Random.Source.PureMT
+                        Data.Random.Source.Std
+                        
+  build-depends:        base >= 3,
+                        bytestring,
+                        mersenne-random-pure64,
+                        monad-loops >= 0.3.0.1,
+                        mtl,
+                        random,
+                        stateref
diff --git a/src/Data/Random.hs b/src/Data/Random.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random.hs
@@ -0,0 +1,58 @@
+{-
+ -      ``Data/Random''
+ -}
+{-# LANGUAGE
+    FlexibleContexts
+  #-}
+
+-- |Random numbers and stuff...
+-- 
+-- Data.Random.Source exports the typeclasses for entropy sources, and
+-- Data.Random.Source.* export various instances and/or functions with which
+-- instances can be defined.
+-- 
+-- Data.Random.Distribution exports the typeclasses for sampling distributions,
+-- and Data.Random.Distribution.* export various specific distributions.
+--
+-- Data.Random.RVar exports the 'RVar' type, which is a probability distribution
+-- monad that allows for concise definitions of random variables, as well as
+-- a couple handy 'RVar's.
+
+module Data.Random
+    ( module Data.Random.Source
+    , module Data.Random.Source.DevRandom
+    , module Data.Random.Source.StdGen
+    , module Data.Random.Source.PureMT
+    , module Data.Random.Source.Std
+    , module Data.Random.Distribution
+    , module Data.Random.Distribution.Bernoulli
+    , module Data.Random.Distribution.Beta
+    , module Data.Random.Distribution.Binomial
+    , module Data.Random.Distribution.Discrete
+    , module Data.Random.Distribution.Gamma
+    , module Data.Random.Distribution.Exponential
+    , module Data.Random.Distribution.Normal
+    , module Data.Random.Distribution.Poisson
+    , module Data.Random.Distribution.Triangular
+    , module Data.Random.Distribution.Uniform
+    , module Data.Random.RVar
+    ) where
+
+import Data.Random.Source
+import Data.Random.Source.DevRandom
+import Data.Random.Source.StdGen
+import Data.Random.Source.PureMT
+import Data.Random.Source.Std
+import Data.Random.Distribution
+import Data.Random.Distribution.Bernoulli
+import Data.Random.Distribution.Beta
+import Data.Random.Distribution.Binomial
+import Data.Random.Distribution.Discrete
+import Data.Random.Distribution.Gamma
+import Data.Random.Distribution.Exponential
+import Data.Random.Distribution.Normal
+import Data.Random.Distribution.Poisson
+import Data.Random.Distribution.Triangular
+import Data.Random.Distribution.Uniform
+import Data.Random.RVar
+
diff --git a/src/Data/Random/Distribution.hs b/src/Data/Random/Distribution.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution.hs
@@ -0,0 +1,32 @@
+{-
+ -      ``Data/Random/Distribution''
+ -}
+{-# LANGUAGE
+    MultiParamTypeClasses, FlexibleContexts
+  #-}
+
+module Data.Random.Distribution where
+
+import {-# SOURCE #-} 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'.
+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)
+
+-- |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
diff --git a/src/Data/Random/Distribution.hs-boot b/src/Data/Random/Distribution.hs-boot
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution.hs-boot
@@ -0,0 +1,10 @@
+{-
+ -      ``Data/Random/Distribution''
+ -}
+{-# LANGUAGE
+    MultiParamTypeClasses, KindSignatures
+  #-}
+
+module Data.Random.Distribution where
+
+class Distribution (d :: * -> *) t
diff --git a/src/Data/Random/Distribution/Bernoulli.hs b/src/Data/Random/Distribution/Bernoulli.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution/Bernoulli.hs
@@ -0,0 +1,47 @@
+{-
+ -      ``Data/Random/Distribution/Bernoulli''
+ -}
+{-# LANGUAGE
+    MultiParamTypeClasses,
+    FlexibleInstances, FlexibleContexts,
+    UndecidableInstances
+  #-}
+
+module Data.Random.Distribution.Bernoulli where
+
+import Data.Random.Internal.Classification
+
+import Data.Random.Source
+import Data.Random.Distribution
+import Data.Random.RVar
+
+import Data.Random.Distribution.Uniform
+
+import Data.Int
+import Data.Word
+
+bernoulli :: (Distribution (Bernoulli b) a) => b -> RVar a
+bernoulli p = sample (Bernoulli p)
+
+boolBernoulli p = do
+    x <- realFloatUniform 0 1
+    return (x <= p)
+
+generalBernoulli t f p = do
+    x <- boolBernoulli p
+    return (if x then t else f)
+
+class (Classification NumericType t c) => BernoulliByClassification c t where
+    bernoulliByClassification :: RealFloat a => a -> RVar t
+
+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)
+
+data Bernoulli b a = Bernoulli b
+
+instance (BernoulliByClassification c t, RealFloat b) => Distribution (Bernoulli b) t where
+    rvar (Bernoulli p) = bernoulliByClassification p
diff --git a/src/Data/Random/Distribution/Beta.hs b/src/Data/Random/Distribution/Beta.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution/Beta.hs
@@ -0,0 +1,39 @@
+{-
+ -      ``Data/Random/Distribution/Beta''
+ -}
+{-# LANGUAGE
+    MultiParamTypeClasses,
+    FlexibleInstances, FlexibleContexts,
+    UndecidableInstances
+  #-}
+
+module Data.Random.Distribution.Beta where
+
+import Data.Random.Source
+import Data.Random.RVar
+import Data.Random.Distribution
+import Data.Random.Distribution.Gamma
+import Data.Random.Distribution.Uniform
+
+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
+    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
+    return (x / (x + y))
+
+beta :: Distribution Beta a => a -> a -> RVar a
+beta a b = sample (Beta a b)
+
+data Beta a = Beta a a
+
+instance (RealFloat a) => Distribution Beta a where
+    rvar (Beta a b) = realFloatBeta a b
diff --git a/src/Data/Random/Distribution/Binomial.hs b/src/Data/Random/Distribution/Binomial.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution/Binomial.hs
@@ -0,0 +1,64 @@
+{-
+ -      ``Data/Random/Distribution/Binomial''
+ -}
+{-# LANGUAGE
+    MultiParamTypeClasses,
+    FlexibleInstances, FlexibleContexts,
+    UndecidableInstances
+  #-}
+
+module Data.Random.Distribution.Binomial where
+
+import Data.Random.Internal.Classification
+
+import Data.Random.Source
+import Data.Random.Distribution
+import Data.Random.RVar
+
+import Data.Random.Distribution.Beta
+import Data.Random.Distribution.Uniform
+
+import Data.Int
+import Data.Word
+import Control.Monad
+
+    -- 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
+integralBinomial t p = bin 0 t p
+    where
+        bin k t p
+            | t > 10    = do
+                let a = 1 + t `div` 2
+                    b = 1 + t - a
+        
+                x <- realFloatBetaFromIntegral a b
+                if x >= p
+                    then bin  k      (a - 1) (p / x)
+                    else bin (k + a) (b - 1) ((p - x) / (1 - x))
+        
+            | otherwise = count k t
+                where
+                    count k  0    = return k
+                    count k (n+1) = do
+                        x <- realFloatStdUniform
+                        (count $! (if x < p then k + 1 else k)) n
+
+
+
+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
+
+data Binomial b a = Binomial a b
+
diff --git a/src/Data/Random/Distribution/Discrete.hs b/src/Data/Random/Distribution/Discrete.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution/Discrete.hs
@@ -0,0 +1,35 @@
+{-
+ -      ``Data/Random/Distribution/Discrete''
+ -}
+{-# LANGUAGE
+    MultiParamTypeClasses,
+    FlexibleInstances, FlexibleContexts
+  #-}
+
+module Data.Random.Distribution.Discrete where
+
+import Data.Random.RVar
+import Data.Random.Distribution
+import Data.Random.Distribution.Uniform
+
+import Control.Monad
+
+discrete :: Distribution (Discrete p) a => [(p,a)] -> RVar a
+discrete ps = rvar (Discrete ps)
+
+data Discrete p a = Discrete [(p, a)]
+
+instance (Num p, Ord p, Distribution Uniform p) => Distribution (Discrete p) a where
+    rvar (Discrete []) = fail "discrete distribution over empty set cannot be sampled"
+    rvar (Discrete ds) = do
+        let (ps, xs) = unzip ds
+            cs = scanl1 (+) ps
+        
+        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
+            ]
diff --git a/src/Data/Random/Distribution/Exponential.hs b/src/Data/Random/Distribution/Exponential.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution/Exponential.hs
@@ -0,0 +1,28 @@
+{-
+ -      ``Data/Random/Distribution/Exponential''
+ -}
+{-# LANGUAGE
+    MultiParamTypeClasses,
+    FlexibleInstances, FlexibleContexts,
+    UndecidableInstances
+  #-}
+
+module Data.Random.Distribution.Exponential where
+
+import Data.Random.Source
+import Data.Random.RVar
+import Data.Random.Distribution
+import Data.Random.Distribution.Uniform
+
+data Exponential a = Exp a
+
+realFloatExponential :: RealFloat a => a -> RVar a
+realFloatExponential lambdaRecip = do
+    x <- realFloatStdUniform
+    return (negate (log x) * lambdaRecip)
+
+exponential :: Distribution Exponential a => a -> RVar a
+exponential = sample . Exp
+
+instance (RealFloat a) => Distribution Exponential a where
+    rvar (Exp lambdaRecip) = realFloatExponential lambdaRecip
diff --git a/src/Data/Random/Distribution/Gamma.hs b/src/Data/Random/Distribution/Gamma.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution/Gamma.hs
@@ -0,0 +1,74 @@
+{-
+ -      ``Data/Random/Distribution/Gamma''
+ -
+ -  needs cleanup, verification, and automagic selection of appropriate
+ -  algorithms, and proper citations.
+ -
+ -  should eliminate spurious 'border crossings' betweer RVars and sampleFrom
+ -  perhaps Distribution class should have as its basis a function of type
+ -  (d t -> RVar t)
+ -}
+{-# LANGUAGE
+    MultiParamTypeClasses,
+    FlexibleInstances, FlexibleContexts,
+    UndecidableInstances
+  #-}
+
+module Data.Random.Distribution.Gamma where
+
+import Data.Random.Source
+import Data.Random.RVar
+import Data.Random.Distribution
+import Data.Random.Distribution.Uniform
+import Data.Random.Distribution.Normal
+
+import Control.Monad
+
+    -- translated from gsl source - seems to be best I've found by far.
+    -- 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 a b
+    | a < 1 
+    = do
+        u <- realFloatStdUniform
+        x <- realFloatGamma (1 + a) b
+        return (x * u ** recip a)
+    | otherwise
+    = go
+        where
+            d = a - (1 / 3)
+            c = recip (3 * sqrt d) -- (1 / 3) / sqrt d
+            
+            go = do
+                x <- realFloatStdNormal
+                let cx = c * x
+                    v = (1 + cx) ^ 3
+                    
+                    x_2 = x * x
+                    x_4 = x_2 * x_2
+                
+                if cx <= (-1)
+                    then go
+                    else do
+                        u <- realFloatStdUniform
+                        
+                        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
+
+realFloatErlang :: (Integral a, RealFloat b) => a -> RVar b
+realFloatErlang a = realFloatGamma (fromIntegral a) 1
+
+gamma :: (Distribution Gamma a) => a -> a -> RVar a
+gamma a b = sample (Gamma a b)
+
+erlang :: (Distribution Gamma a, Integral b, Num a) => b -> a -> RVar a
+erlang a b = sample (Gamma (fromIntegral a) b)
+
+data Gamma a = Gamma a a
+
+instance RealFloat a => Distribution Gamma a where
+    rvar (Gamma a b) = realFloatGamma a b
diff --git a/src/Data/Random/Distribution/Normal.hs b/src/Data/Random/Distribution/Normal.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution/Normal.hs
@@ -0,0 +1,70 @@
+{-
+ -      ``Data/Random/Distribution/Normal''
+ -  
+ -  Quick and dirty implementation - eventually something faster probably 
+ -  ought to be implemented instead.
+ -}
+{-# LANGUAGE
+    MultiParamTypeClasses, FlexibleInstances, FlexibleContexts,
+    UndecidableInstances
+  #-}
+
+module Data.Random.Distribution.Normal where
+
+import Data.Random.Source
+import Data.Random.Distribution
+import Data.Random.Distribution.Uniform
+import Data.Random.RVar
+
+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)
+    let r = sqrt (-2 * log u)
+        
+        x = r * cos t
+        y = r * sin t
+    return (x,y)
+
+-- slightly slower
+knuthPolarNormalPair :: (Floating a, Ord a, Distribution Uniform a) => RVar (a,a)
+knuthPolarNormalPair = do
+    v1 <- uniform (-1) 1
+    v2 <- uniform (-1) 1
+    
+    let s = v1*v1 + v2*v2
+    if s >= 1
+        then knuthPolarNormalPair
+        else return $ if s == 0
+            then (0,0)
+            else let scale = sqrt (-2 * log s / s) 
+                  in (v1 * scale, v2 * scale)
+
+realFloatStdNormal :: RealFloat a => RVar a
+realFloatStdNormal = do
+    u <- realFloatStdUniform
+    t <- realFloatStdUniform
+    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
+    rvar (Normal m s) = do
+        x <- liftM fst normalPair
+        return (x * s + m)
+
+stdNormal :: Distribution Normal a => RVar a
+stdNormal = rvar StdNormal
+
+normal :: Distribution Normal a => a -> a -> RVar a
+normal m s = rvar (Normal m s)
diff --git a/src/Data/Random/Distribution/Poisson.hs b/src/Data/Random/Distribution/Poisson.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution/Poisson.hs
@@ -0,0 +1,67 @@
+{-
+ -      ``Data/Random/Distribution/Poisson''
+ -}
+{-# LANGUAGE
+    MultiParamTypeClasses,
+    FlexibleInstances, FlexibleContexts, UndecidableInstances
+  #-}
+
+module Data.Random.Distribution.Poisson where
+
+import Data.Random.Source
+import Data.Random.Distribution
+import Data.Random.RVar
+
+import Data.Random.Distribution.Uniform
+import Data.Random.Distribution.Gamma
+import Data.Random.Distribution.Binomial
+
+import Data.Int
+import Data.Word
+
+import Control.Monad
+
+-- from Knuth, with interpretation help from gsl sources
+integralPoisson :: (Integral a, RealFloat b) => b -> RVar a
+integralPoisson mu = psn 0 mu
+    where
+        psn k mu
+            | mu > 10   = do
+                let m = floor (mu * (7/8))
+            
+                x <- realFloatErlang m
+                if x >= mu
+                    then do
+                        b <- integralBinomial (m - 1) (mu / x)
+                        return (k + b)
+                    else psn (k + m) (mu - x)
+            
+            | otherwise = prod 1 k
+                where
+                    emu = exp (-mu)
+                
+                    prod p k = do
+                        u <- realFloatStdUniform
+                        if p * u > emu
+                            then prod (p * u) (k + 1)
+                            else return k
+
+
+poisson :: (Distribution (Poisson b) a) => b -> RVar a
+poisson mu = sample (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
+
+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)
diff --git a/src/Data/Random/Distribution/Triangular.hs b/src/Data/Random/Distribution/Triangular.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution/Triangular.hs
@@ -0,0 +1,36 @@
+{-
+ -      ``Data/Random/Distribution/Triangular''
+ -}
+{-# LANGUAGE
+    MultiParamTypeClasses,
+    FlexibleInstances
+  #-}
+
+module Data.Random.Distribution.Triangular where
+
+import Data.Random.RVar
+import Data.Random.Distribution
+import Data.Random.Distribution.Uniform
+
+data Triangular a = Triangular
+    { triLower  :: a
+    , triMid    :: a
+    , triUpper  :: a
+    } deriving (Eq, Show)
+
+realFloatTriangular :: (RealFloat a) => a -> a -> a -> RVar a
+realFloatTriangular a b c
+    | a <= b && b <= c
+    = do
+        let p = (c-b)/(c-a)
+        u <- realFloatStdUniform
+        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
+        return (b - ((1 - sqrt x) * (b-d)))
+
+instance RealFloat a => Distribution Triangular a where
+    rvar (Triangular a b c) = realFloatTriangular a b c
diff --git a/src/Data/Random/Distribution/Uniform.hs b/src/Data/Random/Distribution/Uniform.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution/Uniform.hs
@@ -0,0 +1,122 @@
+{-
+ -      ``Data/Random/Distribution/Uniform''
+ -}
+{-# LANGUAGE
+    MultiParamTypeClasses, FunctionalDependencies,
+    FlexibleContexts, FlexibleInstances, 
+    UndecidableInstances
+  #-}
+
+module Data.Random.Distribution.Uniform
+    ( Uniform(..)
+	, UniformByClassification(..)
+	, uniform
+	
+    , StdUniform(..)
+    , StdUniformByClassification(..)
+    , stdUniform
+    
+    , integralUniform
+    , realFloatUniform
+    
+    , boundedStdUniform
+    , boundedEnumStdUniform
+    , realFloatStdUniform
+    ) where
+
+import Data.Random.Internal.Classification
+
+import Data.Random.Source
+import Data.Random.Distribution
+import Data.Random.RVar
+
+import Data.Word
+import Data.Int
+import Data.Bits
+import Data.List
+
+import Control.Monad.Loops
+
+integralUniform a b
+    | a > b     = compute b a
+    | otherwise = compute a b
+    where
+        compute a b = do
+            let m = 1 + toInteger b - toInteger a
+            
+            let bytes = bytesNeeded m
+                maxXpossible = (powersOf256 !! bytes) - 1
+            
+            x <- iterateUntil (maxXpossible - maxXpossible `mod` m >) (nByteInteger bytes)
+            return (a + fromInteger (x `mod` m))
+
+
+bytesNeeded x = case findIndex (> x) powersOf256 of
+    Just x -> x
+powersOf256 = iterate (256 *) 1
+
+boundedStdUniform :: (Distribution Uniform a, Bounded a) => RVar a
+boundedStdUniform = uniform minBound maxBound
+
+boundedEnumStdUniform :: (Enum a, Bounded a) => RVar a
+boundedEnumStdUniform = enumUniform minBound maxBound
+
+-- (0,1]
+realFloatStdUniform :: RealFloat a => RVar a
+realFloatStdUniform | False     = return one
+                    | otherwise = do
+    let bitsNeeded  = floatDigits one
+        (_, e) = decodeFloat one
+    
+    x <- nBitInteger bitsNeeded
+    if x == 0
+        then return 1
+        else return (encodeFloat x (e-1))
+    
+    where one = 1
+
+realFloatUniform :: RealFloat a => a -> a -> RVar a
+realFloatUniform 0 1 = realFloatStdUniform
+realFloatUniform a b = do
+    x <- realFloatStdUniform
+    return (a + x * (b - a))
+
+enumUniform :: Enum a => a -> a -> RVar a
+enumUniform a b = do
+    x <- integralUniform (fromEnum a) (fromEnum b)
+    return (toEnum x)
+
+uniform :: Distribution Uniform a => a -> a -> RVar a
+uniform a b = rvar (Uniform a b)
+
+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
+
+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 StdUniformByClassification c t => Distribution StdUniform t
+    where rvar _ = stdUniformByClassification
+
+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 (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
diff --git a/src/Data/Random/Internal/Classification.hs b/src/Data/Random/Internal/Classification.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Internal/Classification.hs
@@ -0,0 +1,102 @@
+{-# 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
diff --git a/src/Data/Random/Internal/Words.hs b/src/Data/Random/Internal/Words.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Internal/Words.hs
@@ -0,0 +1,36 @@
+{-
+ -      ``Data/Random/Internal/Words''
+ -}
+
+-- |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.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]
+
+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
+
+
+bytesToWord :: [Word8] -> Word64
+bytesToWord bs = unsafePerformIO . allocaBytes 8 $ \p -> do
+    zipWithM (pokeElemOff p) [0..7] (bs ++ repeat 0)
+    peek (castPtr p)
diff --git a/src/Data/Random/RVar.hs b/src/Data/Random/RVar.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/RVar.hs
@@ -0,0 +1,100 @@
+{-
+ -      ``Data/Random/RVar''
+ -}
+{-# LANGUAGE
+    RankNTypes,
+    MultiParamTypeClasses,
+    FlexibleInstances
+  #-}
+
+-- |Random variables.  An 'RVar' is a sampleable random variable.  Because
+-- probability distributions form a monad, they are quite easy to work with
+-- in the standard Haskell monadic styles.  For examples, see the source for
+-- any of the 'Distribution' instances - they all are defined in terms of
+-- 'RVar's.
+module Data.Random.RVar
+    ( RVar
+    
+    , nByteInteger
+    , nBitInteger
+    ) where
+
+import Data.Random.Distribution
+import Data.Random.Source
+
+import Data.Word
+import Data.Bits
+
+import Control.Applicative
+import Control.Monad
+
+-- |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 }
+
+instance Functor RVar 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 Applicative RVar where
+    pure  = return
+    (<*>) = ap
+
+instance Distribution RVar a where
+    rvar = id
+    sampleFrom src x = runDistM x src
+
+instance MonadRandom RVar where
+    getRandomBytes n = RVar (\s -> getRandomBytesFrom s n)
+    getRandomWords n = RVar (\s -> getRandomWordsFrom s n)
+
+-- some 'fundamental' RVars
+-- this maybe ought to even be a part of the RandomSource class...
+-- |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
+
+-- |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)
diff --git a/src/Data/Random/RVar.hs-boot b/src/Data/Random/RVar.hs-boot
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/RVar.hs-boot
@@ -0,0 +1,16 @@
+{-
+ -      ``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
diff --git a/src/Data/Random/Source.hs b/src/Data/Random/Source.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Source.hs
@@ -0,0 +1,75 @@
+{-
+ -      ``Data/Random/Source''
+ -}
+{-# LANGUAGE
+    MultiParamTypeClasses, FlexibleInstances
+  #-}
+
+module Data.Random.Source
+    ( MonadRandom(..)
+    , RandomSource(..)
+    ) where
+
+import Data.Word
+import Data.Bits
+import Data.List
+
+import Data.Random.Internal.Words
+
+-- |A typeclass for monads with a chosen source of entropy.  For example,
+-- 'RVar' is such a monad - the source from which it is (eventually) sampled
+-- is the only source from which a random variable is permitted to draw, so
+-- when directly requesting entropy for a random variable these functions
+-- are used.
+-- 
+-- The minimal definition is either 'getRandomBytes' or 'getRandomWords'.
+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
+        
+    -- |alternate basis function, providing access to larger chunks
+    getRandomWords :: Int -> m [Word64]
+    getRandomWords n = do
+        bs <- getRandomBytes (n `shiftL` 3)
+        return (bytesToWords bs)
+
+-- |A source of entropy which can be used in the given monad.
+--
+-- The minimal definition is either 'getRandomBytesFrom' or 'getRandomWordsFrom'
+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
+        
+    
+    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
+
+instance Monad m => RandomSource m (Int -> m [Word64]) where
+    getRandomWordsFrom = id
+
diff --git a/src/Data/Random/Source/DevRandom.hs b/src/Data/Random/Source/DevRandom.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Source/DevRandom.hs
@@ -0,0 +1,24 @@
+{-
+ -      ``Data/Random/Source/DevRandom''
+ -}
+{-# LANGUAGE
+    MultiParamTypeClasses
+  #-}
+
+module Data.Random.Source.DevRandom where
+
+import Data.Random.Source
+
+import GHC.IOBase (unsafePerformIO)
+import Data.ByteString (hGet, unpack)
+import System.IO (openBinaryFile, IOMode(..))
+
+-- |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)
+
+instance RandomSource IO DevRandom where
+    getRandomBytesFrom DevRandom n = fmap unpack (hGet devRandom n)
+
diff --git a/src/Data/Random/Source/PureMT.hs b/src/Data/Random/Source/PureMT.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Source/PureMT.hs
@@ -0,0 +1,55 @@
+{-
+ -      ``Data/Random/Source/PureMT''
+ -}
+{-# LANGUAGE
+    MultiParamTypeClasses,
+    FlexibleContexts, FlexibleInstances
+  #-}
+
+module Data.Random.Source.PureMT where
+
+import Data.Random.Source
+import System.Random.Mersenne.Pure64
+
+import Data.StateRef
+import Data.Word
+
+import Control.Monad.State
+
+-- |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
+-- 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 [])
+    
+    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\"
+-- 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
+    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
+    put newMt
+    return ws
+
+instance MonadRandom (State PureMT) where
+    getRandomWords = getRandomWordsFromMTState
+
+instance Monad m => MonadRandom (StateT PureMT m) where
+    getRandomWords = getRandomWordsFromMTState
diff --git a/src/Data/Random/Source/Std.hs b/src/Data/Random/Source/Std.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Source/Std.hs
@@ -0,0 +1,20 @@
+{-
+ -      ``Data/Random/Source/Std''
+ -}
+{-# LANGUAGE
+    MultiParamTypeClasses, FlexibleInstances
+  #-}
+
+module Data.Random.Source.Std where
+
+import Data.Random.Source
+
+-- |A token representing the \"standard\" entropy source in a 'MonadRandom'
+-- monad.  Its sole purpose is to make the following true (when the types check):
+--
+-- > sampleFrom StdRandom === sample
+data StdRandom = StdRandom
+
+instance MonadRandom m => RandomSource m StdRandom where
+    getRandomBytesFrom StdRandom = getRandomBytes
+    getRandomWordsFrom StdRandom = getRandomWords
diff --git a/src/Data/Random/Source/StdGen.hs b/src/Data/Random/Source/StdGen.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Source/StdGen.hs
@@ -0,0 +1,88 @@
+{-
+ -      ``Data/Random/Source/StdGen''
+ -}
+{-# LANGUAGE
+    MultiParamTypeClasses, FlexibleInstances, UndecidableInstances
+  #-}
+
+module Data.Random.Source.StdGen where
+
+import Data.Random.Source
+import System.Random
+import Control.Monad
+import Control.Monad.State
+import Data.StateRef
+import Data.Word
+
+instance (ModifyRef (IORef   StdGen) m StdGen) => RandomSource m (IORef   StdGen) where
+    getRandomBytesFrom = getRandomBytesFromRandomGenRef
+    getRandomWordsFrom = getRandomWordsFromRandomGenRef
+instance (ModifyRef (TVar    StdGen) m StdGen) => RandomSource m (TVar    StdGen) where
+    getRandomBytesFrom = getRandomBytesFromRandomGenRef
+    getRandomWordsFrom = getRandomWordsFromRandomGenRef
+instance (ModifyRef (STRef s StdGen) m StdGen) => RandomSource m (STRef s StdGen) where
+    getRandomBytesFrom = getRandomBytesFromRandomGenRef
+    getRandomWordsFrom = getRandomWordsFromRandomGenRef
+
+getRandomBytesFromStdGenIO :: Int -> IO [Word8]
+getRandomBytesFromStdGenIO n = do
+    ints <- replicateM n (randomRIO (0, 255))
+    let bytes = map fromIntegral (ints :: [Int])
+    return bytes
+
+-- |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
+-- 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\"
+-- 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
+    g <- get
+    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
+    g <- get
+    case randomR (0,2^64-1) g of
+        (i,g) -> do
+            put g
+            return (fromInteger i)
+
+instance MonadRandom (State StdGen) where
+    getRandomBytes = getRandomBytesFromRandomGenState
+    getRandomWords = getRandomWordsFromRandomGenState
+
+instance Monad m => MonadRandom (StateT StdGen m) where
+    getRandomBytes = getRandomBytesFromRandomGenState
+    getRandomWords = getRandomWordsFromRandomGenState