diff --git a/random-fu.cabal b/random-fu.cabal
--- a/random-fu.cabal
+++ b/random-fu.cabal
@@ -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,
diff --git a/src/Data/Random.hs b/src/Data/Random.hs
--- a/src/Data/Random.hs
+++ b/src/Data/Random.hs
@@ -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
 
diff --git a/src/Data/Random/Distribution.hs b/src/Data/Random/Distribution.hs
--- a/src/Data/Random/Distribution.hs
+++ b/src/Data/Random/Distribution.hs
@@ -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)
+
diff --git a/src/Data/Random/Distribution.hs-boot b/src/Data/Random/Distribution.hs-boot
deleted file mode 100644
--- a/src/Data/Random/Distribution.hs-boot
+++ /dev/null
@@ -1,10 +0,0 @@
-{-
- -      ``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
--- a/src/Data/Random/Distribution/Bernoulli.hs
+++ b/src/Data/Random/Distribution/Bernoulli.hs
@@ -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
+
+
diff --git a/src/Data/Random/Distribution/Beta.hs b/src/Data/Random/Distribution/Beta.hs
--- a/src/Data/Random/Distribution/Beta.hs
+++ b/src/Data/Random/Distribution/Beta.hs
@@ -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
diff --git a/src/Data/Random/Distribution/Binomial.hs b/src/Data/Random/Distribution/Binomial.hs
--- a/src/Data/Random/Distribution/Binomial.hs
+++ b/src/Data/Random/Distribution/Binomial.hs
@@ -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
diff --git a/src/Data/Random/Distribution/Discrete.hs b/src/Data/Random/Distribution/Discrete.hs
--- a/src/Data/Random/Distribution/Discrete.hs
+++ b/src/Data/Random/Distribution/Discrete.hs
@@ -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
diff --git a/src/Data/Random/Distribution/Exponential.hs b/src/Data/Random/Distribution/Exponential.hs
--- a/src/Data/Random/Distribution/Exponential.hs
+++ b/src/Data/Random/Distribution/Exponential.hs
@@ -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
diff --git a/src/Data/Random/Distribution/Gamma.hs b/src/Data/Random/Distribution/Gamma.hs
--- a/src/Data/Random/Distribution/Gamma.hs
+++ b/src/Data/Random/Distribution/Gamma.hs
@@ -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
diff --git a/src/Data/Random/Distribution/Normal.hs b/src/Data/Random/Distribution/Normal.hs
--- a/src/Data/Random/Distribution/Normal.hs
+++ b/src/Data/Random/Distribution/Normal.hs
@@ -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
diff --git a/src/Data/Random/Distribution/Poisson.hs b/src/Data/Random/Distribution/Poisson.hs
--- a/src/Data/Random/Distribution/Poisson.hs
+++ b/src/Data/Random/Distribution/Poisson.hs
@@ -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
diff --git a/src/Data/Random/Distribution/Triangular.hs b/src/Data/Random/Distribution/Triangular.hs
--- a/src/Data/Random/Distribution/Triangular.hs
+++ b/src/Data/Random/Distribution/Triangular.hs
@@ -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
diff --git a/src/Data/Random/Distribution/Uniform.hs b/src/Data/Random/Distribution/Uniform.hs
--- a/src/Data/Random/Distribution/Uniform.hs
+++ b/src/Data/Random/Distribution/Uniform.hs
@@ -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
+
diff --git a/src/Data/Random/Internal/Classification.hs b/src/Data/Random/Internal/Classification.hs
deleted file mode 100644
--- a/src/Data/Random/Internal/Classification.hs
+++ /dev/null
@@ -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
diff --git a/src/Data/Random/Internal/Words.hs b/src/Data/Random/Internal/Words.hs
--- a/src/Data/Random/Internal/Words.hs
+++ b/src/Data/Random/Internal/Words.hs
@@ -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)
diff --git a/src/Data/Random/Lift.hs b/src/Data/Random/Lift.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Lift.hs
@@ -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
diff --git a/src/Data/Random/List.hs b/src/Data/Random/List.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/List.hs
@@ -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)
diff --git a/src/Data/Random/RVar.hs b/src/Data/Random/RVar.hs
--- a/src/Data/Random/RVar.hs
+++ b/src/Data/Random/RVar.hs
@@ -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))
diff --git a/src/Data/Random/RVar.hs-boot b/src/Data/Random/RVar.hs-boot
deleted file mode 100644
--- a/src/Data/Random/RVar.hs-boot
+++ /dev/null
@@ -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
diff --git a/src/Data/Random/Sample.hs b/src/Data/Random/Sample.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Sample.hs
@@ -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
diff --git a/src/Data/Random/Source.hs b/src/Data/Random/Source.hs
--- a/src/Data/Random/Source.hs
+++ b/src/Data/Random/Source.hs
@@ -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
diff --git a/src/Data/Random/Source/DevRandom.hs b/src/Data/Random/Source/DevRandom.hs
--- a/src/Data/Random/Source/DevRandom.hs
+++ b/src/Data/Random/Source/DevRandom.hs
@@ -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)
diff --git a/src/Data/Random/Source/PureMT.hs b/src/Data/Random/Source/PureMT.hs
--- a/src/Data/Random/Source/PureMT.hs
+++ b/src/Data/Random/Source/PureMT.hs
@@ -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
+
diff --git a/src/Data/Random/Source/Std.hs b/src/Data/Random/Source/Std.hs
--- a/src/Data/Random/Source/Std.hs
+++ b/src/Data/Random/Source/Std.hs
@@ -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
diff --git a/src/Data/Random/Source/StdGen.hs b/src/Data/Random/Source/StdGen.hs
--- a/src/Data/Random/Source/StdGen.hs
+++ b/src/Data/Random/Source/StdGen.hs
@@ -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