diff --git a/random-fu/changelog.md b/random-fu/changelog.md
index a22da78..6776e24 100644
--- a/random-fu/changelog.md
+++ b/random-fu/changelog.md
@@ -1,3 +1,7 @@
+* Chnages in 0.3.0.0:
+
+  * Drop usage of `random-source` in favor of `random`
+
 * Changes in 0.2.7.7: Update to random-1.2. Revert 0.2.7.6 changes (which added an extra constraint to `Data.Random.Sample.sampleState` and `Data.Random.Sample.sampleStateT`).
 
 * Changes in 0.2.7.4: Compatibility with ghc 8.8.
diff --git a/random-fu/random-fu.cabal b/random-fu/random-fu.cabal
index e7e3db9..9abeace 100644
--- a/random-fu/random-fu.cabal
+++ b/random-fu/random-fu.cabal
@@ -1,5 +1,5 @@
 name:                   random-fu
-version:                0.2.7.7
+version:                1.0.0.0
 stability:              provisional
 
 cabal-version:          >= 1.10
@@ -12,21 +12,21 @@ homepage:               https://github.com/mokus0/random-fu
 
 category:               Math
 synopsis:               Random number generation
-description:            Random number generation based on modeling random 
+description:            Random number generation based on modeling random
                         variables in two complementary ways: first, by the
                         parameters of standard mathematical distributions and,
                         second, by an abstract type ('RVar') which can be
                         composed and manipulated monadically and sampled in
                         either monadic or \"pure\" styles.
                         .
-                        The primary purpose of this library is to support 
+                        The primary purpose of this library is to support
                         defining and sampling a wide variety of high quality
                         random variables.  Quality is prioritized over speed,
                         but performance is an important goal too.
                         .
-                        In my testing, I have found it capable of speed 
+                        In my testing, I have found it capable of speed
                         comparable to other Haskell libraries, but still
-                        a fair bit slower than straight C implementations of 
+                        a fair bit slower than straight C implementations of
                         the same algorithms.
 
 tested-with:            GHC == 7.10.3
@@ -83,25 +83,24 @@ Library
   else
     cpp-options:        -Dold_Fixed
     build-depends:      base >= 4 && <4.2
-  
+
   if flag(mtl2)
     build-depends:      mtl == 2.*
     cpp-options:        -DMTL2
   else
     build-depends:      mtl == 1.*
-  
+
   build-depends:        math-functions,
                         monad-loops >= 0.3.0.1,
                         random >= 1.2 && < 1.3,
                         random-shuffle,
-                        random-source == 0.3.*,
-                        rvar == 0.2.*,
+                        rvar >= 0.3,
                         syb,
                         template-haskell,
                         transformers,
                         vector >= 0.7,
                         erf
-  
+
   if impl(ghc == 7.2.1)
     -- Doesn't work under GHC 7.2.1 due to
     -- http://hackage.haskell.org/trac/ghc/ticket/5410
diff --git a/random-fu/src/Data/Random.hs b/random-fu/src/Data/Random.hs
index 58bcc0d..1218132 100644
--- a/random-fu/src/Data/Random.hs
+++ b/random-fu/src/Data/Random.hs
@@ -1,39 +1,39 @@
 -- |Flexible modeling and sampling of random variables.
 --
--- The central abstraction in this library is the concept of a random 
--- variable.  It is not fully formalized in the standard measure-theoretic 
--- language, but rather is informally defined as a \"thing you can get random 
--- values out of\".  Different random variables may have different types of 
+-- The central abstraction in this library is the concept of a random
+-- variable.  It is not fully formalized in the standard measure-theoretic
+-- language, but rather is informally defined as a \"thing you can get random
+-- values out of\".  Different random variables may have different types of
 -- values they can return or the same types but different probabilities for
 -- each value they can return.  The random values you get out of them are
 -- traditionally called \"random variates\".
--- 
--- Most imperative-language random number libraries are all about obtaining 
--- and manipulating random variates.  This one is about defining, manipulating 
--- and sampling random variables.  Computationally, the distinction is small 
--- and mostly just a matter of perspective, but from a program design 
+--
+-- Most imperative-language random number libraries are all about obtaining
+-- and manipulating random variates.  This one is about defining, manipulating
+-- and sampling random variables.  Computationally, the distinction is small
+-- and mostly just a matter of perspective, but from a program design
 -- perspective it provides both a powerfully composable abstraction and a
 -- very useful separation of concerns.
--- 
+--
 -- Abstract random variables as implemented by 'RVar' are composable.  They can
 -- be defined in a monadic / \"imperative\" style that amounts to manipulating
 -- variates, but with strict type-level isolation.  Concrete random variables
 -- are also provided, but they do not compose as generically.  The 'Distribution'
--- type class allows concrete random variables to \"forget\" their concreteness 
--- so that they can be composed.  For examples of both, see the documentation 
--- for 'RVar' and 'Distribution', as well as the code for any of the concrete 
+-- type class allows concrete random variables to \"forget\" their concreteness
+-- so that they can be composed.  For examples of both, see the documentation
+-- for 'RVar' and 'Distribution', as well as the code for any of the concrete
 -- distributions such as 'Uniform', 'Gamma', etc.
--- 
+--
 -- Both abstract and concrete random variables can be sampled (despite the
 -- types GHCi may list for the functions) by the functions in "Data.Random.Sample".
--- 
+--
 -- Random variable sampling is done with regard to a generic basis of primitive
--- random variables defined in "Data.Random.Internal.Primitives".  This basis 
+-- random variables defined in "Data.Random.Internal.Primitives".  This basis
 -- is very low-level and the actual set of primitives is still fairly experimental,
 -- which is why it is in the \"Internal\" sub-heirarchy.  User-defined variables
 -- should use the existing high-level variables such as 'Uniform' and 'Normal'
 -- rather than these basis variables.  "Data.Random.Source" defines classes for
--- entropy sources that provide implementations of these primitive variables. 
+-- entropy sources that provide implementations of these primitive variables.
 -- Several implementations are available in the Data.Random.Source.* modules.
 module Data.Random
     ( -- * Random variables
@@ -43,32 +43,26 @@ module Data.Random
 
       -- ** Concrete ('Distribution')
       Distribution(..), CDF(..), PDF(..),
-      
+
       -- * Sampling random variables
-      Sampleable(..), sample, sampleState, sampleStateT,
-      
+      Sampleable(..), sample, sampleState, samplePure,
+
       -- * A few very common distributions
       Uniform(..), uniform, uniformT,
       StdUniform(..), stdUniform, stdUniformT,
       Normal(..), normal, stdNormal, normalT, stdNormalT,
       Gamma(..), gamma, gammaT,
-      
+
       -- * Entropy Sources
-      MonadRandom, RandomSource, StdRandom(..),
-      
+      StatefulGen, RandomGen,
+
       -- * Useful list-based operations
       randomElement,
       shuffle, shuffleN, shuffleNofM
-      
+
     ) where
 
 import Data.Random.Sample
-import Data.Random.Source (MonadRandom, RandomSource)
-import Data.Random.Source.IO ()
-import Data.Random.Source.MWC ()
-import Data.Random.Source.StdGen ()
-import Data.Random.Source.PureMT ()
-import Data.Random.Source.Std
 import Data.Random.Distribution
 import Data.Random.Distribution.Gamma
 import Data.Random.Distribution.Normal
@@ -78,3 +72,4 @@ import Data.Random.Lift ()
 import Data.Random.List
 import Data.Random.RVar
 
+import System.Random.Stateful (StatefulGen, RandomGen)
diff --git a/random-fu/src/Data/Random/Distribution.hs b/random-fu/src/Data/Random/Distribution.hs
index 0c460f3..acc8093 100644
--- a/random-fu/src/Data/Random/Distribution.hs
+++ b/random-fu/src/Data/Random/Distribution.hs
@@ -57,7 +57,7 @@ class Distribution d t where
     -- |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 :: d t -> RVarT n t
+    rvarT :: Functor n => d t -> RVarT n t
     rvarT d = lift (rvar d)
 
 -- FIXME: I am not sure about giving default instances
diff --git a/random-fu/src/Data/Random/Distribution/Bernoulli.hs b/random-fu/src/Data/Random/Distribution/Bernoulli.hs
index 1027fd8..ccedcee 100644
--- a/random-fu/src/Data/Random/Distribution/Bernoulli.hs
+++ b/random-fu/src/Data/Random/Distribution/Bernoulli.hs
@@ -27,12 +27,12 @@ bernoulli p = rvar (Bernoulli p)
 -- |Generate a Bernoulli process with the given probability.  For @Bool@ results,
 -- @bernoulli p@ will return True (p*100)% of the time and False otherwise.
 -- For numerical types, True is replaced by 1 and False by 0.
-bernoulliT :: Distribution (Bernoulli b) a => b -> RVarT m a
+bernoulliT :: (Distribution (Bernoulli b) a, Functor m) => b -> RVarT m a
 bernoulliT p = rvarT (Bernoulli p)
 
 -- |A random variable whose value is 'True' the given fraction of the time
 -- and 'False' the rest.
-boolBernoulli :: (Fractional a, Ord a, Distribution StdUniform a) => a -> RVarT m Bool
+boolBernoulli :: (Fractional a, Ord a, Distribution StdUniform a, Functor m) => a -> RVarT m Bool
 boolBernoulli p = do
     x <- stdUniformT
     return (x <= p)
