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.2.0
+version:                0.0.2.1
 stability:              experimental
 
 cabal-version:          >= 1.2
@@ -16,9 +16,7 @@
                         for entropy sources and random variable distributions, all
                         served up on a monadic platter.  Aspires to be useful
                         in an idiomatic way in both \"pure\" and \"impure\" styles,
-                        as well as reasonably fast.  May not yet meet the latter
-                        goal, but I think the former is starting to shape up
-                        nicely.
+                        as well as reasonably fast.
 
 Flag base4
 
@@ -51,17 +49,19 @@
                         Data.Random.Source.PureMT
                         Data.Random.Source.Std
   if flag(base4)
-    build-depends:      base >= 4,
+    build-depends:      base >= 4 && <5,
                         syb
   else
     build-depends:      base >= 3 && < 4
     
   build-depends:        array,
                         containers,
+                        erf,
                         mersenne-random-pure64,
                         monad-loops >= 0.3.0.1,
                         mtl,
                         random,
+                        random-shuffle,
                         stateref,
                         storablevector,
                         template-haskell
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
@@ -9,9 +9,6 @@
 
 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 using 
@@ -23,7 +20,7 @@
 
 class Distribution d t => CDF d t where
     -- |Return the cumulative distribution function of this distribution.
-    -- That is, a function taking 'x :: t' to the probability that the next
+    -- That is, a function taking @x :: t@ to the probability that the next
     -- sample will return a value less than or equal to x, according to some
     -- order or partial order (not necessarily an obvious one).
     --
@@ -31,7 +28,7 @@
     -- to the CDF with respect to that order.
     -- 
     -- In other cases, 'cdf' is only required to satisfy the following law:
-    -- > fmap (cdf d) (rvar d)
+    -- @fmap (cdf d) (rvar d)@
     -- must be uniformly distributed over (0,1).  Inclusion of either endpoint is optional,
     -- though the preferred range is (0,1].
     -- 
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
@@ -12,10 +12,8 @@
 
 import Data.Random.Internal.TH
 
-import Data.Random.Source
-import Data.Random.Distribution
 import Data.Random.RVar
-
+import Data.Random.Distribution
 import Data.Random.Distribution.Uniform
 
 import Data.Int
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
@@ -9,7 +9,6 @@
 
 module Data.Random.Distribution.Beta where
 
-import Data.Random.Source
 import Data.Random.RVar
 import Data.Random.Distribution
 import Data.Random.Distribution.Gamma
@@ -29,12 +28,6 @@
     x <- erlang a
     y <- erlang b
     return (x / (x + y))
-
-fractionalBetaMinusHalf :: (Fractional a, Distribution Gamma a) => a -> a -> RVar a
-fractionalBetaMinusHalf a b = do
-    x <- gamma a 1
-    y <- gamma b 1
-    return ((x-y)/(2*(x+y)))
 
 beta :: Distribution Beta a => a -> a -> RVar a
 beta a b = rvar (Beta 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
@@ -11,11 +11,8 @@
 
 import Data.Random.Internal.TH
 
-import Data.Random.Source
-import Data.Random.Distribution
 import Data.Random.RVar
-
-import Data.Random.Distribution.Bernoulli
+import Data.Random.Distribution
 import Data.Random.Distribution.Beta
 import Data.Random.Distribution.Uniform
 
@@ -43,15 +40,14 @@
         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
-                a <- selectA (fromInteger (floatRadix undefined) ^ floatDigits undefined) t
-                    
-                let b = 1 + t - a
-                
+                let a = 1 + t `div` 2
+                    b = 1 + t - a
+        
                 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))
-            
+        
             | otherwise = count k t
                 where
                     count k  0    = return k
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
@@ -9,7 +9,6 @@
 
 module Data.Random.Distribution.Exponential where
 
-import Data.Random.Source
 import Data.Random.RVar
 import Data.Random.Distribution
 import Data.Random.Distribution.Uniform
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
@@ -16,7 +16,6 @@
 
 module Data.Random.Distribution.Gamma where
 
-import Data.Random.Source
 import Data.Random.RVar
 import Data.Random.Distribution
 import Data.Random.Distribution.Uniform
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
@@ -34,11 +34,7 @@
 import Control.Monad
 import Foreign.Storable
 
-foreign import ccall "math.h erf" erf :: Double -> Double
-foreign import ccall "math.h erfc" erfc :: Double -> Double
-foreign import ccall "math.h erff" erff :: Float -> Float
-erfg :: RealFrac a => a -> a
-erfg = realToFrac . erf . realToFrac
+import Data.Number.Erf
 
 normalPair :: (Floating a, Distribution StdUniform a) => RVar (a,a)
 normalPair = boxMullerNormalPair
@@ -86,11 +82,11 @@
                 else return (r - x)
 
 -- |Construct a 'Ziggurat' for sampling a normal distribution, given
--- a suitable error function, logBase 2 c, and the 'zGetIU' implementation.
+-- logBase 2 c, and the 'zGetIU' implementation.
 normalZ ::
-  (RealFloat a, Storable a, Distribution Uniform a, Integral b) =>
-  (a -> a) -> b -> RVar (Int, a) -> Ziggurat a
-normalZ erf p = mkZigguratRec True normalF normalFInv (normalFInt erf) normalFVol (2^p)
+  (RealFloat a, Erf a, Storable a, Distribution Uniform a, Integral b) =>
+  b -> RVar (Int, a) -> Ziggurat a
+normalZ p = mkZigguratRec True normalF normalFInv normalFInt normalFVol (2^p)
 
 -- | Ziggurat target function
 normalF :: (Floating a, Ord a) => a -> a
@@ -100,17 +96,17 @@
 -- | inverse of 'normalF'
 normalFInv :: Floating a => a -> a
 normalFInv y  = sqrt ((-2) * log y)
