diff --git a/NEWS b/NEWS
--- a/NEWS
+++ b/NEWS
@@ -1,3 +1,35 @@
+Changes in 0.6:
+
+* Major overhaul. A lot of client code will break.
+
+* Replaced the old MC type with a monad transformer,
+  requiring the base monad to be an instance of PrimMonad.
+
+* Removed the MCT type.
+
+* Removed the MonadMC class.
+
+* Removed unsafe operations relying on lazy IO (unsafeInterleaveIO),
+  or replaced them by safe, strict versions.
+
+* Added new fold operation foldMC.
+
+* Refactored and simplified sampling and shuffling functions.  The
+  "sampleSubset" functions are strict now.
+
+* Hide Data.Summary.Utils, and remove Data.Summary.
+
+* Change Data.Summary.Bool and Data.Summary.Double interfaces to mimic
+  the Data.Set functions.  Like Data.Set, these modules should now
+  be used with qualified imports.
+
+* Removed NFData instances from Summary types; add Eq instances.
+
+* Add Eq, Show, Data, Typeable instances.
+
+* Updated examples.
+
+
 Changes in 0.5:
 
 * Clark Gaebel added Monoid and NFData instances to Summary types
diff --git a/examples/Binomial.hs b/examples/Binomial.hs
--- a/examples/Binomial.hs
+++ b/examples/Binomial.hs
@@ -2,15 +2,15 @@
     where
 
 import Control.Monad