@@ -43,7 +43,7 @@ boolBernoulliCDF p False = (1 - realToFrac 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 -> RVarT m a
+generalBernoulli :: (Distribution (Bernoulli b) Bool, Functor m) => a -> a -> b -> RVarT m a
 generalBernoulli f t p = do
     x <- bernoulliT p
     return (if x then t else f)
diff --git a/random-fu/src/Data/Random/Distribution/Beta.hs b/random-fu/src/Data/Random/Distribution/Beta.hs
index 754f388..56422a9 100644
--- a/random-fu/src/Data/Random/Distribution/Beta.hs
+++ b/random-fu/src/Data/Random/Distribution/Beta.hs
@@ -18,9 +18,9 @@ import Data.Random.Distribution.Uniform
 
 import Numeric.SpecFunctions
 
-{-# SPECIALIZE fractionalBeta :: Float  -> Float  -> RVarT m Float #-}
-{-# SPECIALIZE fractionalBeta :: Double -> Double -> RVarT m Double #-}
-fractionalBeta :: (Fractional a, Eq a, Distribution Gamma a, Distribution StdUniform a) => a -> a -> RVarT m a
+{-# SPECIALIZE fractionalBeta :: Functor m => Float  -> Float  -> RVarT m Float #-}
+{-# SPECIALIZE fractionalBeta :: Functor m => Double -> Double -> RVarT m Double #-}
+fractionalBeta :: (Fractional a, Eq a, Distribution Gamma a, Distribution StdUniform a, Functor m) => a -> a -> RVarT m a
 fractionalBeta 1 1 = stdUniformT
 fractionalBeta a b = do
     x <- gammaT a 1
@@ -32,9 +32,9 @@ fractionalBeta a b = do
 beta :: Distribution Beta a => a -> a -> RVar a
 beta a b = rvar (Beta a b)
 
-{-# SPECIALIZE betaT :: Float  -> Float  -> RVarT m Float #-}
-{-# SPECIALIZE betaT :: Double -> Double -> RVarT m Double #-}
-betaT :: Distribution Beta a => a -> a -> RVarT m a
+{-# SPECIALIZE betaT :: Functor m => Float  -> Float  -> RVarT m Float #-}
+{-# SPECIALIZE betaT :: Functor m => Double -> Double -> RVarT m Double #-}
+betaT :: (Distribution Beta a, Functor m) => a -> a -> RVarT m a
 betaT a b = rvarT (Beta a b)
 
 data Beta a = Beta a a
diff --git a/random-fu/src/Data/Random/Distribution/Binomial.hs b/random-fu/src/Data/Random/Distribution/Binomial.hs
index 802a6c8..efb0e15 100644
--- a/random-fu/src/Data/Random/Distribution/Binomial.hs
+++ b/random-fu/src/Data/Random/Distribution/Binomial.hs
@@ -25,10 +25,10 @@ import Numeric ( log1p )
     -- 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 -> RVarT m a
+integralBinomial :: (Integral a, Floating b, Ord b, Distribution Beta b, Distribution StdUniform b, Functor m) => a -> b -> RVarT m a
 integralBinomial = bin 0
     where
-        bin :: (Integral a, Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => a -> a -> b -> RVarT m a
+        bin :: (Integral a, Floating b, Ord b, Distribution Beta b, Distribution StdUniform b, Functor m) => a -> a -> b -> RVarT m a
         bin !k !t !p
             | t > 10    = do
                 let a = 1 + t `div` 2
@@ -118,15 +118,15 @@ floatingBinomialLogPDF t p x = logPdf (Binomial (truncate t :: Integer) p) (floo
 binomial :: Distribution (Binomial b) a => a -> b -> RVar a
 binomial t p = rvar (Binomial t p)
 
-{-# SPECIALIZE binomialT :: Int     -> Float  -> RVarT m Int #-}
-{-# SPECIALIZE binomialT :: Int     -> Double -> RVarT m Int #-}
-{-# SPECIALIZE binomialT :: Integer -> Float  -> RVarT m Integer #-}
-{-# SPECIALIZE binomialT :: Integer -> Double -> RVarT m Integer #-}
-{-# SPECIALIZE binomialT :: Float   -> Float  -> RVarT m Float  #-}
-{-# SPECIALIZE binomialT :: Float   -> Double -> RVarT m Float  #-}
-{-# SPECIALIZE binomialT :: Double  -> Float  -> RVarT m Double #-}
-{-# SPECIALIZE binomialT :: Double  -> Double -> RVarT m Double #-}
-binomialT :: Distribution (Binomial b) a => a -> b -> RVarT m a
+{-# SPECIALIZE binomialT :: Functor m => Int     -> Float  -> RVarT m Int #-}
+{-# SPECIALIZE binomialT :: Functor m => Int     -> Double -> RVarT m Int #-}
+{-# SPECIALIZE binomialT :: Functor m => Integer -> Float  -> RVarT m Integer #-}
+{-# SPECIALIZE binomialT :: Functor m => Integer -> Double -> RVarT m Integer #-}
+{-# SPECIALIZE binomialT :: Functor m => Float   -> Float  -> RVarT m Float  #-}
+{-# SPECIALIZE binomialT :: Functor m => Float   -> Double -> RVarT m Float  #-}
+{-# SPECIALIZE binomialT :: Functor m => Double  -> Float  -> RVarT m Double #-}
+{-# SPECIALIZE binomialT :: Functor m => Double  -> Double -> RVarT m Double #-}
+binomialT :: (Distribution (Binomial b) a, Functor m) => a -> b -> RVarT m a
 binomialT t p = rvarT (Binomial t p)
 
 data Binomial b a = Binomial a b
diff --git a/random-fu/src/Data/Random/Distribution/Categorical.hs b/random-fu/src/Data/Random/Distribution/Categorical.hs
index 575a0a2..c56cf21 100644
--- a/random-fu/src/Data/Random/Distribution/Categorical.hs
+++ b/random-fu/src/Data/Random/Distribution/Categorical.hs
@@ -23,9 +23,7 @@ import Data.Random.Distribution.Uniform
 import Control.Arrow
 import Control.Monad
 import Control.Monad.ST
-import Data.Foldable (Foldable(foldMap))
 import Data.STRef
-import Data.Traversable (Traversable(traverse, sequenceA))
 
 import Data.List
 import Data.Function
@@ -39,7 +37,7 @@ categorical = rvar . fromList
 
 -- |Construct a 'Categorical' random process from a list of probabilities
 -- and categories, where the probabilities all sum to 1.
-categoricalT :: (Num p, Distribution (Categorical p) a) => [(p,a)] -> RVarT m a
+categoricalT :: (Num p, Distribution (Categorical p) a, Functor m) => [(p,a)] -> RVarT m a
 categoricalT = rvarT . fromList
 
 -- |Construct a 'Categorical' random variable from a list of weights
@@ -49,7 +47,7 @@ weightedCategorical = rvar . fromWeightedList
 
 -- |Construct a 'Categorical' random process from a list of weights
 -- and categories. The weights do /not/ have to sum to 1.
-weightedCategoricalT :: (Fractional p, Eq p, Distribution (Categorical p) a) => [(p,a)] -> RVarT m a
+weightedCategoricalT :: (Fractional p, Eq p, Distribution (Categorical p) a, Functor m) => [(p,a)] -> RVarT m a
 weightedCategoricalT = rvarT . fromWeightedList
 
 -- | Construct a 'Categorical' distribution from a list of weighted categories.
diff --git a/random-fu/src/Data/Random/Distribution/ChiSquare.hs b/random-fu/src/Data/Random/Distribution/ChiSquare.hs
index 63eb9c4..9126b3c 100644
--- a/random-fu/src/Data/Random/Distribution/ChiSquare.hs
+++ b/random-fu/src/Data/Random/Distribution/ChiSquare.hs
@@ -16,7 +16,7 @@ import Numeric.SpecFunctions
 chiSquare :: Distribution ChiSquare t => Integer -> RVar t
 chiSquare = rvar . ChiSquare
 
-chiSquareT :: Distribution ChiSquare t => Integer -> RVarT m t
+chiSquareT :: (Distribution ChiSquare t, Functor m) => Integer -> RVarT m t
 chiSquareT = rvarT . ChiSquare
 
 newtype ChiSquare b = ChiSquare Integer
diff --git a/random-fu/src/Data/Random/Distribution/Dirichlet.hs b/random-fu/src/Data/Random/Distribution/Dirichlet.hs
index 7dbc4cc..68269c9 100644
--- a/random-fu/src/Data/Random/Distribution/Dirichlet.hs
+++ b/random-fu/src/Data/Random/Distribution/Dirichlet.hs
@@ -12,7 +12,7 @@ import Data.Random.Distribution.Gamma
 
 import Data.List
 
-fractionalDirichlet :: (Fractional a, Distribution Gamma a) => [a] -> RVarT m [a]
+fractionalDirichlet :: (Fractional a, Distribution Gamma a, Functor m) => [a] -> RVarT m [a]
 fractionalDirichlet []  = return []
 fractionalDirichlet [_] = return [1]
 fractionalDirichlet as = do
@@ -24,7 +24,7 @@ fractionalDirichlet as = do
 dirichlet :: Distribution Dirichlet [a] => [a] -> RVar [a]
 dirichlet as = rvar (Dirichlet as)
 
-dirichletT :: Distribution Dirichlet [a] => [a] -> RVarT m [a]
+dirichletT :: (Distribution Dirichlet [a], Functor m) => [a] -> RVarT m [a]
 dirichletT as = rvarT (Dirichlet as)
 
 newtype Dirichlet a = Dirichlet a deriving (Eq, Show)
diff --git a/random-fu/src/Data/Random/Distribution/Exponential.hs b/random-fu/src/Data/Random/Distribution/Exponential.hs
index df5f5b0..d2a3d8f 100644
--- a/random-fu/src/Data/Random/Distribution/Exponential.hs
+++ b/random-fu/src/Data/Random/Distribution/Exponential.hs
@@ -23,7 +23,7 @@ See also 'exponential'.
 -}
 newtype Exponential a = Exp a
 
-floatingExponential :: (Floating a, Distribution StdUniform a) => a -> RVarT m a
+floatingExponential :: (Floating a, Distribution StdUniform a, Functor m) => a -> RVarT m a
 floatingExponential lambdaRecip = do
     x <- stdUniformT
     return (negate (log x) * lambdaRecip)
@@ -48,7 +48,7 @@ A random variable transformer which samples from the exponential distribution.
 alternatively be viewed as an exponential random variable with parameter @lambda
 = 1 / mu@.
 -}
-exponentialT :: Distribution Exponential a => a -> RVarT m a
+exponentialT :: (Distribution Exponential a, Functor m) => a -> RVarT m a
 exponentialT = rvarT . Exp
 
 instance (Floating a, Distribution StdUniform a) => Distribution Exponential a where
diff --git a/random-fu/src/Data/Random/Distribution/Gamma.hs b/random-fu/src/Data/Random/Distribution/Gamma.hs
index bd27050..de9aa99 100644
--- a/random-fu/src/Data/Random/Distribution/Gamma.hs
+++ b/random-fu/src/Data/Random/Distribution/Gamma.hs
@@ -27,12 +27,12 @@ import Numeric.SpecFunctions
 
 -- |derived from  Marsaglia & Tang, "A Simple Method for generating gamma
 -- variables", ACM Transactions on Mathematical Software, Vol 26, No 3 (2000), p363-372.
-{-# SPECIALIZE mtGamma :: Double -> Double -> RVarT m Double #-}
-{-# SPECIALIZE mtGamma :: Float  -> Float  -> RVarT m Float  #-}
+{-# SPECIALIZE mtGamma :: Functor m => Double -> Double -> RVarT m Double #-}
+{-# SPECIALIZE mtGamma :: Functor m => Float  -> Float  -> RVarT m Float  #-}
 mtGamma
     :: (Floating a, Ord a,
         Distribution StdUniform a,
-        Distribution Normal a)
+        Distribution Normal a, Functor m)
     => a -> a -> RVarT m a
 mtGamma a b
     | a < 1     = do
@@ -64,13 +64,13 @@ mtGamma a b
 gamma :: (Distribution Gamma a) => a -> a -> RVar a
 gamma a b = rvar (Gamma a b)
 
-gammaT :: (Distribution Gamma a) => a -> a -> RVarT m a
+gammaT :: (Distribution Gamma a, Functor m) => a -> a -> RVarT m a
 gammaT a b = rvarT (Gamma a b)
 
 erlang :: (Distribution (Erlang a) b) => a -> RVar b
 erlang a = rvar (Erlang a)
 
-erlangT :: (Distribution (Erlang a) b) => a -> RVarT m b
+erlangT :: (Distribution (Erlang a) b, Functor m) => a -> RVarT m b
 erlangT a = rvarT (Erlang a)
 
 data    Gamma a    = Gamma a a
diff --git a/random-fu/src/Data/Random/Distribution/Multinomial.hs b/random-fu/src/Data/Random/Distribution/Multinomial.hs
index af4c9fa..11e6f9f 100644
--- a/random-fu/src/Data/Random/Distribution/Multinomial.hs
+++ b/random-fu/src/Data/Random/Distribution/Multinomial.hs
@@ -8,7 +8,7 @@ import Data.Random.Distribution.Binomial
 multinomial :: Distribution (Multinomial p) [a] => [p] -> a -> RVar [a]
 multinomial ps n = rvar (Multinomial ps n)
 
-multinomialT :: Distribution (Multinomial p) [a] => [p] -> a -> RVarT m [a]
+multinomialT :: (Distribution (Multinomial p) [a], Functor m) => [p] -> a -> RVarT m [a]
 multinomialT ps n = rvarT (Multinomial ps n)
 
 data Multinomial p a where
diff --git a/random-fu/src/Data/Random/Distribution/Normal.hs b/random-fu/src/Data/Random/Distribution/Normal.hs
index 90b9e11..efe77e6 100644
--- a/random-fu/src/Data/Random/Distribution/Normal.hs
+++ b/random-fu/src/Data/Random/Distribution/Normal.hs
@@ -22,14 +22,13 @@ module Data.Random.Distribution.Normal
     , knuthPolarNormalPair
     ) where
 
-import Data.Random.Internal.Words
 import Data.Bits
 
-import Data.Random.Source
 import Data.Random.Distribution
 import Data.Random.Distribution.Uniform
 import Data.Random.Distribution.Ziggurat
 import Data.Random.RVar
+import Data.Word
 
 import Data.Vector.Generic (Vector)
 import qualified Data.Vector as V
@@ -37,6 +36,8 @@ import qualified Data.Vector.Unboxed as UV
 
 import Data.Number.Erf
 
+import qualified System.Random.Stateful as Random
+
 -- |A random variable that produces a pair of independent
 -- normally-distributed values.
 normalPair :: (Floating a, Distribution StdUniform a) => RVar (a,a)
@@ -80,7 +81,7 @@ knuthPolarNormalPair = do
 
 -- |Draw from the tail of a normal distribution (the region beyond the provided value)
 {-# INLINE normalTail #-}
-normalTail :: (Distribution StdUniform a, Floating a, Ord a) =>
+normalTail :: (Distribution StdUniform a, Floating a, Ord a, Functor m) =>
               a -> RVarT m a
 normalTail r = go
     where
@@ -97,7 +98,7 @@ normalTail r = go
 -- @logBase 2 c@ and the 'zGetIU' implementation.
 normalZ ::
   (RealFloat a, Erf a, Vector v a, Distribution Uniform a, Integral b) =>
-  b -> (forall m. RVarT m (Int, a)) -> Ziggurat v a
+  b -> (forall m. Functor m => RVarT m (Int, a)) -> Ziggurat v a
 normalZ p = mkZigguratRec True normalF normalFInv normalFInt normalFVol (2^p)
 
 -- | Ziggurat target function (upper half of a non-normalized gaussian PDF)
@@ -133,21 +134,21 @@ normalFVol = sqrt (0.5 * pi)
 --
 -- As far as I know, this should be safe to use in any monomorphic
 -- @Distribution Normal@ instance declaration.
-realFloatStdNormal :: (RealFloat a, Erf a, Distribution Uniform a) => RVarT m a
+realFloatStdNormal :: (RealFloat a, Erf a, Distribution Uniform a, Functor m) => RVarT m a
 realFloatStdNormal = runZiggurat (normalZ p getIU `asTypeOf` (undefined :: Ziggurat V.Vector a))
     where
         p :: Int
         p = 6
 
-        getIU :: (Num a, Distribution Uniform a) => RVarT m (Int, a)
+        getIU :: (Num a, Distribution Uniform a, Functor m) => RVarT m (Int, a)
         getIU = do
-            i <- getRandomWord8
+            i <- Random.uniformWord8 RGen
             u <- uniformT (-1) 1
             return (fromIntegral i .&. (2^p-1), u)
 
 -- |A random variable sampling from the standard normal distribution
 -- over the 'Double' type.
-doubleStdNormal :: RVarT m Double
+doubleStdNormal :: Functor m => RVarT m Double
 doubleStdNormal = runZiggurat doubleStdNormalZ
 
 -- doubleStdNormalC must not be over 2^12 if using wordToDoubleWithExcess
@@ -165,15 +166,29 @@ doubleStdNormalZ = mkZiggurat_ True
         getIU
         (normalTail doubleStdNormalR)
     where
-        getIU :: RVarT m (Int, Double)
+        getIU :: Functor m => RVarT m (Int, Double)
         getIU = do
-            !w <- getRandomWord64
+            !w <- Random.uniformWord64 RGen
             let (u,i) = wordToDoubleWithExcess w
             return $! (fromIntegral i .&. (doubleStdNormalC-1), u+u-1)
 
+-- NOTE: inlined from random-source
+{-# INLINE wordToDouble #-}
+-- |Pack the low 52 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 .&. 0x000fffffffffffff {- 2^52-1 -})) $ (-52)
+
+{-# INLINE wordToDoubleWithExcess #-}
+-- |Same as wordToDouble, but also return the unused bits (as the 12
+-- least significant bits of a 'Word64')
+wordToDoubleWithExcess :: Word64 -> (Double, Word64)
+wordToDoubleWithExcess x = (wordToDouble x, x `shiftR` 52)
+
+
 -- |A random variable sampling from the standard normal distribution
 -- over the 'Float' type.
-floatStdNormal :: RVarT m Float
+floatStdNormal :: Functor m => RVarT m Float
 floatStdNormal = runZiggurat floatStdNormalZ
 
 -- floatStdNormalC must not be over 2^9 if using word32ToFloatWithExcess
@@ -191,12 +206,26 @@ floatStdNormalZ = mkZiggurat_ True
         getIU
         (normalTail floatStdNormalR)
     where
-        getIU :: RVarT m (Int, Float)
+        getIU :: Functor m => RVarT m (Int, Float)
         getIU = do
-            !w <- getRandomWord32
+            !w <- Random.uniformWord32 RGen
             let (u,i) = word32ToFloatWithExcess w
             return (fromIntegral i .&. (floatStdNormalC-1), u+u-1)
 
+-- NOTE: inlined from random-source
+{-# INLINE word32ToFloat #-}
+-- |Pack the low 23 bits from a 'Word32' into a 'Float' in the range [0,1).
+-- Used to convert a 'stdUniform' 'Word32' to a 'stdUniform' 'Double'.
+word32ToFloat :: Word32 -> Float
+word32ToFloat x = (encodeFloat $! toInteger (x .&. 0x007fffff {- 2^23-1 -} )) $ (-23)
+
+{-# INLINE word32ToFloatWithExcess #-}
+-- |Same as word32ToFloat, but also return the unused bits (as the 9
+-- least significant bits of a 'Word32')
+word32ToFloatWithExcess :: Word32 -> (Float, Word32)
+word32ToFloatWithExcess x = (word32ToFloat x, x `shiftR` 23)
+
+
 normalCdf :: (Real a) => a -> a -> a -> Double
 normalCdf m s x = normcdf ((realToFrac x - realToFrac m) / realToFrac s)
 
@@ -249,7 +278,7 @@ stdNormal :: Distribution Normal a => RVar a
 stdNormal = rvar StdNormal
 
 -- |'stdNormalT' is a normal process with distribution 'StdNormal'.
-stdNormalT :: Distribution Normal a => RVarT m a
+stdNormalT :: (Distribution Normal a, Functor m) => RVarT m a
 stdNormalT = rvarT StdNormal
 
 -- |@normal m s@ is a random variable with distribution @'Normal' m s@.
@@ -257,5 +286,5 @@ normal :: Distribution Normal a => a -> a -> RVar a
 normal m s = rvar (Normal m s)
 
 -- |@normalT m s@ is a random process with distribution @'Normal' m s@.
-normalT :: Distribution Normal a => a -> a -> RVarT m a
+normalT :: (Distribution Normal a, Functor m) => a -> a -> RVarT m a
 normalT m s = rvarT (Normal m s)
diff --git a/random-fu/src/Data/Random/Distribution/Pareto.hs b/random-fu/src/Data/Random/Distribution/Pareto.hs
index 0ad85e7..d81078f 100644
--- a/random-fu/src/Data/Random/Distribution/Pareto.hs
+++ b/random-fu/src/Data/Random/Distribution/Pareto.hs
@@ -12,7 +12,7 @@ import Data.Random
 pareto :: Distribution Pareto a => a -> a -> RVar a
 pareto xM a = rvar (Pareto xM a)
 
-paretoT :: Distribution Pareto a => a -> a -> RVarT m a
+paretoT :: (Distribution Pareto a, Functor m) => a -> a -> RVarT m a
 paretoT xM a = rvarT (Pareto xM a)
 
 data Pareto a = Pareto !a !a
diff --git a/random-fu/src/Data/Random/Distribution/Poisson.hs b/random-fu/src/Data/Random/Distribution/Poisson.hs
index 81bbec6..fa7b828 100644
--- a/random-fu/src/Data/Random/Distribution/Poisson.hs
+++ b/random-fu/src/Data/Random/Distribution/Poisson.hs
@@ -19,10 +19,10 @@ import Data.Random.Distribution.Binomial
 import Control.Monad
 
 -- from Knuth, with interpretation help from gsl sources
-integralPoisson :: (Integral a, RealFloat b, Distribution StdUniform b, Distribution (Erlang a) b, Distribution (Binomial b) a) => b -> RVarT m a
+integralPoisson :: (Integral a, RealFloat b, Distribution StdUniform b, Distribution (Erlang a) b, Distribution (Binomial b) a, Functor m) => b -> RVarT m a
 integralPoisson = psn 0
     where
-        psn :: (Integral a, RealFloat b, Distribution StdUniform b, Distribution (Erlang a) b, Distribution (Binomial b) a) => a -> b -> RVarT m a
+        psn :: (Integral a, RealFloat b, Distribution StdUniform b, Distribution (Erlang a) b, Distribution (Binomial b) a, Functor m) => a -> b -> RVarT m a
         psn j mu
             | mu > 10   = do
                 let m = floor (mu * (7/8))
@@ -70,7 +70,7 @@ integralPoissonPDF mu k = exp (negate lambda) *
     k_fac_ln = foldl (+) 0 (map (log . fromIntegral) [1..k])
     lambda   = realToFrac mu
 
-fractionalPoisson :: (Num a, Distribution (Poisson b) Integer) => b -> RVarT m a
+fractionalPoisson :: (Num a, Distribution (Poisson b) Integer, Functor m) => b -> RVarT m a
 fractionalPoisson mu = liftM fromInteger (poissonT mu)
 
 fractionalPoissonCDF :: (CDF (Poisson b) Integer, RealFrac a) => b -> a -> Double
@@ -82,7 +82,7 @@ fractionalPoissonPDF mu k = pdf (Poisson mu) (floor k :: Integer)
 poisson :: (Distribution (Poisson b) a) => b -> RVar a
 poisson mu = rvar (Poisson mu)
 
-poissonT :: (Distribution (Poisson b) a) => b -> RVarT m a
+poissonT :: (Distribution (Poisson b) a, Functor m) => b -> RVarT m a
 poissonT mu = rvarT (Poisson mu)
 
 newtype Poisson b a = Poisson b
diff --git a/random-fu/src/Data/Random/Distribution/Rayleigh.hs b/random-fu/src/Data/Random/Distribution/Rayleigh.hs
index 94ab61f..0eab82b 100644
--- a/random-fu/src/Data/Random/Distribution/Rayleigh.hs
+++ b/random-fu/src/Data/Random/Distribution/Rayleigh.hs
@@ -12,7 +12,7 @@ import Data.Random.RVar
 import Data.Random.Distribution
 import Data.Random.Distribution.Uniform
 
-floatingRayleigh :: (Floating a, Eq a, Distribution StdUniform a) => a -> RVarT m a
+floatingRayleigh :: (Floating a, Eq a, Distribution StdUniform a, Functor m) => a -> RVarT m a
 floatingRayleigh s = do
     u <- stdUniformPosT
     return (s * sqrt (-2 * log u))
@@ -26,7 +26,7 @@ newtype Rayleigh a = Rayleigh a
 rayleigh :: Distribution Rayleigh a => a -> RVar a
 rayleigh = rvar . Rayleigh
 
-rayleighT :: Distribution Rayleigh a => a -> RVarT m a
+rayleighT :: (Distribution Rayleigh a, Functor m) => a -> RVarT m a
 rayleighT = rvarT . Rayleigh
 
 rayleighCDF :: Real a => a -> a -> Double
diff --git a/random-fu/src/Data/Random/Distribution/Simplex.hs b/random-fu/src/Data/Random/Distribution/Simplex.hs
index 3e1c0fd..bf5044b 100644
--- a/random-fu/src/Data/Random/Distribution/Simplex.hs
+++ b/random-fu/src/Data/Random/Distribution/Simplex.hs
@@ -35,7 +35,7 @@ instance (Ord a, Fractional a, Distribution StdUniform a) => Distribution StdSim
 stdSimplex :: Distribution StdSimplex [a] => Int -> RVar [a]
 stdSimplex k = rvar (StdSimplex k)
 
-stdSimplexT :: Distribution StdSimplex [a] => Int -> RVarT m [a]
+stdSimplexT :: (Distribution StdSimplex [a], Functor m) => Int -> RVarT m [a]
 stdSimplexT k = rvarT (StdSimplex k)
 
 -- |An algorithm proposed by Rubinstein & Melamed (1998).
diff --git a/random-fu/src/Data/Random/Distribution/StretchedExponential.hs b/random-fu/src/Data/Random/Distribution/StretchedExponential.hs
index f897dd2..55c96ad 100644
--- a/random-fu/src/Data/Random/Distribution/StretchedExponential.hs
+++ b/random-fu/src/Data/Random/Distribution/StretchedExponential.hs
@@ -14,7 +14,7 @@ import Data.Random.Distribution.Uniform
 
 newtype StretchedExponential a = StretchedExp (a,a)
 
-floatingStretchedExponential :: (Floating a, Distribution StdUniform a) => a -> a -> RVarT m a
+floatingStretchedExponential :: (Floating a, Distribution StdUniform a, Functor m) => a -> a -> RVarT m a
 floatingStretchedExponential beta lambdaRecip = do
     x <- stdUniformT
     return (negate (log (1-x))**(1/beta) * lambdaRecip)
@@ -25,7 +25,7 @@ floatingStretchedExponentialCDF beta lambdaRecip x = 1 - exp (negate (realToFrac
 stretchedExponential :: Distribution StretchedExponential a => a -> a -> RVar a
 stretchedExponential beta lambdaRecip = rvar $ StretchedExp (beta, lambdaRecip)
 
-stretchedExponentialT :: Distribution StretchedExponential a => a -> a -> RVarT m a
+stretchedExponentialT :: (Distribution StretchedExponential a, Functor m) => a -> a -> RVarT m a
 stretchedExponentialT beta lambdaRecip = rvarT $ StretchedExp (beta, lambdaRecip)
 
 instance (Floating a, Distribution StdUniform a) => Distribution StretchedExponential a where
diff --git a/random-fu/src/Data/Random/Distribution/T.hs b/random-fu/src/Data/Random/Distribution/T.hs
index 75cb171..47efc12 100644
--- a/random-fu/src/Data/Random/Distribution/T.hs
+++ b/random-fu/src/Data/Random/Distribution/T.hs
@@ -24,7 +24,7 @@ import Numeric.SpecFunctions
 t :: Distribution T a => Integer -> RVar a
 t = rvar . T
 
-tT :: Distribution T a => Integer -> RVarT m a
+tT :: (Distribution T a, Functor m) => Integer -> RVarT m a
 tT = rvarT . T
 
 newtype T a = T Integer
diff --git a/random-fu/src/Data/Random/Distribution/Triangular.hs b/random-fu/src/Data/Random/Distribution/Triangular.hs
index 272c5a7..88fa96a 100644
--- a/random-fu/src/Data/Random/Distribution/Triangular.hs
+++ b/random-fu/src/Data/Random/Distribution/Triangular.hs
@@ -26,7 +26,7 @@ data Triangular a = Triangular {
     deriving (Eq, Show)
 
 -- |Compute a triangular distribution for a 'Floating' type.
-floatingTriangular :: (Floating a, Ord a, Distribution StdUniform a) => a -> a -> a -> RVarT m a
+floatingTriangular :: (Floating a, Ord a, Distribution StdUniform a, Functor m) => a -> a -> a -> RVarT m a
 floatingTriangular a b c
     | a > b     = floatingTriangular b a c
     | b > c     = floatingTriangular a c b
diff --git a/random-fu/src/Data/Random/Distribution/Uniform.hs b/random-fu/src/Data/Random/Distribution/Uniform.hs
index a580fe8..29b414a 100644
--- a/random-fu/src/Data/Random/Distribution/Uniform.hs
+++ b/random-fu/src/Data/Random/Distribution/Uniform.hs
@@ -38,10 +38,8 @@ module Data.Random.Distribution.Uniform
     ) where
 
 import Data.Random.Internal.TH
-import Data.Random.Internal.Words
 import Data.Random.Internal.Fixed
 
-import Data.Random.Source
 import Data.Random.Distribution
 import Data.Random.RVar
 
@@ -51,39 +49,14 @@ import Data.Int
 
 import Control.Monad.Loops
 
+import qualified System.Random.Stateful as Random
+
 -- |Compute a random 'Integral' value between the 2 values provided (inclusive).
 {-# INLINE integralUniform #-}
-integralUniform :: (Integral a) => a -> a -> RVarT m a
-integralUniform !x !y = if x < y then integralUniform' x y else integralUniform' y x
-
-{-# SPECIALIZE integralUniform' :: Int     -> Int     -> RVarT m Int   #-}
-{-# SPECIALIZE integralUniform' :: Int8    -> Int8    -> RVarT m Int8  #-}
-{-# SPECIALIZE integralUniform' :: Int16   -> Int16   -> RVarT m Int16 #-}
-{-# SPECIALIZE integralUniform' :: Int32   -> Int32   -> RVarT m Int32 #-}
-{-# SPECIALIZE integralUniform' :: Int64   -> Int64   -> RVarT m Int64 #-}
-{-# SPECIALIZE integralUniform' :: Word    -> Word    -> RVarT m Word   #-}
-{-# SPECIALIZE integralUniform' :: Word8   -> Word8   -> RVarT m Word8  #-}
-{-# SPECIALIZE integralUniform' :: Word16  -> Word16  -> RVarT m Word16 #-}
-{-# SPECIALIZE integralUniform' :: Word32  -> Word32  -> RVarT m Word32 #-}
-{-# SPECIALIZE integralUniform' :: Word64  -> Word64  -> RVarT m Word64 #-}
-{-# SPECIALIZE integralUniform' :: Integer -> Integer -> RVarT m Integer #-}
-integralUniform' :: (Integral a) => a -> a -> RVarT m a
-integralUniform' !l !u
-    | nReject == 0  = fmap shift prim
-    | otherwise     = fmap shift loop
-    where
-        m = 1 + toInteger u - toInteger l
-        (bytes, nPossible) = bytesNeeded m
-        nReject = nPossible `mod` m
-
-        !prim = getRandomNByteInteger bytes
-        !shift = \(!z) -> l + (fromInteger $! (z `mod` m))
-
-        loop = do
-            z <- prim
-            if z < nReject
-                then loop
-                else return z
+integralUniform :: Functor m => Random.UniformRange a => a -> a -> RVarT m a
+integralUniform !x !y = Random.uniformRM (x, y) RGen
+  -- Maybe switch to uniformIntegralM (requires exposing from `random` internals):
+  -- Random.uniformIntegralM (x, y) RGen
 
 integralUniformCDF :: (Integral a, Fractional b) => a -> a -> a -> b
 integralUniformCDF a b x
@@ -92,15 +65,6 @@ integralUniformCDF a b x
     | x > b     = 1
     | otherwise = (fromIntegral x - fromIntegral a) / (fromIntegral b - fromIntegral a)
 
--- TODO: come up with a decent, fast heuristic to decide whether to return an extra
--- byte.  May involve moving calculation of nReject into this function, and then
--- accepting first if 4*nReject < nPossible or something similar.
-bytesNeeded :: Integer -> (Int, Integer)
-bytesNeeded x = head (dropWhile ((<= x).snd) powersOf256)
-
-powersOf256 :: [(Int, Integer)]
-powersOf256 = zip [0..] (iterate (256 *) 1)
-
 -- |Compute a random value for a 'Bounded' type, between 'minBound' and 'maxBound'
 -- (inclusive for 'Integral' or 'Enum' types, in ['minBound', 'maxBound') for Fractional types.)
 boundedStdUniform :: (Distribution Uniform a, Bounded a) => RVar a
@@ -111,25 +75,29 @@ boundedStdUniformCDF = cdf (Uniform minBound maxBound)
 
 -- |Compute a random value for a 'Bounded' 'Enum' type, between 'minBound' and
 -- 'maxBound' (inclusive)
-boundedEnumStdUniform :: (Enum a, Bounded a) => RVarT m a
+boundedEnumStdUniform :: Functor m => (Enum a, Bounded a) => RVarT m a
 boundedEnumStdUniform = enumUniform minBound maxBound
 
 boundedEnumStdUniformCDF :: (Enum a, Bounded a, Ord a) => a -> Double
 boundedEnumStdUniformCDF = enumUniformCDF minBound maxBound
 
 -- |Compute a uniform random 'Float' value in the range [0,1)
-floatStdUniform :: RVarT m Float
+floatStdUniform :: Functor m => RVarT m Float
 floatStdUniform = do
-    x <- getRandomWord32
-    return (word32ToFloat x)
+    x <- uniformRangeRVarT (0, 1)
+    -- exclude 1. TODO: come up with something smarter
+    if x == 1 then floatStdUniform else pure x
 
 -- |Compute a uniform random 'Double' value in the range [0,1)
 {-# INLINE doubleStdUniform #-}
-doubleStdUniform :: RVarT m Double
-doubleStdUniform = getRandomDouble
+doubleStdUniform :: Functor m => RVarT m Double
+doubleStdUniform = do
+    x <- uniformRangeRVarT (0, 1)
+    -- exclude 1. TODO: come up with something smarter
+    if x == 1 then doubleStdUniform else pure x
 
 -- |Compute a uniform random value in the range [0,1) for any 'RealFloat' type
-realFloatStdUniform :: RealFloat a => RVarT m a
+realFloatStdUniform :: Functor m => RealFloat a => RVarT m a
 realFloatStdUniform = do
     let (b, e) = decodeFloat one
 
@@ -142,7 +110,7 @@ realFloatStdUniform = do
 
 -- |Compute a uniform random 'Fixed' value in the range [0,1), with any
 -- desired precision.
-fixedStdUniform :: HasResolution r => RVarT m (Fixed r)
+fixedStdUniform :: Functor m => HasResolution r => RVarT m (Fixed r)
 fixedStdUniform = x
     where
         res = resolutionOf2 x
@@ -170,7 +138,7 @@ lerp :: Num a => a -> a -> a -> a
 lerp x y a = (1-a)*x + a*y
 
 -- |@floatUniform a b@ computes a uniform random 'Float' value in the range [a,b)
-floatUniform :: Float -> Float -> RVarT m Float
+floatUniform :: Functor m => Float -> Float -> RVarT m Float
 floatUniform 0 1 = floatStdUniform
 floatUniform a b = do
     x <- floatStdUniform
@@ -178,7 +146,7 @@ floatUniform a b = do
 
 -- |@doubleUniform a b@ computes a uniform random 'Double' value in the range [a,b)
 {-# INLINE doubleUniform #-}
-doubleUniform :: Double -> Double -> RVarT m Double
+doubleUniform :: Functor m => Double -> Double -> RVarT m Double
 doubleUniform 0 1 = doubleStdUniform
 doubleUniform a b = do
     x <- doubleStdUniform
@@ -186,7 +154,7 @@ doubleUniform a b = do
 
 -- |@realFloatUniform a b@ computes a uniform random value in the range [a,b) for
 -- any 'RealFloat' type
-realFloatUniform :: RealFloat a => a -> a -> RVarT m a
+realFloatUniform :: Functor m => RealFloat a => a -> a -> RVarT m a
 realFloatUniform 0 1 = realFloatStdUniform
 realFloatUniform a b = do
     x <- realFloatStdUniform
@@ -194,7 +162,7 @@ realFloatUniform a b = do
 
 -- |@fixedUniform a b@ computes a uniform random 'Fixed' value in the range
 -- [a,b), with any desired precision.
-fixedUniform :: HasResolution r => Fixed r -> Fixed r -> RVarT m (Fixed r)
+fixedUniform :: Functor m => HasResolution r => Fixed r -> Fixed r -> RVarT m (Fixed r)
 fixedUniform a b = do
     u <- integralUniform (unMkFixed a) (unMkFixed b)
     return (mkFixed u)
@@ -209,7 +177,7 @@ realUniformCDF a b x
 
 -- |@realFloatUniform a b@ computes a uniform random value in the range [a,b) for
 -- any 'Enum' type
-enumUniform :: Enum a => a -> a -> RVarT m a
+enumUniform :: Functor m => Enum a => a -> a -> RVarT m a
 enumUniform a b = do
     x <- integralUniform (fromEnum a) (fromEnum b)
     return (toEnum x)
@@ -232,7 +200,7 @@ uniform a b = rvar (Uniform a b)
 -- @uniformT a b@ is a uniformly distributed random process in the range
 -- [a,b] for 'Integral' or 'Enum' types and in the range [a,b) for 'Fractional'
 -- types.  Requires a @Distribution Uniform@ instance for the type.
-uniformT :: Distribution Uniform a => a -> a -> RVarT m a
+uniformT :: (Distribution Uniform a, Functor m) => a -> a -> RVarT m a
 uniformT a b = rvarT (Uniform a b)
 
 -- |Get a \"standard\" uniformly distributed variable.
@@ -248,15 +216,15 @@ stdUniform = rvar StdUniform
 -- For integral types, this means uniformly distributed over the full range
 -- of the type (there is no support for 'Integer').  For fractional
 -- types, this means uniformly distributed on the interval [0,1).
-{-# SPECIALIZE stdUniformT :: RVarT m Double #-}
-{-# SPECIALIZE stdUniformT :: RVarT m Float #-}
-stdUniformT :: (Distribution StdUniform a) => RVarT m a
+{-# SPECIALIZE stdUniformT :: Functor m => RVarT m Double #-}
+{-# SPECIALIZE stdUniformT :: Functor m => RVarT m Float #-}
+stdUniformT :: (Distribution StdUniform a, Functor m) => RVarT m a
 stdUniformT = rvarT StdUniform
 
 -- |Like 'stdUniform', but returns only positive or zero values.  Not
 -- exported because it is not truly uniform: nonzero values are twice
 -- as likely as zero on signed types.
-stdUniformNonneg :: (Distribution StdUniform a, Num a, Eq a) => RVarT m a
+stdUniformNonneg :: (Distribution StdUniform a, Num a, Eq a, Functor m) => RVarT m a
 stdUniformNonneg = fmap abs stdUniformT
 
 -- |Like 'stdUniform' but only returns positive values.
@@ -264,7 +232,7 @@ stdUniformPos :: (Distribution StdUniform a, Num a, Eq a) => RVar a
 stdUniformPos = stdUniformPosT
 
 -- |Like 'stdUniform' but only returns positive values.
-stdUniformPosT :: (Distribution StdUniform a, Num a, Eq a) => RVarT m a
+stdUniformPosT :: (Distribution StdUniform a, Num a, Eq a, Functor m) => RVarT m a
 stdUniformPosT = iterateUntil (/= 0) stdUniformNonneg
 
 -- |A definition of a uniform distribution over the type @t@.  See also 'uniform'.
@@ -289,27 +257,19 @@ $( replicateInstances ''Int integralTypes [d|
         instance CDF Uniform Int            where cdf   (Uniform a b) = integralUniformCDF a b
     |])
 
-instance Distribution StdUniform Word8      where rvarT _ = getRandomWord8
-instance Distribution StdUniform Word16     where rvarT _ = getRandomWord16
-instance Distribution StdUniform Word32     where rvarT _ = getRandomWord32
-instance Distribution StdUniform Word64     where rvarT _ = getRandomWord64
+instance Distribution StdUniform Word8      where rvarT _ = Random.uniformWord8 RGen
+instance Distribution StdUniform Word16     where rvarT _ = Random.uniformWord16 RGen
+instance Distribution StdUniform Word32     where rvarT _ = Random.uniformWord32 RGen
+instance Distribution StdUniform Word64     where rvarT _ = Random.uniformWord64 RGen
+instance Distribution StdUniform Word       where rvarT _ = uniformRVarT
 
-instance Distribution StdUniform Int8       where rvarT _ = fromIntegral `fmap` getRandomWord8
-instance Distribution StdUniform Int16      where rvarT _ = fromIntegral `fmap` getRandomWord16
-instance Distribution StdUniform Int32      where rvarT _ = fromIntegral `fmap` getRandomWord32
-instance Distribution StdUniform Int64      where rvarT _ = fromIntegral `fmap` getRandomWord64
+instance Distribution StdUniform Int8       where rvarT _ = uniformRVarT
+instance Distribution StdUniform Int16      where rvarT _ = uniformRVarT
+instance Distribution StdUniform Int32      where rvarT _ = uniformRVarT
+instance Distribution StdUniform Int64      where rvarT _ = uniformRVarT
 
-instance Distribution StdUniform Int where
-    rvar _ =
-        $(if toInteger (maxBound :: Int) > toInteger (maxBound :: Int32)
-            then [|fromIntegral `fmap` getRandomWord64 :: RVar Int|]
-            else [|fromIntegral `fmap` getRandomWord32 :: RVar Int|])
+instance Distribution StdUniform Int        where rvarT _ = uniformRVarT
 
-instance Distribution StdUniform Word where
-    rvar _ =
-        $(if toInteger (maxBound :: Word) > toInteger (maxBound :: Word32)
-            then [|fromIntegral `fmap` getRandomWord64 :: RVar Word|]
-            else [|fromIntegral `fmap` getRandomWord32 :: RVar Word|])
 
 -- Integer has no StdUniform...
 
@@ -331,7 +291,7 @@ instance CDF Uniform Float                  where cdf   (Uniform a b) = realUnif
 instance CDF Uniform Double                 where cdf   (Uniform a b) = realUniformCDF a b
 
 instance Distribution StdUniform Float      where rvarT _ = floatStdUniform
-instance Distribution StdUniform Double     where rvarT _ = getRandomDouble
+instance Distribution StdUniform Double     where rvarT _ = uniformRangeRVarT (0, 1)
 instance CDF StdUniform Float               where cdf   _ = realStdUniformCDF
 instance CDF StdUniform Double              where cdf   _ = realStdUniformCDF
 instance PDF StdUniform Float               where pdf   _ = realStdUniformPDF
@@ -357,7 +317,7 @@ $( replicateInstances ''Char [''Char, ''Bool, ''Ordering] [d|
 
 instance Distribution StdUniform ()         where rvarT ~StdUniform = return ()
 instance CDF StdUniform ()                  where cdf   ~StdUniform = return 1
-instance Distribution StdUniform Bool       where rvarT ~StdUniform = fmap even (getRandomWord8)
+instance Distribution StdUniform Bool       where rvarT ~StdUniform = uniformRVarT
 instance CDF StdUniform Bool                where cdf   ~StdUniform = boundedEnumStdUniformCDF
 
 instance Distribution StdUniform Char       where rvarT ~StdUniform = boundedEnumStdUniform
diff --git a/random-fu/src/Data/Random/Distribution/Ziggurat.hs b/random-fu/src/Data/Random/Distribution/Ziggurat.hs
index 3d680d6..3d41451 100644
--- a/random-fu/src/Data/Random/Distribution/Ziggurat.hs
+++ b/random-fu/src/Data/Random/Distribution/Ziggurat.hs
@@ -73,18 +73,18 @@ data Ziggurat v t = Ziggurat {
         -- single random word (64 bits) can be efficiently converted to
         -- a double (using 52 bits) and a bin number (using up to 12 bits),
         -- for example.
-        zGetIU            :: !(forall m. RVarT m (Int, t)),
+        zGetIU            :: !(forall m. Functor m => RVarT m (Int, t)),
 
         -- |The distribution for the final \"virtual\" bin
         -- (the ziggurat algorithm does not handle distributions
         -- that wander off to infinity, so another distribution is needed
         -- to handle the last \"bin\" that stretches to infinity)
-        zTailDist         :: (forall m. RVarT m t),
+        zTailDist         :: (forall m. Functor m => RVarT m t),
 
         -- |A copy of the uniform RVar generator for the base type,
         -- so that @Distribution Uniform t@ is not needed when sampling
         -- from a Ziggurat (makes it a bit more self-contained).
-        zUniform          :: !(forall m. t -> t -> RVarT m t),
+        zUniform          :: !(forall m. Functor m => t -> t -> RVarT m t),
 
         -- |The (one-sided antitone) PDF, not necessarily normalized
         zFunc             :: !(t -> t),
@@ -99,11 +99,11 @@ data Ziggurat v t = Ziggurat {
 
 -- |Sample from the distribution encoded in a 'Ziggurat' data structure.
 {-# INLINE runZiggurat #-}
-{-# SPECIALIZE runZiggurat :: Ziggurat UV.Vector Float  -> RVarT m Float #-}
-{-# SPECIALIZE runZiggurat :: Ziggurat UV.Vector Double -> RVarT m Double #-}
-{-# SPECIALIZE runZiggurat :: Ziggurat  V.Vector Float  -> RVarT m Float #-}
-{-# SPECIALIZE runZiggurat :: Ziggurat  V.Vector Double -> RVarT m Double #-}
-runZiggurat :: (Num a, Ord a, Vector v a) =>
+{-# SPECIALIZE runZiggurat :: Functor m => Ziggurat UV.Vector Float  -> RVarT m Float #-}
+{-# SPECIALIZE runZiggurat :: Functor m => Ziggurat UV.Vector Double -> RVarT m Double #-}
+{-# SPECIALIZE runZiggurat :: Functor m => Ziggurat  V.Vector Float  -> RVarT m Float #-}
+{-# SPECIALIZE runZiggurat :: Functor m => Ziggurat  V.Vector Double -> RVarT m Double #-}
+runZiggurat :: (Num a, Ord a, Vector v a, Functor m) =>
                Ziggurat v a -> RVarT m a
 runZiggurat !Ziggurat{..} = go
     where
@@ -166,10 +166,10 @@ runZiggurat !Ziggurat{..} = go
 --  * an RVar sampling from the tail (the region where x > R)
 --
 {-# INLINE mkZiggurat_ #-}
-{-# SPECIALIZE mkZiggurat_ :: Bool -> (Float  ->  Float) -> (Float  ->  Float) -> Int -> Float  -> Float  -> (forall m. RVarT m (Int,  Float)) -> (forall m. RVarT m Float ) -> Ziggurat UV.Vector Float #-}
-{-# SPECIALIZE mkZiggurat_ :: Bool -> (Double -> Double) -> (Double -> Double) -> Int -> Double -> Double -> (forall m. RVarT m (Int, Double)) -> (forall m. RVarT m Double) -> Ziggurat UV.Vector Double #-}
-{-# SPECIALIZE mkZiggurat_ :: Bool -> (Float  ->  Float) -> (Float  ->  Float) -> Int -> Float  -> Float  -> (forall m. RVarT m (Int,  Float)) -> (forall m. RVarT m Float ) -> Ziggurat V.Vector Float #-}
-{-# SPECIALIZE mkZiggurat_ :: Bool -> (Double -> Double) -> (Double -> Double) -> Int -> Double -> Double -> (forall m. RVarT m (Int, Double)) -> (forall m. RVarT m Double) -> Ziggurat V.Vector Double #-}
+{-# SPECIALIZE mkZiggurat_ :: Bool -> (Float  ->  Float) -> (Float  ->  Float) -> Int -> Float  -> Float  -> (forall m. Functor m => RVarT m (Int,  Float)) -> (forall m. Functor m => RVarT m Float ) -> Ziggurat UV.Vector Float #-}
+{-# SPECIALIZE mkZiggurat_ :: Bool -> (Double -> Double) -> (Double -> Double) -> Int -> Double -> Double -> (forall m. Functor m => RVarT m (Int, Double)) -> (forall m. Functor m => RVarT m Double) -> Ziggurat UV.Vector Double #-}
+{-# SPECIALIZE mkZiggurat_ :: Bool -> (Float  ->  Float) -> (Float  ->  Float) -> Int -> Float  -> Float  -> (forall m. Functor m => RVarT m (Int,  Float)) -> (forall m. Functor m => RVarT m Float ) -> Ziggurat V.Vector Float #-}
+{-# SPECIALIZE mkZiggurat_ :: Bool -> (Double -> Double) -> (Double -> Double) -> Int -> Double -> Double -> (forall m. Functor m => RVarT m (Int, Double)) -> (forall m. Functor m => RVarT m Double) -> Ziggurat V.Vector Double #-}
 mkZiggurat_ :: (RealFloat t, Vector v t,
                Distribution Uniform t) =>
               Bool
@@ -178,8 +178,8 @@ mkZiggurat_ :: (RealFloat t, Vector v t,
               -> Int
               -> t
               -> t
-              -> (forall m. RVarT m (Int, t))
-              -> (forall m. RVarT m t)
+              -> (forall m. Functor m => RVarT m (Int, t))
+              -> (forall m. Functor m => RVarT m t)
               -> Ziggurat v t
 mkZiggurat_ m f fInv c r v getIU tailDist = Ziggurat
     { zTable_xs         = xs
@@ -209,8 +209,8 @@ mkZiggurat :: (RealFloat t, Vector v t,
               -> (t -> t)
               -> t
               -> Int
-              -> (forall m. RVarT m (Int, t))
-              -> (forall m. t -> RVarT m t)
+              -> (forall m. Functor m => RVarT m (Int, t))
+              -> (forall m. Functor m => t -> RVarT m t)
               -> Ziggurat v t
 mkZiggurat m f fInv fInt fVol c getIU tailDist =
     mkZiggurat_ m f fInv c r v getIU (tailDist r)
@@ -246,11 +246,11 @@ mkZigguratRec ::
   -> (t -> t)
   -> t
   -> Int
-  -> (forall m. RVarT m (Int, t))
+  -> (forall m. Functor m => RVarT m (Int, t))
   -> Ziggurat v t
 mkZigguratRec m f fInv fInt fVol c getIU = z
         where
-            fix :: ((forall m. a -> RVarT m a) -> (forall m. a -> RVarT m a)) -> (forall m. a -> RVarT m a)
+            fix :: ((forall m. Functor m => a -> RVarT m a) -> (forall m. Functor m => a -> RVarT m a)) -> (forall m. Functor m => a -> RVarT m a)
             fix g = g (fix g)
             z = mkZiggurat m f fInv fInt fVol c getIU (fix (mkTail m f fInv fInt fVol c getIU z))
 
@@ -260,10 +260,10 @@ mkTail ::
     -> (a -> a) -> (a -> a) -> (a -> a)
     -> a
     -> Int
-    -> (forall m. RVarT m (Int, a))
+    -> (forall m. Functor m => RVarT m (Int, a))
     -> Ziggurat v a
-    -> (forall m. a -> RVarT m a)
-    -> (forall m. a -> RVarT m a)
+    -> (forall m. Functor m => a -> RVarT m a)
+    -> (forall m. Functor m => a -> RVarT m a)
 mkTail m f fInv fInt fVol c getIU typeRep nextTail r = do
      x <- rvarT (mkZiggurat m f' fInv' fInt' fVol' c getIU nextTail `asTypeOf` typeRep)
      return (x + r * signum x)
diff --git a/random-fu/src/Data/Random/Lift.hs b/random-fu/src/Data/Random/Lift.hs
index 658ac72..f765607 100644
--- a/random-fu/src/Data/Random/Lift.hs
+++ b/random-fu/src/Data/Random/Lift.hs
@@ -5,7 +5,6 @@ module Data.Random.Lift where
 import Data.RVar
 import qualified Data.Functor.Identity as T
 import qualified Control.Monad.Trans.Class as T
-import Data.Random.Source.Std
 
 #ifndef MTL2
 import qualified Control.Monad.Identity as MTL
@@ -40,8 +39,8 @@ instance Lift m m where
 instance Monad m => Lift T.Identity m where
     lift = return . T.runIdentity
 
-instance Lift (RVarT T.Identity) (RVarT m) where
-    lift x = runRVar x StdRandom
+instance Functor m => Lift (RVarT T.Identity) (RVarT m) where
+    lift x = runRVar x RGen
 
 -- | This instance is again incoherent with the others, but provides a
 -- more-specific instance to resolve the overlap between the
@@ -49,21 +48,3 @@ instance Lift (RVarT T.Identity) (RVarT m) where
 instance T.MonadTrans t => Lift T.Identity (t T.Identity) where
     lift = T.lift
 
-#ifndef MTL2
-
--- | This instance is incoherent with the others. However,
--- by the law @lift (return x) == return x@, the results
--- must always be the same.
-instance Monad m => Lift MTL.Identity m where
-    lift = return . MTL.runIdentity
-
-instance Lift (RVarT MTL.Identity) (RVarT m) where
-    lift x = runRVarTWith (return . MTL.runIdentity) x StdRandom
-
--- | This instance is again incoherent with the others, but provides a
--- more-specific instance to resolve the overlap between the
--- @Lift m (t m)@ and @Lift Identity m@ instances.
-instance T.MonadTrans t => Lift MTL.Identity (t MTL.Identity) where
-    lift = T.lift
-
-#endif
diff --git a/random-fu/src/Data/Random/List.hs b/random-fu/src/Data/Random/List.hs
index 1770fac..607157b 100644
--- a/random-fu/src/Data/Random/List.hs
+++ b/random-fu/src/Data/Random/List.hs
@@ -13,7 +13,7 @@ randomElement :: [a] -> RVar a
 randomElement = randomElementT
 
 
-randomElementT :: [a] -> RVarT m a
+randomElementT :: Functor m => [a] -> RVarT m a
 randomElementT [] = error "randomElementT: empty list!"
 randomElementT xs = do
     n <- uniformT 0 (length xs - 1)
@@ -24,7 +24,7 @@ randomElementT xs = do
 shuffle :: [a] -> RVar [a]
 shuffle = shuffleT
 
-shuffleT :: [a] -> RVarT m [a]
+shuffleT :: Functor m => [a] -> RVarT m [a]
 shuffleT [] = return []
 shuffleT xs = do
     is <- zipWithM (\_ i -> uniformT 0 i) (tail xs) [1..]
@@ -38,14 +38,14 @@ shuffleT xs = do
 shuffleN :: Int -> [a] -> RVar [a]
 shuffleN = shuffleNT
 
-shuffleNT :: Int -> [a] -> RVarT m [a]
+shuffleNT :: Functor m => Int -> [a] -> RVarT m [a]
 shuffleNT n xs = shuffleNofMT n n xs
 
 -- | A random variable that selects N arbitrary elements of a list of known length M.
 shuffleNofM :: Int -> Int -> [a] -> RVar [a]
 shuffleNofM = shuffleNofMT
 
-shuffleNofMT :: Int -> Int -> [a] -> RVarT m [a]
+shuffleNofMT :: Functor m => Int -> Int -> [a] -> RVarT m [a]
 shuffleNofMT 0 _ _ = return []
 shuffleNofMT n m xs
     | n > m     = error "shuffleNofMT: n > m"
diff --git a/random-fu/src/Data/Random/RVar.hs b/random-fu/src/Data/Random/RVar.hs
index 1e69672..029ff49 100644
--- a/random-fu/src/Data/Random/RVar.hs
+++ b/random-fu/src/Data/Random/RVar.hs
@@ -2,13 +2,14 @@
 module Data.Random.RVar
     ( RVar, runRVar
     , RVarT, runRVarT, runRVarTWith
+    , RGen(..), uniformRVarT, uniformRangeRVarT
     ) where
 
 import Data.Random.Lift
-import Data.Random.Internal.Source
 import Data.RVar hiding (runRVarT)
+import System.Random.Stateful
 
 -- |Like 'runRVarTWith', but using an implicit lifting (provided by the
 -- 'Lift' class)
-runRVarT :: (Lift n m, RandomSource m s) => RVarT n a -> s -> m a
+runRVarT :: (Lift n m, StatefulGen g m) => RVarT n a -> g -> m a
 runRVarT = runRVarTWith lift
diff --git a/random-fu/src/Data/Random/Sample.hs b/random-fu/src/Data/Random/Sample.hs
index c1c730c..b7c0f93 100644
--- a/random-fu/src/Data/Random/Sample.hs
+++ b/random-fu/src/Data/Random/Sample.hs
@@ -9,11 +9,12 @@
 module Data.Random.Sample where
 
 import Control.Monad.State
+import Control.Monad.Reader
 import Data.Random.Distribution
 import Data.Random.Lift
 import Data.Random.RVar
-import Data.Random.Source
-import Data.Random.Source.Std
+
+import System.Random.Stateful
 
 -- |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
@@ -21,26 +22,29 @@ import Data.Random.Source.Std
 -- separate abstractions for one base concept: a random variable.
 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
+    sampleFrom :: StatefulGen g m => g -> d t -> m t
 
 instance Distribution d t => Sampleable d m t where
-    sampleFrom src d = runRVarT (rvar d) src
+    sampleFrom gen d = runRVarT (rvar d) gen
 
 -- This instance overlaps 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
+    sampleFrom gen x = runRVarT x gen
 
 -- |Sample a random variable 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
+sample :: (Sampleable d m t, StatefulGen g m, MonadReader g m) => d t -> m t
+sample thing = ask >>= \gen -> sampleFrom gen thing
 
 -- |Sample a random variable in a \"functional\" style.  Typical instantiations
 -- of @s@ are @System.Random.StdGen@ or @System.Random.Mersenne.Pure64.PureMT@.
-sampleState :: (Sampleable d (State s) t, MonadRandom (State s)) => d t -> s -> (t, s)
-sampleState thing = runState (sample thing)
+-- sample :: (Distribution d a, StatefulGen g m, MonadReader g m) => d t -> m t
+-- sample thing gen = runStateGen gen (\stateGen -> sampleFrom stateGen thing)
 
--- |Sample a random variable in a \"semi-functional\" style.  Typical instantiations
--- of @s@ are @System.Random.StdGen@ or @System.Random.Mersenne.Pure64.PureMT@.
-sampleStateT :: (Sampleable d (StateT s m) t, MonadRandom (StateT s m)) => d t -> s -> m (t, s)
-sampleStateT thing = runStateT (sample thing)
+sampleState :: (Distribution d t, RandomGen g, MonadState g m) => d t -> m t
+sampleState thing = sampleFrom StateGenM thing
+
+-- |Sample a random variable in a \"functional\" style.  Typical instantiations
+-- of @g@ are @System.Random.StdGen@ or @System.Random.Mersenne.Pure64.PureMT@.
+samplePure :: (Distribution d t, RandomGen g) => d t -> g -> (t, g)
+samplePure thing gen = runStateGen gen (\stateGen -> sampleFrom stateGen thing)
diff --git a/random-source/src/Data/Random/Source/Internal/Prim.hs b/random-source/src/Data/Random/Source/Internal/Prim.hs
index e4bc67e..8d060d4 100644
--- a/random-source/src/Data/Random/Source/Internal/Prim.hs
+++ b/random-source/src/Data/Random/Source/Internal/Prim.hs
@@ -9,7 +9,7 @@ import Data.Typeable
 
 -- |A 'Prompt' GADT describing a request for a primitive random variate.
 -- Random variable definitions will request their entropy via these prompts,
--- and entropy sources will satisfy those requests.  The functions in
+-- and entropy sources will satisfy those requests.  The 'Integer' 0 <= U < 256^n'Integer' 0 <= U < 256^n'Integer' 0 <= U < 256^nfunctions in
 -- "Data.Random.Source.Internal.TH" extend incomplete entropy-source definitions
 -- to complete ones, essentially defining a very flexible 
 -- implementation-defaulting system.
diff --git a/rvar/changelog.md b/rvar/changelog.md
index 8d02357..e9b65c2 100644
--- a/rvar/changelog.md
+++ b/rvar/changelog.md
@@ -1,3 +1,9 @@
+* Changes in 0.3.0.0:
+
+  * Drop usage of `random-source` in favor of `random`
+  * Add `Prim` type that resembles one from `random-source`
+  * Add `RGen` type that serves the same purpose as `StdRandom` in `random-source`
+
 * Changes in 0.2.0.6: None. (Pacify Hackage.)
 
 * Changes in 0.2.0.4: Update for GHC 8.8.
diff --git a/rvar/rvar.cabal b/rvar/rvar.cabal
index 34d110b..e590788 100644
--- a/rvar/rvar.cabal
+++ b/rvar/rvar.cabal
@@ -1,5 +1,5 @@
 name:                   rvar
-version:                0.2.0.6
+version:                0.3.0.0
 stability:              stable
 
 cabal-version:          >= 1.10
@@ -45,6 +45,7 @@ Library
   hs-source-dirs:       src
   default-language:     Haskell2010
   exposed-modules:      Data.RVar
+  other-modules:        Data.RVar.Prim
 
   if flag(mtl2)
     build-depends:      mtl == 2.*
@@ -53,6 +54,14 @@ Library
     build-depends:      mtl == 1.1.*
 
   build-depends:        base            >= 3 && <5,
+
+  -- Begin for testing
+                        vector,
+                        mwc-random,
+  -- End for testing
+
+                        bytestring,
+                        free,
                         MonadPrompt     == 1.0.*,
-                        random-source   == 0.3.*,
-                        transformers    >= 0.2 && < 0.6
+                        transformers    >= 0.2 && < 0.6,
+                        random          >= 1.2.0
\ No newline at end of file
diff --git a/rvar/src/Data/RVar.hs b/rvar/src/Data/RVar.hs
index cdcba6c..0eb43fa 100644
--- a/rvar/src/Data/RVar.hs
+++ b/rvar/src/Data/RVar.hs
@@ -1,49 +1,54 @@
 {-
  -      ``Data/Random/RVar''
  -}
-{-# LANGUAGE
-    RankNTypes,
-    MultiParamTypeClasses,
-    FlexibleInstances,
-    GADTs,
-    ScopedTypeVariables,
-    CPP
-  #-}
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeOperators #-}
 
 -- |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.
+{-# LANGUAGE FlexibleContexts #-}
+
 module Data.RVar
-    ( RandomSource
-    , MonadRandom
-        ( getRandomWord8
-        , getRandomWord16
-        , getRandomWord32
-        , getRandomWord64
-        , getRandomDouble
-        , getRandomNByteInteger
-        )
-
-    , RVar
-    , runRVar, sampleRVar
+    ( RVar
+    , runRVar
+    , sampleReaderRVar
+    , sampleStateRVar
+    , pureRVar
 
     , RVarT
-    , runRVarT, sampleRVarT
-    , runRVarTWith, sampleRVarTWith
-    ) where
+    , runRVarT, sampleReaderRVarT, sampleStateRVarT
+    , runRVarTWith
+    , sampleReaderRVarTWith
+    , sampleStateRVarTWith
 
+    , RGen(..)
+    , uniformRVarT
+    , uniformRangeRVarT
+
+    , Prim(..)
+    ) where
 
-import Data.Random.Internal.Source (Prim(..), MonadRandom(..), RandomSource(..))
-import Data.Random.Source ({-instances-})
 
-import qualified Control.Monad.Trans.Class as T
-import Control.Monad (liftM, ap)
-import Control.Monad.Prompt (MonadPrompt(..), PromptT, runPromptT)
 import qualified Control.Monad.IO.Class as T
-import qualified Control.Monad.Trans as MTL
+import Control.Monad.Reader as MTL
+import Control.Monad.State as MTL
+import qualified Control.Monad.Trans.Class as T
 import qualified Data.Functor.Identity as T
+import Data.RVar.Prim
+import System.Random.Stateful
+
+import Control.Monad.Free.Church
+import Data.Functor.Sum
+
+import Data.Word
 
 -- |An opaque type modeling a \"random variable\" - a value
 -- which depends on the outcome of some random event.  'RVar's
@@ -65,29 +70,42 @@ import qualified Data.Functor.Identity as T
 --
 -- Once defined (in any style), there are several ways to sample 'RVar's:
 --
--- * In a monad, using a 'RandomSource':
---
--- > runRVar (uniform 1 100) DevRandom :: IO Int
+-- * Using an immutable pseudo-random number generator that has an instance for `RandomGen` with
+--   `StateT` monad:
 --
--- * In a monad, using a 'MonadRandom' instance:
+-- >>> import qualified Data.Random as Fu (uniform)
+-- >>> import System.Random (mkStdGen)
+-- >>> import Control.Monad.State (runState)
+-- >>> runState (sampleStateRVar (Fu.uniform 1 (100 :: Integer))) (mkStdGen 2021)
+-- (79,StdGen {unStdGen = SMGen 4687568268719557181 4805600293067301895})
 --
--- > sampleRVar (uniform 1 100) :: State PureMT Int
+-- * Using a mutable pseud-random number generator that has an instance for `StatefulGen` with
+--   `ReaderT` monad.
 --
--- * As a pure function transforming a functional RNG:
+-- >>> import qualified Data.Random as Fu (uniform)
+-- >>> import System.Random.MWC (create)
+-- >>> import Control.Monad.Reader (runReaderT)
+-- >>> import qualified Data.Vector.Storable as VS
+-- >>> initialize (VS.singleton 2021) >>= runReaderT (sampleReaderRVar (uniform 1 (100 :: Integer)))
+-- 8
 --
--- > sampleState (uniform 1 100) :: StdGen -> (Int, StdGen)
---
--- (where @sampleState = runState . sampleRVar@)
 type RVar = RVarT T.Identity
 
+-- | Sample random variable using `RandomGen` generator as source of entropy
+pureRVar :: RandomGen g => RVar a -> g -> (a, g)
+pureRVar rvar g = runStateGen g (runRVar rvar)
+
 -- |\"Run\" an 'RVar' - samples the random variable from the provided
 -- source of entropy.
-runRVar :: RandomSource m s => RVar a -> s -> m a
+runRVar :: StatefulGen g m => RVar a -> g -> m a
 runRVar = runRVarTWith (return . T.runIdentity)
 
 -- |@sampleRVar x@ is equivalent to @runRVar x 'StdRandom'@.
-sampleRVar :: MonadRandom m => RVar a -> m a
-sampleRVar = sampleRVarTWith (return . T.runIdentity)
+sampleReaderRVar :: (StatefulGen g m, MonadReader g m) => RVar a -> m a
+sampleReaderRVar = sampleReaderRVarTWith (return . T.runIdentity)
+
+sampleStateRVar :: (RandomGen g, MonadState g m) => RVar a -> m a
+sampleStateRVar = sampleStateRVarTWith (return . T.runIdentity)
 
 -- |A random variable with access to operations in an underlying monad.  Useful
 -- examples include any form of state for implementing random processes with hysteresis,
@@ -137,13 +155,17 @@ sampleRVar = sampleRVarTWith (return . T.runIdentity)
 -- >         flip runStateT gen .
 -- >             sampleRVarTWith MTL.lift $
 -- >                 replicateM count rwalkState
-newtype RVarT m a = RVarT { unRVarT :: PromptT Prim m a }
+newtype RVarT m a = RVarT { unRVarT :: F (Sum Prim m) a }
 
-runRVarT :: RandomSource m s => RVarT m a -> s -> m a
+runRVarT :: StatefulGen g m => RVarT m a -> g -> m a
 runRVarT = runRVarTWith id
 
-sampleRVarT :: MonadRandom m => RVarT m a -> m a
-sampleRVarT = sampleRVarTWith id
+
+sampleStateRVarT :: (RandomGen g, MonadState g m) => RVarT m a -> m a
+sampleStateRVarT rvar = runRVarT rvar StateGenM
+
+sampleReaderRVarT :: (StatefulGen g m, MonadReader g m) => RVarT m a -> m a
+sampleReaderRVarT rvar = ask >>= runRVarT rvar
 
 -- | \"Runs\" an 'RVarT', sampling the random variable it defines.
 --
@@ -171,24 +193,40 @@ sampleRVarT = sampleRVarTWith id
 -- or functions manipulating 'RVar's would have to use higher-ranked
 -- types to enforce the same kind of isolation and polymorphism.
 {-# INLINE runRVarTWith #-}
-runRVarTWith :: forall m n s a. RandomSource m s => (forall t. n t -> m t) -> RVarT n a -> s -> m a
-runRVarTWith liftN (RVarT m) src = runPromptT return bindP bindN m
-    where
-        bindP :: forall t. (Prim t -> (t -> m a) -> m a)
-        bindP prim cont = getRandomPrimFrom src prim >>= cont
+runRVarTWith :: forall m n g a . StatefulGen g m =>
+           (n :-> m) -> RVarT n a -> g -> m a
+runRVarTWith liftN (RVarT m) gen = runPromptT' (\p -> uniformPrimM p gen) liftN m
 
-        bindN :: forall t. n t -> (t -> m a) -> m a
-        bindN nExp cont = liftN nExp >>= cont
+type f :-> g = forall i. f i -> g i
+
+runPromptT' :: forall m n p r . (Monad n) => (p :-> n) -> (m :-> n) -> F (Sum p m) r -> n r
+runPromptT' prm lft e = foldF alg e
+  where
+    alg :: Sum p m i -> n i
+    alg (InL x) = prm x
+    alg (InR y) = lft y
+
+{-# INLINE uniformPrimM #-}
+uniformPrimM :: StatefulGen g m => Prim t -> g -> m t
+uniformPrimM (Prim f) g = uniformWord32 g >>= return . f
 
 -- |@sampleRVarTWith lift x@ is equivalent to @runRVarTWith lift x 'StdRandom'@.
-sampleRVarTWith :: forall m n a. MonadRandom m => (forall t. n t -> m t) -> RVarT n a -> m a
-sampleRVarTWith liftN (RVarT m) = runPromptT return bindP bindN m
-    where
-        bindP :: forall t. (Prim t -> (t -> m a) -> m a)
-        bindP prim cont = getRandomPrim prim >>= cont
+{-# INLINE sampleReaderRVarTWith #-}
+sampleReaderRVarTWith ::
+       forall m n a g. (StatefulGen g m, MonadReader g m)
+    => (forall t. n t -> m t)
+    -> RVarT n a
+    -> m a
+sampleReaderRVarTWith liftN (RVarT m) = runPromptT' (\p -> ask >>= uniformPrimM p) liftN m
 
-        bindN :: forall t. n t -> (t -> m a) -> m a
-        bindN nExp cont = liftN nExp >>= cont
+-- |@sampleRVarTWith lift x@ is equivalent to @runRVarTWith lift x 'StdRandom'@.
+{-# INLINE sampleStateRVarTWith #-}
+sampleStateRVarTWith ::
+       forall m n a g. (RandomGen g, MonadState g m)
+    => (forall t. n t -> m t)
+    -> RVarT n a
+    -> m a
+sampleStateRVarTWith liftN (RVarT m) = runPromptT' (\p -> uniformPrimM p StateGenM) liftN m
 
 instance Functor (RVarT n) where
     fmap = liftM
@@ -197,28 +235,31 @@ instance Monad (RVarT n) where
     return x = RVarT (return $! x)
     (RVarT m) >>= k = RVarT (m >>= \x -> x `seq` unRVarT (k x))
 
-instance MonadRandom (RVarT n) where
-    getRandomPrim = RVarT . prompt
-
 instance Applicative (RVarT n) where
     pure  = return
     (<*>) = ap
 
-instance MonadPrompt Prim (RVarT n) where
-    prompt = RVarT . prompt
+instance Functor n => MonadFree Prim (RVarT n) where
+  wrap = RVarT . wrap . InL . fmap unRVarT
 
 instance T.MonadTrans RVarT where
-    lift m = RVarT (MTL.lift m)
+  lift = RVarT . wrap . InR . fmap unRVarT . fmap return
 
 instance T.MonadIO m => T.MonadIO (RVarT m) where
     liftIO = T.lift . T.liftIO
 
-#ifndef MTL2
+data RGen = RGen
 
-instance MTL.MonadTrans RVarT where
-    lift m = RVarT (MTL.lift m)
+instance Functor m => StatefulGen RGen (RVarT m) where
+    uniformWord32 RGen             = liftF (Prim id)
+    {-# INLINE uniformWord32 #-}
+    uniformShortByteString _n RGen = RVarT $ error "ShortByteString"
+    {-# INLINE uniformShortByteString #-}
 
-instance MTL.MonadIO m => MTL.MonadIO (RVarT m) where
-    liftIO = MTL.lift . MTL.liftIO
+uniformRVarT :: (Functor m, Uniform a) => RVarT m a
+uniformRVarT = uniformM RGen
+{-# INLINE uniformRVarT #-}
 
-#endif
+uniformRangeRVarT :: (Functor m, UniformRange a) => (a, a) -> RVarT m a
+uniformRangeRVarT r = uniformRM r RGen
+{-# INLINE uniformRangeRVarT #-}
diff --git a/rvar/src/Data/RVar/Prim.hs b/rvar/src/Data/RVar/Prim.hs
new file mode 100644
index 0000000..f5a80c7
--- /dev/null
+++ b/rvar/src/Data/RVar/Prim.hs
@@ -0,0 +1,18 @@
+{-# LANGUAGE DeriveDataTypeable #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE RankNTypes #-}
+-- |This is an internal interface to support the 'RVar' abstraction.  It
+-- reifies the operations provided by `System.Random.Stateful.StatefulGen` in a
+-- uniform and efficient way, as functions of type @Prim a -> m a@.
+module Data.RVar.Prim (Prim(..)) where
+
+import Data.Word
+
+-- |A type describing a request for a primitive random variate. This
+-- data type is needed for creating
+-- `System.Random.Stateful.StatefulGen` instance for `Data.RVar.RVarT`
+
+data Prim a = Prim (Word32 -> a)
+
+instance Functor Prim where
+  fmap f (Prim k) = Prim (f . k)
diff --git a/stack.yaml b/stack.yaml
index d19a64c..8b446ce 100644
--- a/stack.yaml
+++ b/stack.yaml
@@ -35,7 +35,7 @@ resolver: lts-16.31
 #  - wai
 packages:
 - random-fu
-- random-source
+# - random-source
 - rvar
 - tests/speed
 # Dependency packages to be pulled from upstream that are not in the resolver
@@ -46,7 +46,6 @@ extra-deps:
 - flexible-defaults-0.0.3@sha256:6a7ab000561e1075003cb1053dfbbb4020ae2b02916776d1479c9c3fc82f5d0d,2508
 - splitmix-0.1.0.3@sha256:fc3aae74c467f4b608050bef53aec17904a618731df9407e655d8f3bf8c32d5c,6049
 - mwc-random-0.15.0.1@sha256:841c86f132c45cb81116e1d3a8a150cecc27820c2b4e38f8cf86e3fe7735c2ab,3370
-
 # Override default flag values for local packages and extra-deps
 # flags: {}
 
diff --git a/tests/speed/Bench.hs b/tests/speed/Bench.hs
index 1ec329e..2465278 100644
--- a/tests/speed/Bench.hs
+++ b/tests/speed/Bench.hs
@@ -1,4 +1,5 @@
-{-# LANGUAGE FlexibleContexts, RankNTypes #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE RankNTypes #-}
 
 {-# OPTIONS_GHC -fno-warn-missing-signatures #-}
 {-# OPTIONS_GHC -fno-warn-type-defaults      #-}
@@ -6,7 +7,6 @@
 module Main where
 
 import Data.Random
-import Data.Random.Source.DevRandom
 import Data.Random.Distribution.Beta
 import Data.Random.Distribution.Binomial
 import Data.Random.Distribution.Dirichlet
@@ -16,35 +16,32 @@ import Data.Random.Distribution.Poisson
 import Data.Random.Distribution.Rayleigh
 import Data.Random.Distribution.Triangular
 
-import Control.Applicative ((<$>))
+import System.Random.Stateful (runStateGen_)
+
 import Control.DeepSeq (NFData)
 import Control.Monad
+import qualified Control.Monad.Random as CMR
 import Control.Monad.ST
 import Control.Monad.State
-import qualified Control.Monad.Random as CMR
-import Control.Monad.Trans (lift)
 import Foreign
 import System.Random hiding (uniform)
 import qualified System.Random.MWC as MWC
 import TestSupport
-import qualified Data.Vector.Unboxed as Vec
 
 import Criterion.Main
 
 main = do
     let count = 64000
 
-    mwcSrc <- newGenIO
-    mtSrc  <- newMTSrc
-    stdSrc <- newStdSrc
+    mwcGen <- newMWCGenIO
+    mtGen  <- newMTGenM
+    stdGen <- newStdGenM
 
     defaultMain
         [ bgroup "dists"
-            [ bgroup "MWC"          (dists mwcSrc       count)
-            , bgroup "PureMT"       (dists mtSrc        count)
-            , bgroup "StdGen"       (dists stdSrc       count)
-            , bgroup "DevRandom"    (dists DevRandom    count)
-            , bgroup "DevURandom"   (dists DevURandom   count)
+            [ bgroup "MWC"          (dists mwcGen       count)
+            , bgroup "PureMT"       (dists mtGen        count)
+            , bgroup "StdGen"       (dists stdGen       count)
             ]
 
         , bgroup "IO StdGen"
@@ -52,39 +49,39 @@ main = do
                 sum' <$> replicateM count (randomRIO (10,50)) :: IO Double
 
             , bench "uniform A" $ nfIO $ do
-                sum' <$> replicateM count (sampleFrom stdSrc (uniform 10 50)) :: IO Double
+                sum' <$> replicateM count (sampleFrom stdGen (uniform 10 50)) :: IO Double
 
             , bench "uniform B" $ nfIO $ do
-                sum' <$> sampleFrom stdSrc (replicateM count (uniform 10 50)) :: IO Double
+                sum' <$> sampleFrom stdGen (replicateM count (uniform 10 50)) :: IO Double
             ]
 
         , bgroup "pure StdGen"
             [ bgroup "Double"
                 [ bench "getRandomRs" $ nfIO $ do
-                    src <- newStdGen
-                    let xs = CMR.evalRand (CMR.getRandomRs (10,50)) src
+                    gen <- newStdGen
+                    let xs = CMR.evalRand (CMR.getRandomRs (10,50)) gen
                     return $ sum' (take count xs :: [Double])
                 , bench "RVarT, State - sample replicateM" $ nfIO $ do
-                    src <- newStdGen
-                    let xs = evalState (runRVar (replicateM count (uniform 10 50)) StdRandom) src
+                    gen <- newStdGen
+                    let xs = runStateGen_ gen $ runRVar (replicateM count (uniform 10 50))
                     return (sum' xs :: Double)
                 , bench "RVarT, State - replicateM sample" $ nfIO $ do
-                    src <- newStdGen
-                    let xs = evalState (replicateM count (runRVar (uniform 10 50) StdRandom)) src
+                    gen <- newStdGen
+                    let xs = runStateGen_ gen $ replicateM count . runRVar (uniform 10 50)
                     return (sum' xs :: Double)
                 ]
             , bgroup "Int"
                 [ bench "getRandomRs" $ nfIO $ do
-                    src <- newStdGen
-                    let xs = CMR.evalRand (CMR.getRandomRs (10,50)) src
+                    gen <- newStdGen
+                    let xs = CMR.evalRand (CMR.getRandomRs (10,50)) gen
                     return $ sum' (take count xs :: [Int])
                 , bench "RVarT, State - sample replicateM" $ nfIO $ do
-                    src <- newStdGen
-                    let xs = evalState (runRVar (replicateM count (uniform 10 50)) StdRandom) src
+                    gen <- newStdGen
+                    let xs = runStateGen_ gen $ runRVar (replicateM count (uniform 10 50))
                     return $ (sum' xs :: Int)
                 , bench "RVarT, State - replicateM sample" $ nfIO $ do
-                    src <- newStdGen
-                    let xs = evalState (replicateM count (runRVar (uniform 10 50) StdRandom)) src
+                    gen <- newStdGen
+                    let xs = runStateGen_ gen $ replicateM count . runRVar (uniform 10 50)
                     return $ (sum' xs :: Int)
                 ]
             ]
@@ -92,31 +89,32 @@ main = do
         , bgroup "MWC"
             [ bgroup "stdUniform"
                 [ bench "Double" $ nfIO $ do
-                    src <- newGenIO
-                    xs <- stToIO $ replicateM count (MWC.uniform src)
+                    gen <- newMWCGenIO
+                    xs <- stToIO $ replicateM count (MWC.uniform gen)
                     return $ sum' (xs :: [Double])
                 , bench "Int" $ nfIO $ do
-                    src <- newGenIO
-                    xs <- stToIO $ replicateM count (MWC.uniform src)
+                    gen <- newMWCGenIO
+                    xs <- stToIO $ replicateM count (MWC.uniform gen)
                     return $ sum' (xs :: [Int])
                 ]
             , bgroup "uniform"
                 [ bench "Double" $ nfIO $ do
-                    src <- newGenIO
-                    us <- stToIO $ replicateM count (MWC.uniform src)
+                    gen <- newMWCGenIO
+                    us <- stToIO $ replicateM count (MWC.uniform gen)
                     let xs = [(u - 0.5) * 20| u <- us]
                     return $ sum' (xs :: [Double])
                 ]
 --            , bgroup "normal"
 --                [ bench "Double" $ nfIO $ do
---                    src <- newGenIO
---                    xs <- stToIO $ replicateM count (MWC.normal src)
+--                    gen <- newGenIO
+--                    xs <- stToIO $ replicateM count (MWC.normal gen)
 --                    return $ sum' (xs :: [Double])
 --                ]
+            -- FIXME: uniformVector no longer works on Double
             -- , bgroup "uniformVector"
             --     [ bench "Double" $ nfIO $ do
-            --         src <- newGenIO
-            --         xs <- stToIO $ MWC.uniformVector src count
+            --         gen <- newGenIO
+            --         xs <- stToIO $ MWC.uniformVector gen count
             --         -- unboxed, so don't need to force it, but we sum it
             --         -- anyway to make the comparison fair between other tests
             --         return $ (Vec.sum xs :: Double)
@@ -127,71 +125,71 @@ main = do
 
         ]
 
-dists src count =
-    [ multiSuite src count "uniform"     (Uniform 10 50)
-    , multiSuite src count "stdUniform"  StdUniform
-    , multiSuite src count "poisson"     (Poisson 3 :: Num t => Poisson Double t)
-    , multiSuite src count "binomial 10" (Binomial 10 (0.5 :: Float))
-    , doubleSuite src count "stdNormal"   StdNormal
-    , doubleSuite src count "exponential" (Exp  2)
-    , doubleSuite src count "beta"        (Beta  2 5)
-    , doubleSuite src count "gamma"       (Gamma 2 5)
-    , doubleSuite src count "triangular"  (Triangular 2 5 14)
-    , doubleSuite src count "rayleigh"    (Rayleigh 1.6)
+dists gen count =
+    [ multiSuite gen count "uniform"     (Uniform 10 50)
+    , multiSuite gen count "stdUniform"  StdUniform
+    , multiSuite gen count "poisson"     (Poisson 3 :: Num t => Poisson Double t)
+    , multiSuite gen count "binomial 10" (Binomial 10 (0.5 :: Float))
+    , doubleSuite gen count "stdNormal"   StdNormal
+    , doubleSuite gen count "exponential" (Exp  2)
+    , doubleSuite gen count "beta"        (Beta  2 5)
+    , doubleSuite gen count "gamma"       (Gamma 2 5)
+    , doubleSuite gen count "triangular"  (Triangular 2 5 14)
+    , doubleSuite gen count "rayleigh"    (Rayleigh 1.6)
 
     , bench "dirichlet" $ nfIO $ do
-        sum' <$> sampleFrom src (dirichlet [1..fromIntegral count :: Double])
+        sum' <$> sampleFrom gen (dirichlet [1..fromIntegral count :: Double])
 
     , bgroup "multinomial"
         [ bgroup "many p"
             [ bench desc $ nfIO $ do
-                sum' <$> sampleFrom src (multinomial [1..1e4 :: Double] (n :: Int))
+                sum' <$> sampleFrom gen (multinomial [1..1e4 :: Double] (n :: Int))
             | (desc, n) <- [("small n", 10), ("medium n", 10^4), ("large n", 10^8)]
             ]
         , bgroup "few p"
             [bench desc $ nfIO $ do
                 replicateM_ 1000 $ do
-                    sum' <$> sampleFrom src (multinomial [1..10 :: Double] (n :: Int))
+                    sum' <$> sampleFrom gen (multinomial [1..10 :: Double] (n :: Int))
             | (desc, n) <- [("small n", 10), ("medium n", 10^4), ("large n", 10^8)]
             ]
         ]
 
     , bench "shuffle" $ nfIO $ do
-        sum' <$> sampleFrom src (shuffle [1..count])
+        sum' <$> sampleFrom gen (shuffle [1..count])
     ]
 
-multiSuite :: (Distribution d Double, Distribution d Int, RandomSource IO s) => s -> Int -> String -> (forall t. Num t => d t) -> Benchmark
-multiSuite src count name dist = bgroup name
-    [ doubleSuite src count "Double" dist
-    , intSuite    src count "Int"    dist
+multiSuite :: (Distribution d Double, Distribution d Int, StatefulGen g IO) => g -> Int -> String -> (forall t. Num t => d t) -> Benchmark
+multiSuite gen count name dist = bgroup name
+    [ doubleSuite gen count "Double" dist
+    , intSuite    gen count "Int"    dist
     ]
 
-doubleSuite :: (Distribution d Double, RandomSource IO s) => s -> Int -> String -> d Double -> Benchmark
+doubleSuite :: (Distribution d Double, StatefulGen g IO) => g -> Int -> String -> d Double -> Benchmark
 doubleSuite = suite
 
-intSuite :: (Distribution d Int, RandomSource IO s) => s -> Int -> String -> d Int -> Benchmark
+intSuite :: (Distribution d Int, StatefulGen g IO) => g -> Int -> String -> d Int -> Benchmark
 intSuite = suite
 
-suite :: (Storable t, Num t, Distribution d t, NFData t, RandomSource IO s) => s -> Int -> String -> d t -> Benchmark
-suite src count name var = bgroup name
+suite :: (Storable t, Num t, Distribution d t, NFData t, StatefulGen g IO) => g -> Int -> String -> d t -> Benchmark
+suite gen count name var = bgroup name
     [ bench "single sample" $ nfIO $ do
-        sampleFrom src var
+        sampleFrom gen var
 
     -- Ideally, these would all be the same speed:
     , bench "sum of samples (implicit rvar)" $ nfIO $ do
-        sumM count (sampleFrom src var)
+        sumM count (sampleFrom gen var)
 
     , bench "sum of samples (explicit rvar)" $ nfIO $ do
-        sumM count (sampleFrom src (rvar var))
+        sumM count (sampleFrom gen (rvar var))
 
     , bench "sample of sum" $ nfIO $ do
-        sampleFrom src (sumM count (rvar var))
+        sampleFrom gen (sumM count (rvar var))
 
     , bench "array of samples" $ nfIO $ do
         allocaArray count $ \ptr -> do
             sequence_
                 [ do
-                    x <- sampleFrom src var
+                    x <- sampleFrom gen var
 
                     pokeElemOff ptr offset x
                 | offset <- [0 .. count-1]
@@ -199,7 +197,7 @@ suite src count name var = bgroup name
             sumBuf count ptr
 
     , bench "RVarT IO arrays" $ nfIO $ do
-        allocaArray count $ \ptr -> flip runRVarT src $ do
+        allocaArray count $ \ptr -> flip runRVarT gen $ do
             sequence_
                 [ do
                     x <- rvarT var
diff --git a/tests/speed/TestSupport.hs b/tests/speed/TestSupport.hs
index 71aae7f..20d5e68 100644
--- a/tests/speed/TestSupport.hs
+++ b/tests/speed/TestSupport.hs
@@ -2,41 +2,25 @@
 
 {-# OPTIONS_GHC -fno-warn-simplifiable-class-constraints #-}
 {-# OPTIONS_GHC -fno-warn-type-defaults                  #-}
-{-# OPTIONS_GHC -fno-warn-missing-signatures              #-}
+{-# OPTIONS_GHC -fno-warn-missing-signatures             #-}
 
 module TestSupport where
 
 import System.Random.Mersenne.Pure64
 import System.Random.MWC
 import Data.List
-import Data.StateRef
 import Foreign
-import System.Random
+import System.Random.Stateful
+import Control.Monad.ST (RealWorld)
 
--- type Src = IORef PureMT
--- getTestSource = newMTSrc
+newMTGenM :: IO (IOGenM PureMT)
+newMTGenM = newIOGenM =<< newPureMT
 
--- type Src = IORef StdGen
--- getTestSource = newStdSrc
+newStdGenM :: IO (IOGenM StdGen)
+newStdGenM = newIOGenM =<< newStdGen
 
-type Src = Gen RealWorld
-getTestSource :: IO Src
-getTestSource = newGenIO
-
-newMTSrc :: IO (IORef PureMT)
-newMTSrc = do
-    mt <- newPureMT
-    newReference mt
-
-newStdSrc :: IO (IORef StdGen)
-newStdSrc = do
-    mt <- newStdGen
-    newReference mt
-
-newGenIO :: IO (Gen RealWorld)
-newGenIO = do
-    seed <- withSystemRandom (save :: Gen RealWorld -> IO Seed)
-    restore seed
+newMWCGenIO :: IO (Gen RealWorld)
+newMWCGenIO = createSystemRandom
 
 
 sum' xs = foldl' (+) 0 xs
diff --git a/tests/speed/speed-tests.cabal b/tests/speed/speed-tests.cabal
index 1d1ea8b..93e7e24 100644
--- a/tests/speed/speed-tests.cabal
+++ b/tests/speed/speed-tests.cabal
@@ -13,15 +13,12 @@ description:            Various benchmarks for the random-fu library.
 flag split-random-fu
 
 benchmark random-fu-bench
-  type:             exitcode-stdio-1.0
+  type:                 exitcode-stdio-1.0
   main-is:              Bench.hs
   build-depends:        base >= 4, criterion, MonadRandom, mtl,
-                        stateref, mersenne-random-pure64, mwc-random >= 0.5,
-                        random, vector, deepseq
-  if flag(split-random-fu)
-    build-depends:      random-fu, random-source
-  else
-    build-depends:      random-fu
+                        stateref, mersenne-random-pure64, mwc-random >= 0.15.0.1,
+                        random, vector, deepseq,
+                        random-fu
  other-modules:         TestSupport