--- | integral of 'normalF', parameterized over 'erf'
-normalFInt :: (Floating a, Ord a) => (a -> a) -> a -> a
-normalFInt erf x 
+-- | integral of 'normalF'
+normalFInt :: (Floating a, Erf a, Ord a) => a -> a
+normalFInt x 
     | x <= 0    = 0
     | otherwise = normalFVol * erf (x * sqrt 0.5)
 -- | volume of 'normalF'
 normalFVol :: Floating a => a
 normalFVol = sqrt (0.5 * pi)
 
-realFloatStdNormal :: (RealFloat a, Storable a, Distribution Uniform a) => RVar a
-realFloatStdNormal = rvar (normalZ erfg p getIU)
+realFloatStdNormal :: (RealFloat a, Erf a, Storable a, Distribution Uniform a) => RVar a
+realFloatStdNormal = runZiggurat (normalZ p getIU)
     where 
         p = 6
         
@@ -120,7 +116,7 @@
             return (fromIntegral i .&. (2^p-1), u)
 
 doubleStdNormal :: RVar Double
-doubleStdNormal = rvar doubleStdNormalZ
+doubleStdNormal = runZiggurat doubleStdNormalZ
 
 -- doubleStdNormalC must not be over 12 if using wordToDoubleWithExcess
 doubleStdNormalC :: Int
@@ -142,7 +138,7 @@
             return (fromIntegral i .&. (doubleStdNormalC-1), u+u-1)
 
 floatStdNormal :: RVar Float
-floatStdNormal = rvar floatStdNormalZ
+floatStdNormal = runZiggurat floatStdNormalZ
 
 floatStdNormalC :: Int
 floatStdNormalC = 512
@@ -168,29 +164,33 @@
 normalPdf :: Real a => a -> a -> a -> Double
 normalPdf m s x = recip (realToFrac s * sqrt (2*pi)) * exp (-0.5 * (realToFrac x - realToFrac m)^2 / (realToFrac s)^2)
 
-normalCdf :: Real a => a -> a -> a -> Double
-normalCdf m s x = 0.5 * (1 + erf ((realToFrac x - realToFrac m) / (realToFrac s * sqrt 2)))
+normalCdf :: (Real a, Erf a) => a -> a -> a -> Double
+normalCdf m s x = realToFrac (normcdf ((x-m) / s))
 
 data Normal a
     = StdNormal
     | Normal a a -- mean, sd
 
 instance Distribution Normal Double where
+    {-# SPECIALIZE instance Distribution Normal Double #-}
     rvar StdNormal = doubleStdNormal
     rvar (Normal m s) = do
         x <- doubleStdNormal
         return (x * s + m)
 
 instance Distribution Normal Float where
+    {-# SPECIALIZE instance Distribution Normal Float #-}
     rvar StdNormal = floatStdNormal
     rvar (Normal m s) = do
         x <- floatStdNormal
         return (x * s + m)
 
-instance (Real a, Distribution Normal a) => CDF Normal a where
+instance (Real a, Erf a, Distribution Normal a) => CDF Normal a where
     cdf StdNormal    = normalCdf 0 1
     cdf (Normal m s) = normalCdf m s
 
+{-# SPECIALIZE stdNormal :: RVar Double #-}
+{-# SPECIALIZE stdNormal :: RVar Float #-}
 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
@@ -11,17 +11,14 @@
 
 import Data.Random.Internal.TH
 
-import Data.Random.Source
-import Data.Random.Distribution
 import Data.Random.RVar
-
+import Data.Random.Distribution
 import Data.Random.Distribution.Uniform
 import Data.Random.Distribution.Gamma
 import Data.Random.Distribution.Binomial
 
 import Data.Int
 import Data.Word
-
 import Control.Monad
 
 -- from Knuth, with interpretation help from gsl sources
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
@@ -37,10 +37,8 @@
 import Data.Random.Distribution
 import Data.Random.RVar
 
-import Data.Ratio
 import Data.Word
 import Data.Int
-import Data.Bits
 import Data.List
 
 import Control.Monad.Loops
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
@@ -10,7 +10,6 @@
 
 import Data.Bits
 import Data.Word
-import Control.Monad
 
 {-# INLINE buildWord #-}
 buildWord :: Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word64
diff --git a/src/Data/Random/List.hs b/src/Data/Random/List.hs
--- a/src/Data/Random/List.hs
+++ b/src/Data/Random/List.hs
@@ -1,19 +1,10 @@
-{-
- -      ``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
+import qualified System.Random.Shuffle as SRS
+import Control.Monad
 
 randomElement :: [a] -> RVar a
 randomElement [] = error "randomElement: empty list!"
@@ -21,47 +12,16 @@
     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 [] = return []
+shuffle xs = do
+    is <- zipWithM (\_ i -> uniform 0 i) (tail xs) [1..]
+    
+    return (SRS.shuffle xs (reverse is))
 
--- |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)
+shuffleN :: Int -> [a] -> RVar [a]
+shuffleN 0 xs = return []
+shuffleN (n+1) xs = do
+    is <- sequence [uniform 0 i | i <- [n,n-1..1]]
+    return (SRS.shuffle xs is)
+    
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
@@ -23,16 +23,13 @@
     ) where
 
 
-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
 
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
@@ -11,8 +11,6 @@
     ) where
 
 import Data.Word
-import Data.Bits
-import Data.List
 import Control.Monad
 
 import Data.Random.Internal.Words
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
@@ -11,7 +11,6 @@
 
 import Data.Random.Source
 
-import GHC.IOBase (unsafePerformIO)
 import System.IO (openBinaryFile, hGetBuf, IOMode(..))
 import Foreign
 
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
@@ -10,7 +10,6 @@
 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