+import Control.Monad.Primitive( PrimMonad )
 import Data.List( foldl' )
 import Text.Printf( printf )
 
 import Control.Monad.MC
-import Data.Summary
-import Data.Summary.Utils( inInterval )
+import qualified Data.Summary.Double as S
 
 -- | Sample from a binomial distribution with the given parameters.
-binomial :: (MonadMC m) => Int -> Double -> m Int
+binomial :: (PrimMonad m) => Int -> Double -> MC m Int
 binomial n p = let
     q     = 1 - p
     probs = map (\i -> (fromIntegral $ n `choose` i) * p^^i * q^^(n-i)) [0..n]
@@ -18,29 +18,31 @@
 
 -- | Get a sample confidence interval for the mean after @reps@ replications of
 -- a binomial with the given parameters.
-binomialMean :: (MonadMC m) => Int -> Double -> Int -> m (Double,Double)
+binomialMean :: (PrimMonad m) => Int -> Double -> Int -> MC m (Double,Double)
 binomialMean n p reps =
-    liftM (sampleCI 0.95 . summary . map fromIntegral) $
-        replicateMC reps (binomial n p)
+    liftM (S.meanCI 0.95) $
+        foldMC (\s x -> return $! S.insertWith fromIntegral x s) S.empty
+               reps (binomial n p)
 
 -- | Compute @reps@ 95% confidence intervals for the mean of an @(n,p)@
 -- binormal based on samples of the given size, and record the number
 -- of intervals that contain the true mean.
-coverage :: (MonadMC m) => Int -> Double -> Int -> Int -> m Int
+coverage :: (PrimMonad m) => Int -> Double -> Int -> Int -> MC m Int
 coverage n p size reps =
-    liftM (length . filter (mu `inInterval`)) $
-        replicateMC reps $
-            binomialMean n p size
+    foldMC (\tot ci -> return $! update tot (mu `inInterval` ci)) 0
+           reps (binomialMean n p size)
   where
     mu = fromIntegral n * p
+    x `inInterval` (a,b) = x > a && x < b
+    update tot b = tot + (if b then 1 else 0)
 
 main =
     let seed = 0
-        reps = 100
+        reps = 10000
         n    = 10
         p    = 0.2
         size = 500
-        c    =  coverage n p size reps `evalMC` mt19937 seed 
+        c    = (coverage n p size reps) `evalMC` (mt19937 seed)
     in
         printf "\nOf %d 95%%-intervals, %d contain the true value.\n" reps c
 
diff --git a/examples/Pi.lhs b/examples/Pi.lhs
--- a/examples/Pi.lhs
+++ b/examples/Pi.lhs
@@ -4,11 +4,11 @@
 of them fall in the unit circle.
 
 \begin{code}
-import Control.Monad( liftM, liftM2 )
-import Control.Monad.MC( MC, uniform, replicateMC, evalMC, mt19937 )
-import Data.Summary.Bool( summary, sampleMean, sampleSE )
-import Data.Summary.Utils( interval )
-import Text.Printf
+import Control.Monad( liftM2 )
+import Control.Monad.MC( STMC, evalMC, foldMC, mt19937, uniform )
+import Data.Monoid( (<>), mempty )
+import qualified Data.Summary.Bool as S
+import Text.Printf( printf )
 \end{code}
 
 First, we need a function to test whether or not a point is in the
@@ -19,75 +19,26 @@
 inUnitCircle (x,y) = x*x + y*y <= 1
 \end{code}
 
-The first line is the type signature, which tells us that "inUnitCircle"
-is a function which takes a pair of `Double`s and returns a `Bool`.  In
-English, "::" means "has type".  The second line is the one that defines
-the function.
-
-\begin{code}
-estimatePi :: [(Double,Double)] -> (Double,Double)
-estimatePi xs =
-    let s       = summary $ map inUnitCircle xs
-        (mu,se) = (sampleMean s, sampleSE s) in
-    (4*mu,4*se)
-\end{code}
-
 Next, we need to generate a random point in the unit box.  In the
 Control.Monad.MC module, there is a function for generating uniform
-values in an interval, called "uniform".  This function has a 
+values in an interval, called "uniform".  This function has a
 funny-looking type, but you can think of it as:
 
-    uniform :: Double -> Double -> MC Double
+    uniform :: Double -> Double -> STMC s Double
 
-This type means that the function takes the two endpoints of the 
+This type means that the function takes the two endpoints of the
 interval as arguments, and returns a Monte-Carlo action which produces
 a Double.
 
-You can think of the type "MC Double" as a random number generator.
-For general types, "MC a" is a generator for values of type "a".  In
-fact, "MC" is one of a general class of objects called a Monad.  We
-use Monads in Haskell for making sure events happen in the right order.
-That is, if we have three parts of a simulation, say A, B, and C, and
-we want them to happen in the order
- 
-    A ====> B ====> C
-
-then we would like to make sure that A is done consuming random numbers
-before B consumes anything.  Likewise, we want B to finish before C
-starts.  Monads are the magic that enable us to ensure this.
-
-There are a number of functions in the standard library for working with
-monads.  The first we will use is
-
-    liftM2 :: (Monad m) => (a1 -> a2 -> r) -> m a1 -> m a2 -> m r
-
-This function works on *any* Monad.  When we use it on the MC monad, it
-will have type
-
-    liftM2 :: (a1 -> a2 -> r) -> MC a1 -> MC a2 -> MC r
-    
-What liftM2 does is it takes a function of two arguments and two Monte
-Carlo actions.  It returns a new Monte Carlo action that does the following:
-
-  1. generate a random value of type a1 using the first action
-  2. generate a random value of type a2 using the second action
-  3. apply a function to the two values and return the result
-  
-We do not need to write liftM2 ourselves, since it is provided in the 
-"Control.Monad" module.  But, if we did have to define it, the code would
-look like:
+You can think of the type "STMC s Double" as a random number generator.
+For general types, "STMC s a" is a generator for values of type "a".
 
-liftM2 f ma1 ma2 = do
-    a1 <- ma1
-    a2 <- ma2
-    return (f a1 a2)
-    
-This code for this uses the "do" notation of Haskell, which allows us
-to specify a series of actions in sequential order.
+The following code generates an x value in the range (-1,1), then a y
+value in the range (-1,1), then returns the pair (x,y):
 
 \begin{code}
-unitBox :: MC (Double,Double)
-unitBox = liftM2 (,) (uniform (-1) 1) 
+unitBox :: STMC s (Double,Double)
+unitBox = liftM2 (,) (uniform (-1) 1)
                      (uniform (-1) 1)
 \end{code}
 
@@ -95,17 +46,20 @@
 -- the sample mean and standard error.
 
 \begin{code}
-simulation :: Int -> MC (Double,Double)
-simulation n = 
-    estimatePi `fmap` replicateMC n unitBox
+estimatePi :: Int -> STMC s (Double,Double)
+estimatePi n = let
+    circ = fmap inUnitCircle unitBox
+    mc   = foldMC (\s x -> return $! S.insert x s) S.empty n circ
+    in fmap (\s -> let (mu,se) = (S.mean s, S.meanSE s)
+                   in (4*mu,4*se)) mc
 \end{code}
 
 \begin{code}
 main =
     let seed    = 0
         n       = 1000000 
-        (mu,se) = simulation n `evalMC` mt19937 seed
-        (l,u)   = interval 0.95 mu se
+        (mu,se) = estimatePi n `evalMC` mt19937 seed
+        (l,u)   = (mu - 2.576 * se, mu + 2.576 * se)
     in do
         printf "\nEstimate: %g" mu
         printf "\n99%% Confidence Interval: (%g, %g)" l u
diff --git a/examples/Poker.hs b/examples/Poker.hs
--- a/examples/Poker.hs
+++ b/examples/Poker.hs
@@ -1,16 +1,18 @@
 module Main where
-    
+
 import Control.Monad
+import Control.Monad.Primitive( PrimMonad ) 
 import Control.Monad.MC
+import Control.Monad.ST
 import Data.List
 import Data.Map( Map )
 import qualified Data.Map as Map
 import Text.Printf
-    
+
 -- | Data types for representing cards.  An Ace has 'number' equal to @1@.
 -- Jack, Queen, and King have numbers @11@, @12@, and @13@, respectively.
 data Suit = Club | Diamond  | Heart | Spade deriving (Eq, Show)
-data Card = Card { number :: Int 
+data Card = Card { number :: Int
                  , suit   :: Suit
                  }
           deriving (Eq, Show)
@@ -24,19 +26,19 @@
 
 -- | Get a list of cards that make up a 52-card deck.
 deck :: [Card]
-deck = [ Card i s 
+deck = [ Card i s
        | i <- [ ace..king ]
        , s <- [ Club, Diamond, Heart, Spade ] ]
 
 -- | A type for the various poker hands.
 data Hand = HighCard  | Pair | TwoPair | ThreeOfAKind | Straight | Flush
-          | FullHouse | FourOfAKind | StraightFlush 
+          | FullHouse | FourOfAKind | StraightFlush
           deriving (Eq, Show, Ord)
 
 -- | Determine the hand corresponding to a list of five cards.
 hand :: [Card] -> Hand
-hand cs = 
-    case matches of 
+hand cs =
+    case matches of
         [1,1,1,1,1] -> case undefined of
                            _ | isStraight && isFlush -> StraightFlush
                            _ | isFlush               -> Flush
@@ -50,7 +52,7 @@
   where
     (x:xs) = (sort . map number) cs
     (s:ss) = map suit cs
-    
+
     isStraight | x == ace && xs == [ 10..king ] = True
                | otherwise                      = xs == [ x+1..x+4 ]
 
@@ -58,9 +60,9 @@
 
     matches = (sort . map length . group) (x:xs)
 
-    
+
 -- | Deal a five-card hand by choosing a random subset of the deck.
-deal :: (MonadMC m) => m [Card]
+deal :: (PrimMonad m) => MC m [Card]
 deal = sampleSubset deck 5
 
 -- | A type for storing the frequencies of the various hands.
@@ -78,15 +80,15 @@
 main =
     let seed   = 0
         reps   = 100000
-        counts = foldl' updateCounts emptyCounts $ 
-                     replicateMC reps deal `evalMC` mt19937 seed 
+        counts = foldl' updateCounts emptyCounts $
+                     replicateMC reps deal (mt19937 seed)
     in do
         printf "\n"
         printf "    Hand       Count    Probability     99%% Interval   \n"
         printf "-------------------------------------------------------\n"
         forM_ ((reverse . Map.toAscList) counts) $ \(h,c) ->
             let n     = fromIntegral reps :: Double
-                p     = fromIntegral c / n 
+                p     = fromIntegral c / n
                 se    = sqrt (p * (1 - p) / n)
                 delta = 2.575829 * se
                 (l,u) = (p-delta, p+delta) in
diff --git a/examples/Queue.hs b/examples/Queue.hs
--- a/examples/Queue.hs
+++ b/examples/Queue.hs
@@ -1,34 +1,44 @@
 
 import Control.Monad
-import Control.Monad.MC
+import Control.Monad.Primitive( PrimMonad )
 import Data.List( foldl' )
-import Data.Summary
 import Text.Printf( printf )
 
+import Control.Monad.MC
+import Data.Summary.Double( Summary )
+import qualified Data.Summary.Double as S
+
+
 -- | There a three items on the menu.
 data Item = Cheeseburger | Fries | Milkshake
 
+
 -- | A customer orders some number of items
 data Customer = Customer { orderOf :: [Item] }
 
+
 -- | The order size is a Poisson random variable with mean 2.
-orderSize :: MC Int
+orderSize :: (PrimMonad m) => MC m Int
 orderSize = liftM (1+) $ poisson 2
 
+
 -- | The items are sampled with the given weights.
-item :: MC Item
+item :: (PrimMonad m) => MC m Item
 item = sampleWithWeights [ (4, Cheeseburger), (2, Fries), (1, Milkshake) ]
 
+
 -- | Generate a random order.
-order :: MC [Item]
+order :: (PrimMonad m) => MC m [Item]
 order = do
     n <- orderSize
     replicateM n item
 
+
 -- | Generate a random customer.
-customer :: MC Customer
+customer :: (PrimMonad m) => MC m Customer
 customer = liftM Customer order
-    
+
+
 -- | A customer event.  The interarrival time is the time that elapeses
 -- between when the previous customer arrives and when the current customer 
 -- arrives.
@@ -36,32 +46,37 @@
                                    , interarrivalTime :: !Double
                                    }
 
+
 -- | Generate a random customer event.  The interarrival time distribution
 -- is exponential with mean 1.
-customerEvent :: MC CustomerEvent
+customerEvent :: (PrimMonad m) => MC m CustomerEvent
 customerEvent = do
     c     <- customer
     delta <- exponential 10
     return $ CustomerEvent c delta
 
+
 -- | The time it takes to make an item.
-cook :: Item -> MC Double
+cook :: (PrimMonad m) => Item -> MC m Double
 cook Cheeseburger = exponential 3
 cook Fries        = exponential 1
 cook Milkshake    = exponential 2
 
+
 -- | The time it takes to cook all of the items in the list is equal
 -- to the maximum time.
-cookAll :: [Item] -> MC Double
+cookAll :: (PrimMonad m) => [Item] -> MC m Double
 cookAll items = do
     ts <- mapM cook items
     return $ foldl' max 0 ts
 
+
 -- | A customer in line, along with how long they have been waiting.
 data Waiting = Waiting { waiting        :: !Customer
                        , hasBeenWaiting :: !Double
                        }
 
+
 -- | A customer, along with how long it takes to prepare the customer's order
 -- and how long the customer has to wait.
 data Service = Service { serving     :: !Customer
@@ -69,37 +84,49 @@
                        , serviceTime :: !Double
                        }
 
+
 -- | Given a customer who has been wating in line, provide them with service.
 -- If the customer has been waiting for longer than 5 minutes, work twice as
 -- fast to cook the food.
-serveWaiting :: Waiting -> MC Service
+serveWaiting :: (PrimMonad m) => Waiting -> MC m Service
 serveWaiting (Waiting c w) = do
     t <- cookAll $ orderOf c
     let t' = if w > 5 then 0.5*t else t
     return $ Service c w t'
 
+
 -- | A resturant has one server, who may be busy. There is a list of
 -- customers wating in line.
 data Restaurant = Restaurant { inProgress  :: Maybe InProgress
                              , waitingLine :: [Waiting]
                              }
 
+
 -- | An in-progress service event.
 data InProgress = InProgress { service      :: !Service
                              , timeToFinish :: !Double
                              }
 
+
+-- | An empty restaurant.
+emptyRestaurant :: Restaurant
+emptyRestaurant = Restaurant Nothing []
+
+
 -- | Update the amount of time the customers have been waiting by adding
 -- the given amount.
 addToWait :: Double -> [Waiting] -> [Waiting]
 addToWait delta = map (\(Waiting w t) -> Waiting w (t+delta))
 
+
 -- | Serve customers in the restaurant for the given amount of time.        
-serveForTime :: Double -> Restaurant -> MC ([Service], Restaurant)
-serveForTime = 
+serveForTime :: (PrimMonad m) => Double
+                              -> Restaurant
+                              -> MC m ([Service], Restaurant)
+serveForTime =
     let serveForTimeHelp ss t r = case r of
             -- When no one is being served and no one is in line, do nothing.
-            Restaurant Nothing  [] -> 
+            Restaurant Nothing  [] ->
                 return $ (ss, r)
 
             -- When no one is being served, take the first person in line
@@ -124,66 +151,74 @@
                                        xs'    = addToWait t xs
                                        r'     = Restaurant y' xs' in
                                    return (ss,r')
-    in serveForTimeHelp []                                   
+    in serveForTimeHelp []
 
+
 -- | Given a new customer arrival event, produce a list of all of the new
 -- service events that happen before the customer gets there, and return
 -- the updated restaurant state at the time immediately after the customer
 -- arrives.
-processEvent :: CustomerEvent    
-             -> Restaurant       
-             -> MC ([Service], Restaurant)
+processEvent :: (PrimMonad m) => CustomerEvent
+                              -> Restaurant
+                              -> MC m ([Service], Restaurant)
 processEvent (CustomerEvent c t) r = do
     (ss,(Restaurant y xs)) <- serveForTime t r
     return $ (ss, (Restaurant y $ xs ++ [Waiting c 0]))
 
+
 -- | Finish serving all of the customers in line.
-finishServing :: Restaurant -> MC [Service]
+finishServing :: (PrimMonad m) => Restaurant -> MC m [Service]
 finishServing r = do
     (ss,_) <- serveForTime infinity r
     return ss
   where
     infinity = 1/0
-    
--- | A restaurant takes a list of customer events and generates a random
--- list of service events.  The reason for the call to "unsafeInterleaveMC"
--- is that we want to make sure that we return a lazy list.   Without it,
--- the function will return only after it has consumed all of the random
--- numbers it needs.  This is problemeatic if the input list is large or
--- or infinite.
-restaurant :: [CustomerEvent] -> MC [Service]
-restaurant = 
-    let restaurantHelp r []     = finishServing r
-        restaurantHelp r (c:cs) = unsafeInterleaveMC $ do
-            (ss,r') <- processEvent c r
-            ss'     <- restaurantHelp r' cs
-            return $ ss ++ ss'
-    in restaurantHelp (Restaurant Nothing [])
 
+
+-- | Run a restaurnt.  Whenever a new set of service events is generated,
+-- update the accumulator.
+foldRestaurant :: (PrimMonad m) => (a -> [Service] -> MC m a)
+                                -> a
+                                -> [CustomerEvent]
+                                -> Restaurant
+                                -> MC m a
+foldRestaurant f a []     r = finishServing r >>= f a
+foldRestaurant f a (c:cs) r = do
+    (ss,r') <- processEvent c r
+    a' <- f a ss
+    foldRestaurant f a' cs r'
+
+
+-- | Compute a summary of the total waiting times for each customer.
+summarizeService :: (PrimMonad m) => [CustomerEvent] -> Restaurant -> MC m Summary
+summarizeService cs r =
+    foldRestaurant (\s ss -> return $!
+                                 foldl' (flip $ S.insertWith totalTime) s ss)
+                   S.empty cs r
+  where
+    totalTime (Service _ w s) = w+s
+
+
 -- | An infinite stream of customerEvents.  This stream uses its own private 
 -- random number generator (mt19937 is the Mersenne-Twister algorithm).
 customerEvents :: Seed -> [CustomerEvent]
-customerEvents seed = repeatMC customerEvent `evalMC` mt19937 seed
+customerEvents seed = customerEvent `repeatMC` (mt19937 seed)
 
+
 -- | Given a seed for the customers and a seed for the restaurant, run the
 -- simulation.
-simulation :: Seed -> Seed -> [Service]
-simulation customerSeed restaurantSeed =
-    restaurant (customerEvents customerSeed) `evalMC` mt19937 restaurantSeed
+simulation :: Seed -> Seed -> Int -> Summary
+simulation customerSeed restaurantSeed n = let
+    cs = take n $ customerEvents customerSeed
+    r  = emptyRestaurant
+    in evalMC (summarizeService cs r) (mt19937 restaurantSeed)
 
--- | Compute a summary of the total waitings time for each customer.
-summarize :: [Service] -> Summary
-summarize = summary . map totalTime
-  where 
-    totalTime (Service _ w s) = w+s
-    
 -- | Run the program
-main = 
+main =
     let customerSeed    = 0
         restaurantSeed  = 100
         numTransactions = 100000
-        results         = summarize $ take numTransactions $ 
-                              simulation  customerSeed restaurantSeed
+        results         = simulation customerSeed restaurantSeed numTransactions
     in do
         putStrLn ""
         putStrLn "Total Service Time:"
diff --git a/lib/Control/Monad/MC.hs b/lib/Control/Monad/MC.hs
--- a/lib/Control/Monad/MC.hs
+++ b/lib/Control/Monad/MC.hs
@@ -6,9 +6,7 @@
 -- Maintainer : Patrick Perry <patperry@gmail.com>
 -- Stability  : experimental
 --
--- A monad and monad transformer for monte carlo computations.  Currently,
--- the default is the GNU Scientific Library-based implementation, but this
--- may change in the future.
+-- A monad and monad transformer for Monte Carlo computations.
 --
 
 module Control.Monad.MC (
diff --git a/lib/Control/Monad/MC/Base.hs b/lib/Control/Monad/MC/Base.hs
deleted file mode 100644
--- a/lib/Control/Monad/MC/Base.hs
+++ /dev/null
@@ -1,181 +0,0 @@
-{-# LANGUAGE TypeFamilies, PolyKinds #-}
------------------------------------------------------------------------------
--- |
--- Module     : Control.Monad.MC.Base
--- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
--- License    : BSD3
--- Maintainer : Patrick Perry <patperry@gmail.com>
--- Stability  : experimental
---
-
-module Control.Monad.MC.Base
-    where
-
-import Control.Monad
-import qualified Control.Monad.MC.GSLBase as GSL
-
-import qualified Data.Vector.Storable as VS
-
-class HasRNG m where
-    -- | The random number generator type for the monad.
-    type RNG m
-
-class (Monad m, HasRNG m) => MonadMC m where
-    -- | Get the current random number generator.
-    getRNG :: m (RNG m)
-
-    -- | Set the current random number generator.
-    setRNG :: RNG m -> m ()
-
-    -- | @uniform a b@ generates a value uniformly distributed in @[a,b)@.
-    uniform :: Double -> Double -> m Double
-
-    -- | @uniformInt n@ generates an integer uniformly in the range @[0,n-1]@.
-    -- It is an error to call this function with a non-positive value.
-    uniformInt :: Int -> m Int
-
-    -- | @normal mu sigma@ generates a Normal random variable with mean
-    -- @mu@ and standard deviation @sigma@.
-    normal :: Double -> Double -> m Double
-
-    -- | @exponential mu@ generates an Exponential variate with mean @mu@.
-    exponential :: Double -> m Double
-
-    -- | @levy c alpha@ gets a Levy alpha-stable variate with scale @c@ and
-    -- exponent @alpha@.  The algorithm only works for @0 < alpha <= 2@.
-    levy :: Double -> Double -> m Double
-
-    -- | @levySkew c alpha beta @ gets a skew Levy alpha-stable variate
-    -- with scale @c@, exponent @alpha@, and skewness @beta@.  The skew
-    -- parameter must lie in the range @[-1,1]@.  The algorithm only works
-    -- for @0 < alpha <= 2@.
-    levySkew :: Double -> Double -> Double -> m Double
-
-    -- | @poisson mu@ generates a Poisson random variable with mean @mu@.
-    poisson :: Double -> m Int
-
-    -- | @cauchy a@ generates a Cauchy random variable with scale
-    -- parameter @a@.
-    cauchy :: Double -> m Double
-
-    -- | @beta a b@ generates a beta random variable with
-    -- parameters @a@ and @b@.
-    beta :: Double -> Double -> m Double
-
-    -- | @logistic a@ generates a logistic random variable with
-    -- parameter @a@.
-    logistic :: Double -> m Double
-
-    -- | @pareto a b@ generates a Pareto random variable with
-    -- exponent @a@ and scale @b@.
-    pareto :: Double -> Double -> m Double
-
-    -- | @weibull a b@ generates a Weibull random variable with
-    -- scale @a@ and exponent @b@.
-    weibull :: Double -> Double -> m Double
-
-    -- | @gamma a b@ generates a gamma random variable with
-    -- parameters @a@ and @b@.
-    gamma :: Double -> Double -> m Double
-
-    -- | @multinomial n ps@ generates a multinomial random
-    -- variable with parameters @ps@ formed by @n@ trials.
-    multinomial :: Int -> VS.Vector Double -> m (VS.Vector Int)
-
-    -- | @dirichlet alphas@ generates a Dirichlet random variable
-    -- with parameters @alphas@.
-    dirichlet :: VS.Vector Double -> m (VS.Vector Double)
-
-    -- | Get the baton from the Monte Carlo monad without performing any
-    -- computations.  Useful but dangerous.
-    unsafeInterleaveMC :: m a -> m a
-
-
--- | Generate 'True' events with the given probability
-bernoulli :: (MonadMC m) => Double -> m Bool
-bernoulli p = liftM (< p) $ uniform 0 1
-{-# INLINE bernoulli #-}
-
-------------------------------- Instances -----------------------------------
-
-instance HasRNG GSL.MC where
-    type RNG GSL.MC = GSL.RNG
-
-instance MonadMC GSL.MC where
-    getRNG = GSL.getRNG
-    {-# INLINE getRNG #-}
-    setRNG = GSL.setRNG
-    {-# INLINE setRNG #-}
-    uniform = GSL.uniform
-    {-# INLINE uniform #-}
-    uniformInt = GSL.uniformInt
-    {-# INLINE uniformInt #-}
-    normal = GSL.normal
-    {-# INLINE normal #-}
-    exponential = GSL.exponential
-    {-# INLINE exponential #-}
-    levy = GSL.levy
-    {-# INLINE levy #-}
-    levySkew = GSL.levySkew
-    {-# INLINE levySkew #-}
-    poisson = GSL.poisson
-    {-# INLINE poisson #-}
-    cauchy = GSL.cauchy
-    {-# INLINE cauchy #-}
-    beta = GSL.beta
-    {-# INLINE beta #-}
-    logistic = GSL.logistic
-    {-# INLINE logistic #-}
-    pareto = GSL.pareto
-    {-# INLINE pareto #-}
-    weibull = GSL.weibull
-    {-# INLINE weibull #-}
-    gamma = GSL.gamma
-    {-# INLINE gamma #-}
-    multinomial = GSL.multinomial
-    {-# INLINE multinomial #-}
-    dirichlet = GSL.dirichlet
-    {-# INLINE dirichlet #-}
-    unsafeInterleaveMC = GSL.unsafeInterleaveMC
-    {-# INLINE unsafeInterleaveMC #-}
-
-instance (Monad m) => HasRNG (GSL.MCT m) where
-    type RNG (GSL.MCT m) = GSL.RNG
-
-instance (Monad m) => MonadMC (GSL.MCT m) where
-    getRNG = GSL.liftMCT GSL.getRNG
-    {-# INLINE getRNG #-}
-    setRNG r = GSL.liftMCT $ GSL.setRNG r
-    {-# INLINE setRNG #-}
-    uniform a b = GSL.liftMCT $ GSL.uniform a b
-    {-# INLINE uniform #-}
-    uniformInt n = GSL.liftMCT $ GSL.uniformInt n
-    {-# INLINE uniformInt #-}
-    normal mu sigma = GSL.liftMCT $ GSL.normal mu sigma
-    {-# INLINE normal #-}
-    exponential mu = GSL.liftMCT $ GSL.exponential mu
-    {-# INLINE exponential #-}
-    levy c alpha = GSL.liftMCT $ GSL.levy c alpha
-    {-# INLINE levy #-}
-    levySkew c alpha beta_ = GSL.liftMCT $ GSL.levySkew c alpha beta_
-    {-# INLINE levySkew #-}
-    poisson mu = GSL.liftMCT $ GSL.poisson mu
-    {-# INLINE poisson #-}
-    cauchy a = GSL.liftMCT $ GSL.cauchy a
-    {-# INLINE cauchy #-}
-    beta a b = GSL.liftMCT $ GSL.beta a b
-    {-# INLINE beta #-}
-    logistic a = GSL.liftMCT $ GSL.logistic a
-    {-# INLINE logistic #-}
-    pareto a b = GSL.liftMCT $ GSL.pareto a b
-    {-# INLINE pareto #-}
-    weibull a b = GSL.liftMCT $ GSL.weibull a b
-    {-# INLINE weibull #-}
-    gamma a b = GSL.liftMCT $ GSL.gamma a b
-    {-# INLINE gamma #-}
-    multinomial n ps = GSL.liftMCT $ GSL.multinomial n ps
-    {-# INLINE multinomial #-}
-    dirichlet alphas = GSL.liftMCT $ GSL.dirichlet alphas
-    {-# INLINE dirichlet #-}
-    unsafeInterleaveMC = GSL.unsafeInterleaveMCT
-    {-# INLINE unsafeInterleaveMC #-}
diff --git a/lib/Control/Monad/MC/Class.hs b/lib/Control/Monad/MC/Class.hs
deleted file mode 100644
--- a/lib/Control/Monad/MC/Class.hs
+++ /dev/null
@@ -1,27 +0,0 @@
------------------------------------------------------------------------------
--- |
--- Module     : Control.Monad.MC.Class
--- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
--- License    : BSD3
--- Maintainer : Patrick Perry <patperry@gmail.com>
--- Stability  : experimental
---
--- The abstract MonadMC interface and utility functions for Monte Carlo
--- computations.
---
-
-module Control.Monad.MC.Class (
-    -- * The Monte Carlo monad type class
-    HasRNG(..),
-    MonadMC(..),
-
-    -- * Random distributions
-    bernoulli,
-
-    module Control.Monad.MC.Sample,
-    module Control.Monad.MC.Repeat,
-    ) where
-
-import Control.Monad.MC.Base
-import Control.Monad.MC.Sample
-import Control.Monad.MC.Repeat
diff --git a/lib/Control/Monad/MC/GSL.hs b/lib/Control/Monad/MC/GSL.hs
--- a/lib/Control/Monad/MC/GSL.hs
+++ b/lib/Control/Monad/MC/GSL.hs
@@ -10,34 +10,11 @@
 -- of the functions in the GNU Scientific Library.
 
 module Control.Monad.MC.GSL (
-    -- * The Monte Carlo monad
-    MC,
-    runMC,
-    evalMC,
-    execMC,
-
-    -- * The Monte Carlo monad transformer
-    MCT,
-    runMCT,
-    evalMCT,
-    execMCT,
-    liftMCT,
-
-    -- * Pure random number generator creation
-    RNG,
-    Seed,
-    mt19937,
-    mt19937WithState,
-    rngName,
-    rngSize,
-    rngState,
-
-    -- * Overloaded Monte Carlo monad interface
-    module Control.Monad.MC.Class,
-
+    module Control.Monad.MC.GSLBase,
+    module Control.Monad.MC.Sample,
+    module Control.Monad.MC.Repeat,
     ) where
 
-import Control.Monad.MC.GSLBase ( MC, runMC, evalMC, execMC,
-    MCT, runMCT, evalMCT, execMCT, liftMCT, RNG, Seed, mt19937, mt19937WithState,
-    rngName, rngSize, rngState )
-import Control.Monad.MC.Class hiding ( RNG )
+import Control.Monad.MC.GSLBase
+import Control.Monad.MC.Sample
+import Control.Monad.MC.Repeat
diff --git a/lib/Control/Monad/MC/GSLBase.hs b/lib/Control/Monad/MC/GSLBase.hs
--- a/lib/Control/Monad/MC/GSLBase.hs
+++ b/lib/Control/Monad/MC/GSLBase.hs
@@ -1,4 +1,4 @@
-{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses, UndecidableInstances #-}
+{-# LANGUAGE DeriveDataTypeable, RankNTypes, TypeFamilies #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module     : Control.Monad.MC.GSLBase
@@ -9,99 +9,99 @@
 --
 
 module Control.Monad.MC.GSLBase (
-    -- * The Monte Carlo monad
+    -- * Monte Carlo monad transformer
     MC(..),
-    runMC,
+    STMC,
+    IOMC,
     evalMC,
-    execMC,
-    unsafeInterleaveMC,
 
-    -- * The Monte Carlo monad transformer
-    MCT(..),
-    runMCT,
-    evalMCT,
-    execMCT,
-    unsafeInterleaveMCT,
-    liftMCT,
-
-    -- * Pure random number generator creation
+    -- * Random number generator
+    -- ** Types
     RNG,
+    IORNG,
+    STRNG,
     Seed,
+    -- ** Creation
     mt19937,
     mt19937WithState,
-    rngName,
-    rngSize,
-    rngState,
-
-    -- * Getting and setting the random number generator
-    getRNG,
-    setRNG,
+    -- ** State
+    getRNGName,
+    getRNGSize,
+    getRNGState,
+    setRNGState,
 
     -- * Random number distributions
+    -- ** Uniform
     uniform,
     uniformInt,
+    -- ** Continuous
     normal,
     exponential,
+    gamma,
+    cauchy,
     levy,
     levySkew,
-    poisson,
-    cauchy,
-    beta,
-    logistic,
     pareto,
     weibull,
-    gamma,
-    multinomial,
+    logistic,
+    beta,
+    -- ** Discrete
+    bernoulli,
+    poisson,
+    -- ** Multivariate
     dirichlet,
+    multinomial,
     ) where
 
-import Control.Applicative      ( Applicative(..), (<$>) )
-import Control.Monad            ( ap, liftM, MonadPlus(..) )
-import Control.Monad.Cont       ( MonadCont(..) )
-import Control.Monad.Error      ( MonadError(..) )
-import Control.Monad.Reader     ( MonadReader(..) )
-import Control.Monad.State      ( MonadState(..) )
-import Control.Monad.Writer     ( MonadWriter(..) )
-import Control.Monad.Trans      ( MonadTrans(..), MonadIO(..) )
-import Data.Word
-import System.IO.Unsafe         ( unsafePerformIO, unsafeInterleaveIO )
+import Control.Applicative       ( Applicative(..) )
+import Control.Monad             ( liftM )
+import Control.Monad.Fix         ( MonadFix(..) )
+import Control.Monad.IO.Class    ( MonadIO(..) )
+import Control.Monad.ST          ( ST, runST )
+import Control.Monad.Primitive   ( PrimMonad(..), unsafePrimToPrim )
+import Control.Monad.Trans.Class ( MonadTrans(..) )
+import Data.Typeable             ( Typeable )
+import Data.Word                 ( Word8, Word64 )
 
 import qualified Data.Vector.Storable as VS
 
 import qualified GSL.Random.Gen as GSL
 import GSL.Random.Dist
 
--- | A Monte Carlo monad with an internal random number generator.
-newtype MC a = MC (GSL.RNG -> IO a)
 
--- | Run this Monte Carlo monad with the given initial random number generator,
--- getting the result and the new random number generator.
-runMC :: MC a -> RNG -> (a, RNG)
-runMC (MC g) (RNG r) = unsafePerformIO $ do
-    r' <- GSL.cloneRNG r
-    a  <- g r'
-    return (a,RNG r')
-{-# NOINLINE runMC #-}
+-- | A Monte Carlo monad transformer.  This type provides access
+-- to a random number generator while allowing operations in a
+-- base monad, @m@.
+newtype MC m a = MC { runMC :: RNG (PrimState m) -> m a }
 
--- | Evaluate this Monte Carlo monad and throw away the final random number
--- generator.  Very much like @fst@ composed with @runMC@.
-evalMC :: MC a -> RNG -> a
-evalMC g r = fst $ runMC g r
 
--- | Exicute this Monte Carlo monad and return the final random number
--- generator.  Very much like @snd@ composed with @runMC@.
-execMC :: MC a -> RNG -> RNG
-execMC g r = snd $ runMC g r
+-- | Type alias for when the base monad is 'ST'.
+type STMC s a = MC (ST s) a
 
-unsafeInterleaveMC :: MC a -> MC a
-unsafeInterleaveMC (MC m) = MC $ \r ->
-    unsafeInterleaveIO (m r)
+-- | Type alias for when the base monad is 'IO'.
+type IOMC a = MC IO a
 
-instance Functor MC where
-    fmap f (MC m) = MC $ \r ->
-        fmap f (m r)
 
-instance Monad MC where
+-- | Evaluate the result of a Monte Carlo computation using the given
+-- random number generator.
+evalMC :: (forall s. STMC s a) -> (forall s. ST s (STRNG s)) -> a
+evalMC ma mr = runST $ do
+    r <- mr
+    runMC ma r
+
+
+instance (Functor m) => Functor (MC m) where
+    fmap f (MC m) = MC $ \r -> fmap f (m r)
+    {-# INLINE fmap #-}
+
+instance (Applicative m) => Applicative (MC m) where
+    pure a = MC $ \_ -> pure a
+    {-# INLINE pure #-}
+
+    (MC gf) <*> (MC ga) = MC $ \r -> (gf r) <*> (ga r)
+    {-# INLINE (<*>) #-}
+
+instance (Monad m) => Monad (MC m) where
     return a = MC $ \_ -> return a
     {-# INLINE return #-}
 
@@ -111,250 +111,163 @@
             in m' r
     {-# INLINE (>>=) #-}
 
-    fail s = MC $ \_ -> fail s
-    {-# INLINE fail #-}
-
-instance Applicative MC where
-    pure  = return
-    (<*>) = ap
-
--- | A parameterizable Monte Carlo monad for encapsulating an inner
--- monad.
-newtype MCT m a = MCT (GSL.RNG -> IO (m a))
-
--- | Similar to 'runMC'.
-runMCT :: (Monad m) => MCT m a -> RNG -> m (a,RNG)
-runMCT (MCT g) (RNG r) = unsafePerformIO $ do
-    r' <- GSL.cloneRNG r
-    ma <- g r'
-    return (ma >>= \a -> return (a, RNG r'))
-{-# NOINLINE runMCT #-}
-
--- | Similar to 'evalMC'.
-evalMCT :: (Monad m) => MCT m a -> RNG -> m a
-evalMCT g r = do
-    ~(a,_) <- runMCT g r
-    return a
-
--- | Similar to 'execMC'.
-execMCT :: (Monad m) => MCT m a -> RNG -> m RNG
-execMCT g r = do
-    ~(_,r') <- runMCT g r
-    return r'
-
--- | Take a Monte Carlo computations and lift it to an MCT computation.
-liftMCT :: (Monad m) => MC a -> MCT m a
-liftMCT (MC g) = MCT $ \r -> do
-    a <- g r
-    return (return a)
-{-# INLINE liftMCT #-}
-
-unsafeInterleaveMCT :: (Monad m) => MCT m a -> MCT m a
-unsafeInterleaveMCT (MCT g) = MCT $ \r ->
-    unsafeInterleaveIO (g r)
-{-# INLINE unsafeInterleaveMCT #-}
-
-instance (Monad m) => Functor (MCT m) where
-    fmap f (MCT g) = MCT $ \r -> do
-        ma <- g r
-        return (ma >>= return . f)
-    {-# INLINE fmap #-}
-
-instance (Monad m) => Monad (MCT m) where
-    return a = MCT $ \_ -> return (return a)
-    {-# INLINE return #-}
-
-    (MCT g) >>= k =
-        MCT $ \r -> do
-            ma <- g r
-            return $ ma >>= \a ->
-                let (MCT m') = k a
-                in unsafePerformIO $ m' r
-    {-# NOINLINE (>>=) #-}
-
-    fail str = MCT $ \_ -> fail str
+    fail msg = MC $ \_ -> fail msg
     {-# INLINE fail #-}
 
-instance (Monad m) => Applicative (MCT m) where
-    pure  = return
-    (<*>) = ap
-
-instance (MonadPlus m) => MonadPlus (MCT m) where
-    mzero = MCT $ \_ -> mzero
-    {-# INLINE mzero #-}
+instance (MonadFix m) => MonadFix (MC m) where
+    mfix f = MC $ \r -> mfix $ flip (runMC . f) r
+    {-# INLINE mfix #-}
 
-    (MCT m) `mplus` (MCT n) =
-        MCT $ \r -> do
-            r' <- GSL.cloneRNG r
-            mr <- m r
-            nr <- n r'
-            return (mr `mplus` nr)
+instance (MonadIO m) => MonadIO (MC m) where
+    liftIO io = MC $ \_ -> liftIO io
+    {-# INLINE liftIO #-}
 
-instance MonadTrans MCT where
-    lift m = MCT $ \_ -> return m
+instance MonadTrans MC where
+    lift m = MC $ \_ -> m
     {-# INLINE lift #-}
 
-instance (MonadCont m) => MonadCont (MCT m) where
-    callCC f = MCT $ \r ->
-        return $ callCC $ \k ->
-            let (MCT m) = f (\a -> MCT $ \_ -> return (k a))
-            in unsafePerformIO (m r)
-    {-# NOINLINE callCC #-}
 
-instance (MonadError e m) => MonadError e (MCT m) where
-    throwError             = lift . throwError
-    {-# INLINE throwError #-}
-
-    (MCT g) `catchError` h = MCT $ \r -> do
-        ma <- g r
-        return $ ma `catchError` \e ->
-            let (MCT m') = h e
-            in unsafePerformIO (m' r)
-    {-# NOINLINE catchError #-}
-
-instance (MonadIO m) => MonadIO (MCT m) where
-    liftIO = lift . liftIO
-    {-# INLINE liftIO #-}
-
-instance (MonadReader r m) => MonadReader r (MCT m) where
-    ask              = lift ask
-    {-# INLINE ask #-}
-
-    local f (MCT g) = MCT $ \r -> do
-        ma <- g r
-        return $ local f ma
-    {-# INLINE local #-}
-
-instance (MonadState s m) => MonadState s (MCT m) where
-    get = lift get
-    {-# INLINE get #-}
-
-    put = lift . put
-    {-# INLINE put #-}
-
-instance (MonadWriter w m) => MonadWriter w (MCT m) where
-    tell           = lift . tell
-    {-# INLINE tell #-}
-
-    listen (MCT g) = MCT $ \r -> do
-        ma <- g r
-        return (listen ma)
-    {-# INLINE listen #-}
+---------------------------- Random Number Generators -----------------------
 
-    pass (MCT g) = MCT $ \r -> do
-        maf <- g r
-        return (pass maf)
-    {-# INLINE pass #-}
+-- | The random number generator type.
+newtype RNG s = RNG GSL.RNG
+    deriving(Eq, Show, Typeable)
 
----------------------------- Random Number Generators -----------------------
+-- | A shorter name for RNG in the 'IO' monad.
+type IORNG = RNG (PrimState IO)
 
--- | The random number generator type associated with 'MC' and 'MCT'.
-newtype RNG = RNG GSL.RNG
+-- | A shorter name for RNG in the 'ST' monad.
+type STRNG s = RNG (PrimState (ST s))
 
 -- | The seed type for the random number generators.
 type Seed = Word64
 
+
 -- | Get the name of the random number generator algorithm.
-rngName :: RNG -> String
-rngName (RNG r) = unsafePerformIO $ GSL.getName r
-{-# NOINLINE rngName #-}
+getRNGName :: (PrimMonad m) => RNG (PrimState m) -> m String
+getRNGName (RNG r) = unsafePrimToPrim $ GSL.getName r
 
 -- | Get the size of the generator state, in bytes.
-rngSize :: RNG -> Int
-rngSize (RNG r) = fromIntegral $ unsafePerformIO $ GSL.getSize r
-{-# NOINLINE rngSize #-}
+getRNGSize :: (PrimMonad m) => RNG (PrimState m) -> m Int
+getRNGSize (RNG r) = liftM fromIntegral $ unsafePrimToPrim $ GSL.getSize r
 
 -- | Get the state of the generator.
-rngState :: RNG -> [Word8]
-rngState (RNG r) = unsafePerformIO $ GSL.getState r
-{-# NOINLINE rngState #-}
-
-getRNG :: MC RNG
-getRNG = MC (\r -> liftM RNG $ GSL.cloneRNG r)
-{-# INLINE getRNG #-}
+getRNGState :: (PrimMonad m) => RNG (PrimState m) -> m [Word8]
+getRNGState (RNG r) = unsafePrimToPrim $ GSL.getState r
 
-setRNG :: RNG -> MC ()
-setRNG (RNG r') = MC $ \r -> GSL.copyRNG r r'
-{-# INLINE setRNG #-}
+-- | Set the state of the generator.
+setRNGState :: (PrimMonad m) => RNG (PrimState m) -> [Word8] -> m ()
+setRNGState (RNG r) xs = unsafePrimToPrim $ GSL.setState r xs
 
--- | Get a Mersenne Twister random number generator seeded with the given
+-- | Create a Mersenne Twister random number generator seeded with the given
 -- value.
-mt19937 :: Seed -> RNG
-mt19937 s = unsafePerformIO $ do
+mt19937 :: (PrimMonad m) => Seed -> m (RNG (PrimState m))
+mt19937 s = unsafePrimToPrim $ do
     r <- GSL.newRNG GSL.mt19937
     GSL.setSeed r s
     return (RNG r)
-{-# NOINLINE mt19937 #-}
 
--- | Get a Mersenne Twister seeded with the given state.
-mt19937WithState :: [Word8] -> RNG
-mt19937WithState xs = unsafePerformIO $ do
+-- | Create a Mersenne Twister seeded with the given state.
+mt19937WithState :: (PrimMonad m) => [Word8] -> m (RNG (PrimState m))
+mt19937WithState xs = unsafePrimToPrim $ do
     r <- GSL.newRNG GSL.mt19937
     GSL.setState r xs
     return (RNG r)
-{-# NOINLINE mt19937WithState #-}
 
 -------------------------- Random Number Distributions ----------------------
 
-uniform :: Double -> Double -> MC Double
+-- | @uniform a b@ generates a value uniformly distributed in @[a,b)@.
+uniform :: (PrimMonad m) => Double -> Double -> MC m Double
 uniform 0 1 = liftRan0 GSL.getUniform
 uniform a b = liftRan2 getFlat a b
 
-uniformInt :: Int -> MC Int
+-- | @uniformInt n@ generates an integer uniformly in the range @[0,n-1]@.
+-- It is an error to call this function with a non-positive value.
+uniformInt :: (PrimMonad m) => Int -> MC m Int
 uniformInt = liftRan1 GSL.getUniformInt
 
-normal :: Double -> Double -> MC Double
-normal 0  1     =            liftRan0 getUGaussianRatioMethod
-normal mu 1     = (mu +) <$> liftRan0 getUGaussianRatioMethod
-normal 0  sigma =            liftRan1 getGaussianRatioMethod sigma
-normal mu sigma = (mu +) <$> liftRan1 getGaussianRatioMethod sigma
+-- | @normal mu sigma@ generates a Normal random variable with mean
+-- @mu@ and standard deviation @sigma@.
+normal :: (PrimMonad m) => Double -> Double -> MC m Double
+normal 0  1     =                liftRan0 getUGaussianRatioMethod
+normal mu 1     = liftM (mu +) $ liftRan0 getUGaussianRatioMethod
+normal 0  sigma =                liftRan1 getGaussianRatioMethod sigma
+normal mu sigma = liftM (mu +) $ liftRan1 getGaussianRatioMethod sigma
 
-exponential :: Double -> MC Double
+-- | @exponential mu@ generates an Exponential variate with mean @mu@.
+exponential :: (PrimMonad m) => Double -> MC m Double
 exponential = liftRan1 getExponential
 
-poisson :: Double -> MC Int
+-- | @poisson mu@ generates a Poisson random variable with mean @mu@.
+poisson :: (PrimMonad m) => Double -> MC m Int
 poisson = liftRan1 getPoisson
 
-levy :: Double -> Double -> MC Double
+-- | @levy c alpha@ gets a Levy alpha-stable variate with scale @c@ and
+-- exponent @alpha@.  The algorithm only works for @0 < alpha <= 2@.
+levy :: (PrimMonad m) => Double -> Double -> MC m Double
 levy = liftRan2 getLevy
 
-levySkew :: Double -> Double -> Double -> MC Double
+-- | @levySkew c alpha beta @ gets a skew Levy alpha-stable variate
+-- with scale @c@, exponent @alpha@, and skewness @beta@.  The skew
+-- parameter must lie in the range @[-1,1]@.  The algorithm only works
+-- for @0 < alpha <= 2@.
+levySkew :: (PrimMonad m) => Double -> Double -> Double -> MC m Double
 levySkew = liftRan3 getLevySkew
 
-cauchy :: Double -> MC Double
+-- | @cauchy a@ generates a Cauchy random variable with scale
+-- parameter @a@.
+cauchy :: (PrimMonad m) => Double -> MC m Double
 cauchy = liftRan1 getCauchy
 
-beta :: Double -> Double -> MC Double
+-- | @beta a b@ generates a beta random variable with
+-- parameters @a@ and @b@.
+beta :: (PrimMonad m) => Double -> Double -> MC m Double
 beta = liftRan2 getBeta
 
-logistic :: Double -> MC Double
+-- | @logistic a@ generates a logistic random variable with
+-- parameter @a@.
+logistic :: (PrimMonad m) => Double -> MC m Double
 logistic = liftRan1 getLogistic
 
-pareto :: Double -> Double -> MC Double
+-- | @pareto a b@ generates a Pareto random variable with
+-- exponent @a@ and scale @b@.
+pareto :: (PrimMonad m) => Double -> Double -> MC m Double
 pareto = liftRan2 getPareto
 
-weibull :: Double -> Double -> MC Double
+-- | @weibull a b@ generates a Weibull random variable with
+-- scale @a@ and exponent @b@.
+weibull :: (PrimMonad m) => Double -> Double -> MC m Double
 weibull = liftRan2 getWeibull
 
-gamma :: Double -> Double -> MC Double
+-- | @gamma a b@ generates a gamma random variable with
+-- parameters @a@ and @b@.
+gamma :: (PrimMonad m) => Double -> Double -> MC m Double
 gamma = liftRan2 getGamma
 
-multinomial :: Int -> VS.Vector Double -> MC (VS.Vector Int)
+-- | @multinomial n ps@ generates a multinomial random
+-- variable with parameters @ps@ formed by @n@ trials.
+multinomial :: (PrimMonad m) => Int -> VS.Vector Double -> MC m (VS.Vector Int)
 multinomial = liftRan2 getMultinomial
 
-dirichlet :: VS.Vector Double -> MC (VS.Vector Double)
+-- | @dirichlet alphas@ generates a Dirichlet random variable
+-- with parameters @alphas@.
+dirichlet :: (PrimMonad m) => VS.Vector Double -> MC m (VS.Vector Double)
 dirichlet = liftRan1 getDirichlet
 
+-- | Generate 'True' events with the given probability.
+bernoulli :: (PrimMonad m ) => Double -> MC m Bool
+bernoulli p = liftM (< p) $ uniform 0 1
+{-# INLINE bernoulli #-}
 
 
-liftRan0 :: (GSL.RNG -> IO a) -> MC a
-liftRan0 = MC
+liftRan0 :: (PrimMonad m) => (GSL.RNG -> IO a) -> MC m a
+liftRan0 f = MC $ \(RNG r) -> unsafePrimToPrim $ f r
 
-liftRan1 :: (GSL.RNG -> a -> IO b) -> a -> MC b
-liftRan1 f a = MC $ \r -> f r a
+liftRan1 :: (PrimMonad m) => (GSL.RNG -> a -> IO b) -> a -> MC m b
+liftRan1 f a = MC $ \(RNG r) -> unsafePrimToPrim $ f r a
 
-liftRan2 :: (GSL.RNG -> a -> b -> IO c) -> a -> b -> MC c
-liftRan2 f a b = MC $ \r -> f r a b
+liftRan2 :: (PrimMonad m) => (GSL.RNG -> a -> b -> IO c) -> a -> b -> MC m c
+liftRan2 f a b = MC $ \(RNG r) -> unsafePrimToPrim $ f r a b
 
-liftRan3 :: (GSL.RNG -> a -> b -> c -> IO d) -> a -> b -> c -> MC d
-liftRan3 f a b c = MC $ \r -> f r a b c
+liftRan3 :: (PrimMonad m) => (GSL.RNG -> a -> b -> c -> IO d) -> a -> b -> c -> MC m d
+liftRan3 f a b c = MC $ \(RNG r) -> unsafePrimToPrim $ f r a b c
diff --git a/lib/Control/Monad/MC/Repeat.hs b/lib/Control/Monad/MC/Repeat.hs
--- a/lib/Control/Monad/MC/Repeat.hs
+++ b/lib/Control/Monad/MC/Repeat.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE RankNTypes #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module     : Control.Monad.MC.Repeat
@@ -9,28 +10,50 @@
 
 module Control.Monad.MC.Repeat (
     -- * Repeating computations
+    foldMC,
     repeatMC,
     replicateMC,
     ) where
 
-import Control.Monad.MC.Base
+import Control.Monad.Primitive( PrimMonad )
+import Control.Monad.MC.GSLBase
+import Control.Monad.ST( ST, runST )
+import Control.Monad.ST.Unsafe( unsafeInterleaveST )
 
--- | Produce a lazy infinite list of values from the given Monte Carlo
--- generator.
-repeatMC :: (MonadMC m) => m a -> m [a]
-repeatMC = interleaveSequence . repeat
-{-# INLINE repeatMC #-}
 
+-- | Generate a sequence of replicates and incrementally consume
+-- them via a left fold.
+--
+-- This fold is /not/ strict. The replicate consumer is responsible for
+-- forcing the evaluation of its result to avoid space leaks.
+foldMC :: (PrimMonad m) => (a -> b -> MC m a) -- ^ Replicate consumer.
+                        -> a                  -- ^ Initial state for replicate consumer.
+                        -> Int                -- ^ Number of replicates.
+                        -> MC m b             -- ^ Generator.
+                        -> MC m a
+foldMC f a n mb | n <= 0    = return a
+                | otherwise = do
+    b <- mb
+    a' <- f a b
+    foldMC f a' (n-1) mb
+{-# INLINE foldMC #-}
+
+
+-- | Produce a lazy infinite list of replicates from the given random
+-- number generator and Monte Carlo procedure.
+repeatMC :: (forall s. STMC s a) -> (forall s. ST s (STRNG s)) -> [a]
+repeatMC mc mrng = runST $ do
+    rng <- mrng
+    go $ runMC mc rng
+  where
+    go m = unsafeInterleaveST $ do
+        a  <- m
+        as <- go m
+        return (a:as)
+
+
 -- | Produce a lazy list of the given length using the specified
--- generator.
-replicateMC :: (MonadMC m) => Int -> m a -> m [a]
-replicateMC n = interleaveSequence . replicate n
-{-# INLINE replicateMC #-}
+-- random number genrator and Monte Carlo procedure.
+replicateMC :: Int -> (forall s. STMC s a) -> (forall s. ST s (STRNG s)) -> [a]
+replicateMC n mc mrng = take n $ repeatMC mc mrng
 
-interleaveSequence :: (MonadMC m) => [m a] -> m [a]
-interleaveSequence []     = return []
-interleaveSequence (m:ms) = unsafeInterleaveMC $ do
-    a  <- m
-    as <- interleaveSequence ms
-    return (a:as)
-{-# INLINE interleaveSequence #-}
diff --git a/lib/Control/Monad/MC/Sample.hs b/lib/Control/Monad/MC/Sample.hs
--- a/lib/Control/Monad/MC/Sample.hs
+++ b/lib/Control/Monad/MC/Sample.hs
@@ -8,295 +8,180 @@
 --
 
 module Control.Monad.MC.Sample (
-    -- * Sampling from lists
+    -- * Sampling
+    -- ** Lists
     sample,
     sampleWithWeights,
     sampleSubset,
-    sampleSubset',
     sampleSubsetWithWeights,
-    sampleSubsetWithWeights',
+    shuffle,
 
-    -- * Sampling @Int@s
+    -- ** Ints
     sampleInt,
     sampleIntWithWeights,
     sampleIntSubset,
-    sampleIntSubset',
     sampleIntSubsetWithWeights,
-    sampleIntSubsetWithWeights',
-
-    -- * Shuffling
-    shuffle,
     shuffleInt,
-    shuffleInt',
     ) where
 
-import Control.Monad
-import Control.Monad.ST hiding (unsafeInterleaveST)
-import Control.Monad.ST.Unsafe (unsafeInterleaveST)
-import Control.Monad.MC.Base
-import Control.Monad.MC.Repeat
-import Control.Monad.MC.Walker
+import Control.Monad( forM_, liftM )
+import Control.Monad.Primitive( PrimMonad )
+import Control.Monad.Trans.Class( lift )
 import Data.List( foldl', sort )
-
-import Data.Vector.Unboxed( MVector, Unbox )
 import qualified Data.Vector as BV
-import qualified Data.Vector.Mutable as BMV
 import qualified Data.Vector.Unboxed as V
-import qualified Data.Vector.Generic.Mutable as MV
+import qualified Data.Vector.Unboxed.Mutable as MV
 
+import Control.Monad.MC.GSLBase
+import Control.Monad.MC.Walker
+
+
 -- | @sample xs@ samples a value uniformly from the elements of @xs@.  The
 -- results are undefined if @length xs@ is zero.
-sample :: (MonadMC m) => [a] -> m a
+sample :: (PrimMonad m) => [a] -> MC m a
 sample xs = let
     n = length xs
     in sampleHelp n xs $ sampleInt n
 {-# INLINE sample #-}
 
+
 -- | @sampleWithWeights wxs@ samples a value from the list with the given
 -- weight.
-sampleWithWeights :: (MonadMC m) => [(Double, a)] -> m a
+sampleWithWeights :: (PrimMonad m) => [(Double, a)] -> MC m a
 sampleWithWeights wxs = let
     (ws,xs) = unzip wxs
     n       = length xs
     in sampleHelp n xs $ sampleIntWithWeights ws n
 {-# INLINE sampleWithWeights #-}
 
+
+sampleHelp :: (PrimMonad m) => Int -> [a] -> MC m Int -> MC m a
+sampleHelp _n xs f = let
+    arr = BV.fromList xs
+    in liftM (BV.unsafeIndex arr) f
+{-# INLINE sampleHelp #-}
+
+
 -- | @sampleSubset xs k@ samples a subset of size @k@ from @xs@ by
 -- sampling without replacement.  The return value is a list of length @k@
--- with the elements in the subset in the order that they were sampled.  Note
--- also that the elements are lazily generated.
-sampleSubset :: (MonadMC m) => [a] -> Int -> m [a]
+-- with the elements in the subset in the order that they were sampled.
+sampleSubset :: (PrimMonad m) => [a] -> Int -> MC m [a]
 sampleSubset xs k = let
     n = length xs
-    in sampleListHelp n xs $ sampleIntSubset n k
+    in sampleSubsetHelp n xs $ sampleIntSubset n k
 {-# INLINE sampleSubset #-}
 
--- | Strict version of 'sampleSubset'.
-sampleSubset' :: (MonadMC m) => [a] -> Int -> m [a]
-sampleSubset' xs k = do
-    s <- sampleSubset xs k
-    length s `seq` return s
-{-# INLINE sampleSubset' #-}
 
 -- | Sample a subset of the elements with the given weights.  Return
--- the elements of the subset lazily in the order they were sampled.
-sampleSubsetWithWeights :: (MonadMC m) => [(Double,a)] -> Int -> m [a]
+-- the elements of the subset in the order they were sampled.
+sampleSubsetWithWeights :: (PrimMonad m) => [(Double,a)] -> Int -> MC m [a]
 sampleSubsetWithWeights wxs k = let
     (ws,xs) = unzip wxs
     n = length ws
-    in sampleListHelp n xs $ sampleIntSubsetWithWeights ws n k
+    in sampleSubsetHelp n xs $ sampleIntSubsetWithWeights ws n k
 {-# INLINE sampleSubsetWithWeights #-}
 
--- | Strict version of 'sampleSubsetWithWeights'.
-sampleSubsetWithWeights' :: (MonadMC m) => [(Double,a)] -> Int -> m [a]
-sampleSubsetWithWeights' wxs k = do
-    s <- sampleSubsetWithWeights wxs k
-    length s `seq` return s
-{-# INLINE sampleSubsetWithWeights' #-}
 
-sampleHelp :: (Monad m) => Int -> [a] -> m Int -> m a
-sampleHelp _n xs f = let
-    arr = BV.fromList xs
-    in liftM (BV.unsafeIndex arr) f
-{-# INLINE sampleHelp #-}
-
-sampleHelpU :: (Unbox a, Monad m) => Int -> [a] -> m Int -> m a
-sampleHelpU _n xs f = let
-    arr = V.fromList xs
-    in liftM (V.unsafeIndex arr) f
-{-# INLINE sampleHelpU #-}
-
-{-# RULES "sampleHelp/Double" forall n xs f.
-              sampleHelp n (xs :: [Double]) f = sampleHelpU n xs f #-}
-{-# RULES "sampleHelp/Int" forall n xs f.
-              sampleHelp n (xs :: [Int]) f = sampleHelpU n xs f #-}
-
-sampleListHelp :: (Monad m) => Int -> [a] -> m [Int] -> m [a]
-sampleListHelp _n xs f = let
+sampleSubsetHelp :: (Monad m) => Int -> [a] -> m [Int] -> m [a]
+sampleSubsetHelp _n xs f = let
     arr = BV.fromList xs
     in liftM (map $ BV.unsafeIndex arr) f
-{-# INLINE sampleListHelp #-}
+{-# INLINE sampleSubsetHelp #-}
 
-sampleListHelpU :: (Unbox a, Monad m) => Int -> [a] -> m [Int] -> m [a]
-sampleListHelpU _n xs f = let
-    arr = V.fromList xs
-    in liftM (map $ V.unsafeIndex arr) f
 
-{-# RULES "sampleListHelp/Double" forall n xs f.
-              sampleListHelp n (xs :: [Double]) f = sampleListHelpU n xs f #-}
-{-# RULES "sampleListHelp/Int" forall n xs f.
-              sampleListHelp n (xs :: [Int]) f = sampleListHelpU n xs f #-}
-
 -- | @sampleInt n@ samples integers uniformly from @[ 0..n-1 ]@.  It is an
 -- error to call this function with a non-positive @n@.
-sampleInt :: (MonadMC m) => Int -> m Int
+sampleInt :: (PrimMonad m) => Int -> MC m Int
 sampleInt n | n < 1     = fail "invalid argument"
             | otherwise = uniformInt n
-{-# INLINE sampleInt #-}
 
+
 -- | @sampleIntWithWeights ws n@ samples integers from @[ 0..n-1 ]@ with the
 -- probability of choosing @i@ proportional to @ws !! i@.  The list @ws@ must
 -- have length equal to @n@.  Also, the elements of @ws@ must be non-negative
 -- with at least one nonzero entry.
-sampleIntWithWeights :: (MonadMC m) => [Double] -> Int -> m Int
+sampleIntWithWeights :: (PrimMonad m) => [Double] -> Int -> MC m Int
 sampleIntWithWeights ws n =
     let qjs = computeTable n ws
     in liftM (indexTable qjs) (uniform 0 1)
-{-# INLINE sampleIntWithWeights #-}
 
+
 -- | @sampleIntSubset n k@ samples a subset of size @k@ by sampling without
 -- replacement from the integers @{ 0, ..., n-1 }@.  The return value is a
 -- list of length @k@ with the elements in the subset in the order that they
--- were sampled.  Note also that the elements are lazily generated.
-sampleIntSubset :: (MonadMC m) => Int -> Int -> m [Int]
+-- were sampled.
+sampleIntSubset :: (PrimMonad m) => Int -> Int -> MC m [Int]
 sampleIntSubset n k | k < 0     = fail "negative subset size"
                     | k > n     = fail "subset size is too big"
                     | otherwise = do
-    us <- randomIndices n k
-    return $ runST $ do
-        ints <- MV.new n :: ST s (MVector s Int)
-        sequence_ [ MV.unsafeWrite ints i i | i <- [0 .. n-1] ]
-        sampleIntSubsetHelp ints us (n-1)
+    xs <- lift $ (V.thaw . V.fromList) [ 0..n-1 ]
+    go xs [] n k
   where
-    randomIndices n' k' | k' == 0   = return []
-                        | otherwise = unsafeInterleaveMC $ do
-        u  <- uniformInt n'
-        us <- randomIndices (n'-1) (k'-1)
-        return (u:us)
+    go xs ys n' k' | k' == 0   = return $ reverse ys
+                   | otherwise = do
+        u <- uniformInt n'
+        y <- lift $ do
+                 i <- MV.unsafeRead xs u
+                 j <- MV.unsafeRead xs (n'-1)
+                 MV.unsafeWrite xs u j
+                 return i
+        go xs (y:ys) (n'-1) (k'-1)
 
-    sampleIntSubsetHelp :: MVector m Int -> [Int] -> Int -> ST m [Int]
-    sampleIntSubsetHelp _    []     _  = return []
-    sampleIntSubsetHelp ints (u:us) n' = unsafeInterleaveST $ do
-        i <- MV.unsafeRead ints u
-        MV.unsafeWrite ints u =<< MV.unsafeRead ints n'
-        is <- sampleIntSubsetHelp ints us (n'-1)
-        return (i:is)
-{-# INLINE sampleIntSubset #-}
 
--- | Strict version of 'sampleIntSubset'.
-sampleIntSubset' :: (MonadMC m) => Int -> Int -> m [Int]
-sampleIntSubset' n k = do
-    s <- sampleIntSubset n k
-    length s `seq` return s
-{-# INLINE sampleIntSubset' #-}
-
 -- | @sampleIntSubsetWithWeights ws n k@ samplea size @k@ subset of
 -- @{ 0, ..., n-1 }@ with the given weights by sampling elements without
--- replacement.  It returns the elements of the subset lazily in the order
+-- replacement.  It returns the elements of the subset in the order
 -- they were sampled.
-sampleIntSubsetWithWeights :: (MonadMC m) => [Double] -> Int -> Int -> m [Int]
-sampleIntSubsetWithWeights ws n k = let
-    w_sum0 = foldl' (+) 0 $ take n ws
-    wjs = [ (w / w_sum0, j) | (w,j) <- reverse $ sort $ zip ws [ 0..n-1 ] ]
+sampleIntSubsetWithWeights :: (PrimMonad m) => [Double] -> Int -> Int -> MC m [Int]
+sampleIntSubsetWithWeights ws n k | k < 0     = fail "negative subset size"
+                                  | k > n     = fail "subset size is too big"
+                                  | otherwise = let
+    wsum = foldl' (+) 0 $ take n ws
+    wjs  = [ (w / wsum, j) | (w,j) <- reverse $ sort $ zip ws [ 0..n-1 ] ]
     in do
-        us <- replicateMC k $ uniform 0 1
-        return $ runST $ do
-            ints <- MV.new n :: ST s (MVector s (Double,Int))
-            sequence_ [ MV.unsafeWrite ints i wj | (i,wj) <- zip [ 0.. ] wjs ]
-            go ints n 1 us
+        xs <- lift $ (V.thaw . V.fromList) wjs
+        go xs wsum [] n k
   where
-    go :: MVector m (Double, Int) -> Int -> Double -> [Double] -> ST m [Int]
-    go ints n' w_sum us | null us   = return []
-                        | otherwise = let
-        target = head us * w_sum
-        in unsafeInterleaveST $ do
-            (i,(w,j)) <- findTarget ints n' target 0 0
-            shiftDown ints (i+1) (n'-1)
-            let w_sum' = w_sum - w
-                n''    = n' - 1
-                us'    = tail us
-            js <- go ints n'' w_sum' us'
-            return $ j:js
+    go xs wsum' ys n' k' | k' == 0   = return $ reverse ys
+                         | otherwise = do
+        target <- uniform 0 wsum'
+        (w,y) <- lift $ do
+                     (i,wj) <- findTarget xs n' target 0 0
+                     shiftDown xs (i+1) (n'-1)
+                     return wj
+        let wsum'' = wsum' - w
+            ys'    = y:ys
+            n''    = n' - 1
+            k''    = k' - 1
+        go xs wsum'' ys' n'' k''
 
-    findTarget :: MVector m (Double, Int)
-                    -> Int -> Double -> Int -> Double -> ST m (Int, (Double, Int))
-    findTarget ints n' target i acc
+    findTarget xs n' target i acc
         | i == n' - 1 = do
-            wj <- MV.unsafeRead ints i
+            wj <- MV.unsafeRead xs i
             return (i,wj)
         | otherwise = do
-            (w,j) <- MV.unsafeRead ints i
+            (w,j) <- MV.unsafeRead xs i
             let acc' = acc + w
             if target <= acc'
                 then return (i,(w,j))
-                else findTarget ints n' target (i+1) acc'
+                else findTarget xs n' target (i+1) acc'
 
-    shiftDown :: MVector m (Double, Int) -> Int -> Int -> ST m ()
-    shiftDown ints from to =
+    shiftDown xs from to =
         forM_ [ from..to ] $ \i -> do
-            wj <- MV.unsafeRead ints i
-            MV.unsafeWrite ints (i-1) wj
+            wj <- MV.unsafeRead xs i
+            MV.unsafeWrite xs (i-1) wj
 
-{-# INLINE sampleIntSubsetWithWeights #-}
 
--- | Strict version of 'sampleIntSubsetWithWeights'.
-sampleIntSubsetWithWeights' :: (MonadMC m) => [Double] -> Int -> Int -> m [Int]
-sampleIntSubsetWithWeights' ws n k = do
-    s <- sampleIntSubsetWithWeights ws n k
-    length s `seq` return s
-{-# INLINE sampleIntSubsetWithWeights' #-}
-
 -- | @shuffle xs@ randomly permutes the list @xs@ and returns
 -- the result.  All permutations of the elements of @xs@ are equally
 -- likely.
-shuffle :: (MonadMC m) => [a] -> m [a]
+shuffle :: (PrimMonad m) => [a] -> MC m [a]
 shuffle xs = let
-    n = length xs
-    in shuffleInt n >>= \swaps -> (return . BV.toList) $ BV.create $ do
-           marr <- MV.new n
-           zipWithM_ (MV.unsafeWrite marr) [0 .. n-1] xs
-           mapM_ (swap marr) swaps
-           return marr
-  where
-    swap :: BMV.MVector s a -> (Int,Int) -> ST s ()
-    swap marr (i,j) | i == j    = return ()
-                    | otherwise = do
-        x <- MV.unsafeRead marr i
-        y <- MV.unsafeRead marr j
-        MV.unsafeWrite marr i y
-        MV.unsafeWrite marr j x
-{-# INLINE shuffle #-}
-
-shuffleU :: (Unbox a, MonadMC m) => [a] -> m [a]
-shuffleU xs = let
-    n = length xs
-    in shuffleInt n >>= \swaps -> (return . V.toList) $ V.create $ do
-           marr <- MV.new n
-           zipWithM_ (MV.unsafeWrite marr) [0 .. n-1] xs
-           mapM_ (swap marr) swaps
-           return marr
-  where
-    swap :: (Unbox a) => MVector s a -> (Int,Int) -> ST s ()
-    swap marr (i,j) | i == j    = return ()
-                    | otherwise = do
-        x <- MV.unsafeRead marr i
-        y <- MV.unsafeRead marr j
-        MV.unsafeWrite marr i y
-        MV.unsafeWrite marr j x
-{-# INLINE shuffleU #-}
-
-{-# RULES "shuffle/Double" forall xs.
-              shuffle (xs :: [Double]) = shuffleU xs #-}
-{-# RULES "shuffle/Int" forall xs.
-              shuffle (xs :: [Int]) = shuffleU xs #-}
+    n   = length xs
+    mis = liftM BV.fromList $ shuffleInt n
+    in liftM (BV.toList . BV.unsafeBackpermute (BV.fromList xs)) mis
 
 
--- | @shuffleInt n@ generates a sequence of swaps equivalent to a
--- uniformly-chosen random permutatation of the integers @{0, ..., n-1}@.
--- For an input of @n@, there are @n-1@ swaps, which are lazily generated.
-shuffleInt :: (MonadMC m) => Int -> m [(Int,Int)]
-shuffleInt n =
-    let shuffleIntHelp i | i <= 1    = return []
-                         | otherwise = unsafeInterleaveMC $ do
-            j   <- uniformInt i
-            ijs <- shuffleIntHelp (i-1)
-            return $ (i-1,j):ijs in
-    shuffleIntHelp n
-{-# INLINE shuffleInt #-}
-
--- | Strict version of 'shuffleInt'.
-shuffleInt' :: (MonadMC m) => Int -> m [(Int,Int)]
-shuffleInt' n = do
-    ss <- shuffleInt n
-    length ss `seq` return ss
-{-# INLINE shuffleInt' #-}
+-- | @shuffleInt n@ randomly permutes the elements of the list @[ 0..n-1 ]@.
+shuffleInt :: (PrimMonad m) => Int -> MC m [Int]
+shuffleInt n = sampleIntSubset n n
diff --git a/lib/Data/Summary.hs b/lib/Data/Summary.hs
deleted file mode 100644
--- a/lib/Data/Summary.hs
+++ /dev/null
@@ -1,15 +0,0 @@
------------------------------------------------------------------------------
--- |
--- Module     : Data.Summary
--- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
--- License    : BSD3
--- Maintainer : Patrick Perry <patperry@gmail.com>
--- Stability  : experimental
---
--- Summary Statistics
---
-module Data.Summary (
-    module Data.Summary.Double
-    ) where
-
-import Data.Summary.Double
diff --git a/lib/Data/Summary/Bool.hs b/lib/Data/Summary/Bool.hs
--- a/lib/Data/Summary/Bool.hs
+++ b/lib/Data/Summary/Bool.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE DeriveDataTypeable #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module     : Data.Summary.Bool
@@ -10,26 +11,48 @@
 --
 
 module Data.Summary.Bool (
-    -- * The @Summary@ data type
+    -- * Summary type
     Summary,
-    summary,
-    update,
 
-    -- * @Summary@ properties
-    sampleSize,
-    count,
-    sampleMean,
-    sampleSE,
-    sampleCI,
+    -- * Properties
+    -- ** Sum
+    size,
+    sum,
+    -- ** Mean
+    mean,
+    meanSE,
+    meanCI,
 
+    -- * Construction
+    empty,
+    singleton,
+
+    -- * Insertion
+    insert,
+    insertWith,
+
+    -- * Combination
+    union,
+    unions,
+
+    -- * Conversion
+    -- ** Lists
+    fromList,
+    fromListWith,
+
+    -- ** Statistics
+    toStats,
+    fromStats,
+
     ) where
 
-import Control.DeepSeq
+import Prelude hiding (sum)
 import Data.List( foldl' )
-import Data.Monoid
-import Text.Printf
+import Data.Monoid( Monoid(..) )
+import Data.Data( Data, Typeable )
+import Text.Printf( printf )
 
-import Data.Summary.Utils
+import Data.Summary.Utils( interval )
 
 
 -- | A type for storing summary statistics for a data set of
@@ -37,63 +60,86 @@
 -- of 'True' events and gives estimates for the success
 -- probability.  'True' is interpreted as a one, and 'False'
 -- is interpreted as a zero.
-data Summary = S {-# UNPACK #-} !Int  -- sample size
-                 {-# UNPACK #-} !Int  -- number of successes
+data Summary = S {-# UNPACK #-} !Int -- number of observations
+                 {-# UNPACK #-} !Int -- number of True values
+    deriving(Eq, Data, Typeable)
 
 instance Show Summary where
-    show s@(S n c) =
-        printf "    sample size: %d" n
-        ++ printf "\n      successes: %d" c
-        ++ printf "\n     proportion: %g" (sampleMean s)
-        ++ printf "\n             SE: %g" (sampleSE s)
+    show s@(S n x) =
+             printf "    sample size: %d" n
+        ++ printf "\n      successes: %d" x
+        ++ printf "\n     proportion: %g" (mean s)
+        ++ printf "\n             SE: %g" (meanSE s)
         ++ printf "\n         99%% CI: (%g, %g)" c1 c2
-      where (c1,c2) = sampleCI 0.99 s
+      where (c1,c2) = meanCI 0.99 s
 
 instance Monoid Summary where
     mempty = empty
     mappend = union
 
-instance NFData Summary
+-- | Number of observations.
+size :: Summary -> Int
+size (S n _) = n
 
+-- | Number of 'True' values.
+sum :: Summary -> Int
+sum (S _ x) = x
 
--- | Get a summary of a list of values.
-summary :: [Bool] -> Summary
-summary = foldl' update empty
+-- | Proportion of 'True' values.
+mean :: Summary -> Double
+mean (S n x) = fromIntegral x / fromIntegral n
 
+-- | Standard error for the mean (proportion of 'True' values).
+meanSE :: Summary -> Double
+meanSE s = sqrt (p*(1-p) / n)
+  where p = mean s
+        n = fromIntegral $ size s
+
+-- | Central Limit Theorem based confidence interval for the
+-- population mean (proportion) at the specified coverage level.  The
+-- level must be in the range @(0,1)@.
+meanCI :: Double -> Summary -> (Double,Double)
+meanCI level s = interval level (mean s) (meanSE s)
+
 -- | Get an empty summary.
 empty :: Summary
 empty = S 0 0
 
--- | Take the union of two summaries.
-union :: Summary -> Summary -> Summary
-union (S na ca) (S nb cb) = S (na + nb) (ca + cb)
+-- | Summarize a single value.
+singleton :: Bool -> Summary
+singleton x = S 1 (if x then 1 else 0)
 
 -- | Update the summary with a data point.
-update :: Summary -> Bool -> Summary
-update (S n c) i =
+insert :: Bool -> Summary -> Summary
+insert y (S n x) =
     let n' = n+1
-        c' = if i then c+1 else c
-    in S n' c'
+        x' = if y then x+1 else x
+    in S n' x'
 
--- | Get the sample size.
-sampleSize :: Summary -> Int
-sampleSize (S n _) = n
+-- | Apply a function and update the summary with the result.
+insertWith :: (a -> Bool) -> a -> Summary -> Summary
+insertWith f a = insert (f a)
 
--- | Get the number of 'True' values
-count :: Summary -> Int
-count (S _ c) = c
+-- | Take the union of two summaries.
+union :: Summary -> Summary -> Summary
+union (S na xa) (S nb xb) = S (na + nb) (xa + xb)
 
--- | Get the proportion of 'True' events.
-sampleMean :: Summary -> Double
-sampleMean (S n c) = fromIntegral c / fromIntegral n
+-- | Take the union of a list of summaries.
+unions :: [Summary] -> Summary
+unions = foldl' union empty
 
--- | Get the standard error for the sample proportion.
-sampleSE :: Summary -> Double
-sampleSE s = sqrt (p*(1-p) / n)
-  where p = sampleMean s
-        n = fromIntegral $ sampleSize s
+-- | Get a summary of a list of values.
+fromList :: [Bool] -> Summary
+fromList = foldl' (flip insert) empty
 
--- | Get a Central Limit Theorem-based confidence interval for the mean
--- with the specified coverage level.  The level must be in the range @(0,1)@.
-sampleCI :: Double -> Summary -> (Double,Double)
-sampleCI level s = interval level (sampleMean s) (sampleSE s)
+-- | Map a function over a list of values and summarize the results.
+fromListWith :: (a -> Bool) -> [a] -> Summary
+fromListWith f = fromList . map f
+
+-- | Convert to (size,sum).
+toStats :: Summary -> (Int,Int)
+toStats (S n x) = (n,x)
+
+-- | Convert from (size,sum).  No validation is performed.
+fromStats :: Int -> Int -> Summary
+fromStats = S
diff --git a/lib/Data/Summary/Double.hs b/lib/Data/Summary/Double.hs
--- a/lib/Data/Summary/Double.hs
+++ b/lib/Data/Summary/Double.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE DeriveDataTypeable #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module     : Data.Summary.Double
@@ -10,29 +11,57 @@
 --
 
 module Data.Summary.Double (
-    -- * The @Summary@ data type
+    -- * Summary type
     Summary,
-    summary,
-    update,
 
-    -- * @Summary@ properties
-    sampleSize,
-    sampleMin,
-    sampleMax,
-    sampleMean,
-    sampleSE,
-    sampleVar,
-    sampleSD,
-    sampleCI,
+    -- * Properties
+    -- ** Sum
+    size,
+    sum,
+    sumSquaredErrors,
+    -- ** Mean
+    mean,
+    meanSE,
+    meanCI,
+    -- ** Scale
+    stddev,
+    variance,
+    -- ** Extremes
+    min,
+    max,
 
+    -- * Construction
+    empty,
+    singleton,
+
+    -- * Insertion
+    insert,
+    insertWith,
+
+    -- * Combination
+    union,
+    unions,
+
+    -- * Conversion
+    -- ** Lists
+    fromList,
+    fromListWith,
+
+    -- ** Statistics
+    toStats,
+    fromStats,
+
     ) where
 
-import Control.DeepSeq
+import Prelude hiding (min, max, sum)
+import qualified Prelude as P
+
+import Data.Data( Data, Typeable )
 import Data.List( foldl' )
 import Data.Monoid
 import Text.Printf
 
-import Data.Summary.Utils
+import Data.Summary.Utils( interval )
 
 
 -- | A type for storing summary statistics for a data set including
@@ -42,38 +71,78 @@
                  {-# UNPACK #-} !Double  -- sum of squares
                  {-# UNPACK #-} !Double  -- sample min
                  {-# UNPACK #-} !Double  -- sample max
+    deriving(Eq, Data, Typeable)
 
 instance Show Summary where
     show s@(S n mu _ l h) =
-        printf "    sample size: %d" n
+             printf "    sample size: %d" n
         ++ printf "\n            min: %g" l
         ++ printf "\n            max: %g" h
         ++ printf "\n           mean: %g" mu
-        ++ printf "\n             SE: %g" (sampleSE s)
+        ++ printf "\n             SE: %g" (meanSE s)
         ++ printf "\n         99%% CI: (%g, %g)" c1 c2
-      where (c1,c2) = sampleCI 0.99 s
+      where (c1,c2) = meanCI 0.99 s
 
 instance Monoid Summary where
     mempty = empty
     mappend = union
 
-instance NFData Summary
+-- | Number of observations.
+size :: Summary -> Int
+size (S n _ _ _ _) = n
 
+-- | Sum of values.
+sum :: Summary -> Double
+sum s = (fromIntegral $ size s) * (mean s)
 
--- | Get a summary of a list of values.
-summary :: [Double] -> Summary
-summary = foldl' update empty
+-- | Mean value.
+mean :: Summary -> Double
+mean (S _ m _ _ _) = m
 
--- | Get an empty summary.
+-- | Standard error of the mean.
+meanSE :: Summary -> Double
+meanSE s = sqrt (variance s / fromIntegral (size s))
+
+-- | Get a Central Limit Theorem based confidence interval for the
+-- population mean with the specified coverage level.  The level must
+-- be in the range @(0,1)@.
+meanCI :: Double -> Summary -> (Double,Double)
+meanCI level s = interval level (mean s) (meanSE s)
+
+-- | Sample standard deviation.
+stddev :: Summary -> Double
+stddev s = sqrt (variance s)
+
+-- | Sample variance.
+variance :: Summary -> Double
+variance s = (sumSquaredErrors s) / fromIntegral (size s - 1)
+
+-- | Sum of squared errors @(x[i] - mean)^2@.
+sumSquaredErrors :: Summary -> Double
+sumSquaredErrors (S _ _ s _ _) = s
+
+-- | Minimum value.
+min :: Summary -> Double
+min (S _ _ _ l _) = l
+
+-- | Maximum value.
+max :: Summary -> Double
+max (S _ _ _ _ h) = h
+
+-- | An empty summary.
 empty :: Summary
 empty = S 0 0 0 (1/0) (-1/0)
 
+-- | Summarize a single value.
+singleton :: Double -> Summary
+singleton x = S 1 x 0 x x
+
 -- | Update the summary with a data point.
 -- Running mean and variance computed as in Knuth, Vol 2, page 232,
--- 3rd edition, see http://www.johndcook.com/standard_deviation.html for
+-- 3rd edition, see <http://www.johndcook.com/standard_deviation.html> for
 -- a description.
-update :: Summary -> Double -> Summary
-update (S n m s l h) x =
+insert :: Double -> Summary -> Summary
+insert x (S n m s l h) =
     let n'    = n+1
         delta = x - m
         m'    = m + delta / fromIntegral n'
@@ -82,10 +151,14 @@
         h'    = if x > h then x else h
     in S n' m' s' l' h'
 
+-- | Apply a function and update the summary with the result.
+insertWith :: (a -> Double) -> a -> Summary -> Summary
+insertWith f a = insert (f a)
+
 -- | Take the union of two summaries.
--- Use the updating rules from Chan et al. "Updating Formulae and a Pairwise
---   Algorithm for Computing Sample Variances," available at
--- ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/79/773/CS-TR-79-773.pdf
+-- Use the updating rules from Chan et al. \"Updating Formulae and a Pairwise
+--   Algorithm for Computing Sample Variances,\" available at
+--   <http://infolab.stanford.edu/pub/cstr/reports/cs/tr/79/773/CS-TR-79-773.pdf>.
 union :: Summary -> Summary -> Summary
 union (S na ma sa la ha) (S nb mb sb lb hb) =
     let delta = mb - ma
@@ -97,39 +170,28 @@
            | otherwise = ma + weightedDelta
         s  | n == 0    = 0
            | otherwise = sa + sb + delta*na'*weightedDelta
-        l  = min la lb
-        h  = max ha hb
+        l  = P.min la lb
+        h  = P.max ha hb
     in S n m s l h
 
--- | Get the sample size.
-sampleSize :: Summary -> Int
-sampleSize (S n _ _ _ _) = n
-
--- | Get the sample mean.
-sampleMean :: Summary -> Double
-sampleMean (S _ m _ _ _) = m
-
--- | Get the sample variance.
-sampleVar :: Summary -> Double
-sampleVar (S n _ s _ _) = s / fromIntegral (n - 1)
+-- | Take the union of a list of summaries.
+unions :: [Summary] -> Summary
+unions = foldl' union empty
 
--- | Get the sample standard deviation.
-sampleSD :: Summary -> Double
-sampleSD s = sqrt (sampleVar s)
+-- | Get a summary of a list of values.
+fromList :: [Double] -> Summary
+fromList = foldl' (flip insert) empty
 
--- | Get the sample standard error.
-sampleSE :: Summary -> Double
-sampleSE s = sqrt (sampleVar s / fromIntegral (sampleSize s))
+-- | Map a function over a list of values and summarize the results.
+fromListWith :: (a -> Double) -> [a] -> Summary
+fromListWith f = fromList . map f
 
--- | Get a Central Limit Theorem-based confidence interval for the mean
--- with the specified coverage level.  The level must be in the range @(0,1)@.
-sampleCI :: Double -> Summary -> (Double,Double)
-sampleCI level s = interval level (sampleMean s) (sampleSE s)
+-- | Convert to (size, mean, sumSquaredErrors, min, max).
+toStats :: Summary -> (Int,Double,Double,Double,Double)
+toStats (S n m s l u) = (n,m,s,l,u)
 
--- | Get the minimum of the sample.
-sampleMin :: Summary -> Double
-sampleMin (S _ _ _ l _) = l
+-- | Convert from (size, mean, sumSquaredErrors, min, max).
+-- No validation is performed.
+fromStats :: Int -> Double -> Double -> Double -> Double -> Summary
+fromStats = S
 
--- | Get the maximum of the sample.
-sampleMax :: Summary -> Double
-sampleMax (S _ _ _ _ h) = h
diff --git a/lib/Data/Summary/Utils.hs b/lib/Data/Summary/Utils.hs
--- a/lib/Data/Summary/Utils.hs
+++ b/lib/Data/Summary/Utils.hs
@@ -16,12 +16,13 @@
 
 import GSL.Random.Dist( ugaussianPInv )
 
--- | Get a Central Limit Theorem-based confidence interval for the
--- population mean with the specified coverage level.  The level must
--- be in the range @(0,1)@.
-interval :: Double -- ^ the confidence level
-         -> Double -- ^ the sample mean
-         -> Double -- ^ the sample standard error
+
+-- | Get a Central Limit Theorem based confidence interval for a
+-- population parameter with the specified coverage level.  The level
+-- must be in the range @(0,1)@.
+interval :: Double -- ^ Confidence level.
+         -> Double -- ^ Estimate.
+         -> Double -- ^ Standard error.
          -> (Double,Double)
 interval level xbar se | not (level > 0 && level < 1) =
                              error "level must be between 0 and 1"
@@ -31,6 +32,7 @@
         delta = z*se
     in (xbar-delta, xbar+delta)
 
--- | Tests if the value is in the open interval (a,b)
+
+-- | Tests if the value is in the open interval (a,b).
 inInterval :: Double -> (Double,Double) -> Bool
 x `inInterval` (a,b) = x > a && x < b
diff --git a/monte-carlo.cabal b/monte-carlo.cabal
--- a/monte-carlo.cabal
+++ b/monte-carlo.cabal
@@ -1,5 +1,5 @@
 name:           monte-carlo
-version:        0.5
+version:        0.6
 license:        BSD3
 license-file:   LICENSE
 author:         Patrick Perry
@@ -7,47 +7,45 @@
 homepage:       http://github.com/patperry/hs-monte-carlo
 category:       Math
 synopsis:       A monad and transformer for Monte Carlo calculations.
-description:    A monad and transformer for Monte Carlo calculations.  The
-                monads carry and provide access to a random number generator.
-                Importantly, they only keep one copy of the generator state,
-                and so are much more efficient than MonadRandom.  Currently,
-                only the generator that comes with the GNU Scientific Library
-                (GSL) is supported.
+description:    A monad and transformer for performing Monte Carlo
+                calculations. This monad carries and provides access to
+                a pseudo-random number generator. Internally, the monad
+                mutates rather than copies the random gnerator state.  By
+                avoiding copies, it can deliver faster performance than
+                many pure random number implementations.  The package is
+                built around the facilities provided by the GNU Scientific
+                Library (GSL).
 build-type:     Simple
 stability:      experimental
 cabal-version:  >= 1.8
 extra-source-files: NEWS examples/Binomial.hs examples/Pi.lhs
-                    examples/Poker.hs  examples/Queue.hs tests/Main.hs
+                    examples/Poker.hs examples/Queue.hs tests/Main.hs
                     tests/Makefile
 
 library
     exposed-modules:
             Control.Monad.MC
-            Control.Monad.MC.Class
-            Control.Monad.MC.GSL
-            Data.Summary
             Data.Summary.Bool
             Data.Summary.Double
-            Data.Summary.Utils
 
     other-modules:
-            Control.Monad.MC.Base
+            Control.Monad.MC.GSL
             Control.Monad.MC.GSLBase
             Control.Monad.MC.Repeat
             Control.Monad.MC.Sample
             Control.Monad.MC.Walker
+            Data.Summary.Utils
 
     extensions:
-            FlexibleInstances,
-            MultiParamTypeClasses,
-            TypeFamilies,
-            UndecidableInstances
+            DeriveDataTypeable
+            RankNTypes
+            TypeFamilies
 
-    build-depends:  base       >= 4     && < 5,
-                    gsl-random >= 0.4.3 && < 0.5,
-                    mtl        >= 1.1   && < 3.0,
-                    vector     >= 0.6   && < 0.11,
-                    deepseq    >= 1.0   && < 2.0
+    build-depends:  base         >= 4     && < 5,
+                    gsl-random   >= 0.5   && < 0.6,
+                    primitive    >= 0.5   && < 0.6,
+                    transformers >= 0.3   && < 0.5,
+                    vector       >= 0.6   && < 0.11
 
     hs-source-dirs: lib
     ghc-options:    -Wall
@@ -66,13 +64,14 @@
   ghc-options:
       -Wall -Werror
   build-depends:
-      base       >= 4     && < 5
-    , gsl-random >= 0.4.3 && < 0.5
-    , mtl        >= 1.1   && < 3.0
-    , vector     >= 0.6   && < 0.11
-    , ieee754 >= 0.7 && < 0.8
-    , random >=1.0 && < 1.1
-    , QuickCheck >= 2.4.0.1 && < 2.6
-    , test-framework-quickcheck2 >= 0.2 && < 0.4
-    , test-framework >= 0.4 && < 0.9
-    , deepseq >= 1.0 && < 2.0
+      base                       >= 4       && < 5,
+      gsl-random                 >= 0.5     && < 0.6,
+      primitive                  >= 0.5     && < 0.6,
+      transformers               >= 0.3     && < 0.5,
+      vector                     >= 0.6     && < 0.11,
+      ieee754                    >= 0.7     && < 0.8,
+      random                     >= 1.0     && < 1.1,
+      QuickCheck                 >= 2.4.0.1 && < 2.8,
+      test-framework             >= 0.4     && < 0.9,
+      test-framework-quickcheck2 >= 0.2     && < 0.4
+
diff --git a/tests/Main.hs b/tests/Main.hs
--- a/tests/Main.hs
+++ b/tests/Main.hs
@@ -4,7 +4,8 @@
 import Control.Monad
 import Data.AEq
 import Data.Monoid
-import Data.Summary
+import Data.Summary.Double( Summary )
+import qualified Data.Summary.Double as S
 import Test.QuickCheck
 import Test.Framework
 import Test.Framework.Providers.QuickCheck2
@@ -40,17 +41,17 @@
 
 prop_monoid_update_equiv :: [Double] -> [Double] -> Bool
 prop_monoid_update_equiv xs ys =
-    approxEqualS (summary $ xs <> ys)
-                 (summary xs <> summary ys)
+    approxEqualS (S.fromList $ xs <> ys)
+                 (S.fromList xs <> S.fromList ys)
 
 prop_monoid_assoc :: [Double] -> [Double] -> [Double] -> Bool
 prop_monoid_assoc xs ys zs =
-    let (sxs, sys, szs) = (summary xs, summary ys, summary zs)
+    let (sxs, sys, szs) = (S.fromList xs, S.fromList ys, S.fromList zs)
      in ((sxs <> sys) <> szs) `approxEqualS` (sxs <> (sys <> szs))
 
 prop_monoid_commute :: [Double] -> [Double] -> Bool
 prop_monoid_commute xs ys =
-    let (sxs, sys) = (summary xs, summary ys)
+    let (sxs, sys) = (S.fromList xs, S.fromList ys)
      in (sxs <> sys) `approxEqualS` (sys <> sxs)
 
 tests_monoid :: Test
@@ -70,8 +71,8 @@
 
 approxEqualS :: Summary -> Summary -> Bool
 approxEqualS a b =
-    sampleSize a == sampleSize b &&
-      all eq [ sampleMin, sampleMax, sampleMean, sampleVar ]
+    S.size a == S.size b &&
+      all eq [ S.min, S.max, S.mean, S.variance ]
     where
         eq f = f a ~== f b