diff --git a/random-fu.cabal b/random-fu.cabal
--- a/random-fu.cabal
+++ b/random-fu.cabal
@@ -1,6 +1,6 @@
 name:                   random-fu
-version:                0.0.3.2
-stability:              experimental
+version:                0.1.0.0
+stability:              provisional
 
 cabal-version:          >= 1.2
 build-type:             Simple
@@ -12,34 +12,65 @@
 
 category:               Math
 synopsis:               Random number generation
-description:            Random number generation based on orthogonal typeclasses
-                        for entropy sources and random variable distributions, all
-                        served up on a monadic platter.  Aspires to be useful
-                        in an idiomatic way in both \"pure\" and \"impure\" styles,
-                        as well as reasonably fast.
+description:            Random number generation based on modeling random 
+                        variables in two complementary ways: first, by the
+                        parameters of standard mathematical distributions and,
+                        second, by an abstract type ('RVar') which can be
+                        composed and manipulated monadically and sampled in
+                        either monadic or \"pure\" styles.
+                        
+                        The primary purpose of this library is to support 
+                        defining and sampling a wide variety of high quality
+                        random variables.  Quality is prioritized over speed,
+                        but performance is an important goal too.
+                        
+                        In my testing, I have found it capable of speed 
+                        comparable to other Haskell libraries, but still
+                        a fair bit slower than straight C implementations of 
+                        the same algorithms.
+                        
+                        Warning to anyone upgrading from \"< 0.1\": 'Discrete'
+                        has been renamed 'Categorical', the entropy source 
+                        classes have been redesigned, and many things are no
+                        longer exported from the root module "Data.Random"
+                        (In particular, DevRandom - this is not available on 
+                        windows, so it will likely move to its own package 
+                        eventually so that client code dependencies on it will 
+                        be made explicit).
+                        
+                        The "Data.Random" module itself should now have a
+                        relatively stable interface, but the other modules
+                        are still subject to change.  Specifically, I am 
+                        considering hiding data constructors for most or all 
+                        of the distributions.
 
 Flag base4
 Flag base4_2
     Description:        base-4.2 has an incompatible change in Data.Fixed (HasResolution)
 
 Library
+  ghc-options:          -Wall -funbox-strict-fields -fno-method-sharing
   hs-source-dirs:       src
   exposed-modules:      Data.Random
                         Data.Random.Distribution
                         Data.Random.Distribution.Bernoulli
                         Data.Random.Distribution.Beta
                         Data.Random.Distribution.Binomial
-                        Data.Random.Distribution.Discrete
+                        Data.Random.Distribution.Categorical
+                        Data.Random.Distribution.Dirichlet
                         Data.Random.Distribution.Exponential
                         Data.Random.Distribution.Gamma
+                        Data.Random.Distribution.Multinomial
                         Data.Random.Distribution.Normal
                         Data.Random.Distribution.Poisson
                         Data.Random.Distribution.Rayleigh
                         Data.Random.Distribution.Triangular
                         Data.Random.Distribution.Uniform
+                        Data.Random.Distribution.Weibull
                         Data.Random.Distribution.Ziggurat
                         Data.Random.Internal.Find
                         Data.Random.Internal.Fixed
+                        Data.Random.Internal.Primitives
                         Data.Random.Internal.TH
                         Data.Random.Internal.Words
                         Data.Random.Lift
@@ -47,7 +78,7 @@
                         Data.Random.RVar
                         Data.Random.Sample
                         Data.Random.Source
-                        Data.Random.Source.DevRandom
+                        Data.Random.Source.MWC
                         Data.Random.Source.PureMT
                         Data.Random.Source.Std
                         Data.Random.Source.StdGen
@@ -64,12 +95,21 @@
     
   build-depends:        array,
                         containers,
-                        erf,
                         mersenne-random-pure64,
                         monad-loops >= 0.3.0.1,
+                        MonadPrompt,
+                        mwc-random,
                         mtl,
                         random,
                         random-shuffle,
                         stateref >= 0.3 && < 0.4,
-                        storablevector,
-                        template-haskell
+                        tagged,
+                        template-haskell,
+                        vector
+
+  if os(Windows)
+    cpp-options:        -Dwindows
+    build-depends:      erf-native
+  else
+    build-depends:      erf
+    exposed-modules:    Data.Random.Source.DevRandom
diff --git a/src/Data/Random.hs b/src/Data/Random.hs
--- a/src/Data/Random.hs
+++ b/src/Data/Random.hs
@@ -1,66 +1,78 @@
-{-
- -      ``Data/Random''
- -}
-{-# LANGUAGE
-    FlexibleContexts
-  #-}
-
--- |Random numbers and stuff...
+-- |Flexible modeling and sampling of random variables.
+--
+-- The central abstraction in this library is the concept of a random 
+-- variable.  It is not fully formalized in the standard measure-theoretic 
+-- language, but rather is informally defined as a \"thing you can get random 
+-- values out of\".  Different random variables may have different types of 
+-- values they can return or the same types but different probabilities for
+-- each value they can return.  The random values you get out of them are
+-- traditionally called \"random variates\".
 -- 
--- "Data.Random.Source" exports the typeclasses for entropy sources, and
--- Data.Random.Source.* export various instances and/or functions with which
--- instances can be defined.
+-- Most imperative-language random number libraries are all about obtaining 
+-- and manipulating random variates.  This one is about defining, manipulating 
+-- and sampling random variables.  Computationally, the distinction is small 
+-- and mostly just a matter of perspective, but from a program design 
+-- perspective it provides both a powerfully composable abstraction and a
+-- very useful separation of concerns.
 -- 
--- "Data.Random.Distribution" exports the typeclasses for sampling distributions,
--- and Data.Random.Distribution.* export various specific distributions.
---
--- "Data.Random.RVar" exports the 'RVar' type, which is a probability distribution
--- monad that allows for concise definitions of random variables, as well as
--- a couple handy 'RVar's.
-
+-- Abstract random variables as implemented by 'RVar' are composable.  They can
+-- be defined in a monadic / \"imperative\" style that amounts to manipulating
+-- variates, but with strict type-level isolation.  Concrete random variables
+-- are also provided, but they do not compose as generically.  The 'Distribution'
+-- type class allows concrete random variables to \"forget\" their concreteness 
+-- so that they can be composed.  For examples of both, see the documentation 
+-- for 'RVar' and 'Distribution', as well as the code for any of the concrete 
+-- distributions such as 'Uniform', 'Gamma', etc.
+-- 
+-- Both abstract and concrete random variables can be sampled (despite the
+-- types GHCi may list for the functions) by the functions in "Data.Random.Sample".
+-- 
+-- Random variable sampling is done with regard to a generic basis of primitive
+-- random variables defined in "Data.Random.Internal.Primitives".  This basis 
+-- is very low-level and the actual set of primitives is still fairly experimental,
+-- which is why it is in the \"Internal\" sub-heirarchy.  User-defined variables
+-- should use the existing high-level variables such as 'Uniform' and 'Normal'
+-- rather than these basis variables.  "Data.Random.Source" defines classes for
+-- entropy sources that provide implementations of these primitive variables. 
+-- Several implementations are available in the Data.Random.Source.* modules.
 module Data.Random
-    ( module Data.Random.Sample
-    , module Data.Random.Source
-    , module Data.Random.Source.DevRandom
-    , module Data.Random.Source.StdGen
-    , module Data.Random.Source.PureMT
-    , module Data.Random.Source.Std
-    , module Data.Random.Distribution
-    , module Data.Random.Distribution.Bernoulli
-    , module Data.Random.Distribution.Beta
-    , module Data.Random.Distribution.Binomial
-    , module Data.Random.Distribution.Discrete
-    , module Data.Random.Distribution.Gamma
-    , module Data.Random.Distribution.Exponential
-    , module Data.Random.Distribution.Normal
-    , module Data.Random.Distribution.Poisson
-    , module Data.Random.Distribution.Rayleigh
-    , module Data.Random.Distribution.Triangular
-    , module Data.Random.Distribution.Uniform
-    , module Data.Random.Distribution.Ziggurat
-    , module Data.Random.List
-    , module Data.Random.RVar
+    ( -- * Random variables
+      -- ** Abstract ('RVar')
+      RVar, RVarT,
+      runRVar, runRVarT, runRVarTWith,
+
+      -- ** Concrete ('Distribution')
+      Distribution(..), CDF(..),
+      
+      -- * Sampling random variables
+      Sampleable(..), sample, sampleState, sampleStateT,
+      
+      -- * A few very common distributions
+      Uniform(..), uniform, 
+      StdUniform(..), stdUniform,
+      Normal(..), normal, stdNormal,
+      Gamma(..), gamma,
+      
+      -- * Entropy Sources
+      MonadRandom, RandomSource, StdRandom(..),
+      
+      -- * Useful list-based operations
+      randomElement,
+      shuffle, shuffleN, shuffleNofM
+      
     ) where
 
 import Data.Random.Sample
-import Data.Random.Source
-import Data.Random.Source.DevRandom
-import Data.Random.Source.StdGen
-import Data.Random.Source.PureMT
+import Data.Random.Source (MonadRandom, RandomSource)
+import Data.Random.Source.MWC ()
+import Data.Random.Source.StdGen ()
+import Data.Random.Source.PureMT ()
 import Data.Random.Source.Std
 import Data.Random.Distribution
-import Data.Random.Distribution.Bernoulli
-import Data.Random.Distribution.Beta
-import Data.Random.Distribution.Binomial
-import Data.Random.Distribution.Discrete
 import Data.Random.Distribution.Gamma
-import Data.Random.Distribution.Exponential
 import Data.Random.Distribution.Normal
-import Data.Random.Distribution.Poisson
-import Data.Random.Distribution.Rayleigh
-import Data.Random.Distribution.Triangular
 import Data.Random.Distribution.Uniform
-import Data.Random.Distribution.Ziggurat
+
 import Data.Random.Lift ()
 import Data.Random.List
 import Data.Random.RVar
diff --git a/src/Data/Random/Distribution.hs b/src/Data/Random/Distribution.hs
--- a/src/Data/Random/Distribution.hs
+++ b/src/Data/Random/Distribution.hs
@@ -1,22 +1,66 @@
-{-
- -      ``Data/Random/Distribution''
- -}
 {-# LANGUAGE
     MultiParamTypeClasses, FlexibleContexts
   #-}
-
 module Data.Random.Distribution where
 
 import Data.Random.Lift
 import Data.Random.RVar
 
--- |A definition of a random variable's distribution.  From the distribution
--- an 'RVar' can be created, or the distribution can be directly sampled using 
--- 'sampleFrom' or 'sample'.
+-- |A 'Distribution' is a data representation of a random variable's probability
+-- structure.  For example, in "Data.Random.Distribution.Normal", the 'Normal'
+-- distribution is defined as:
+--
+-- > data Normal a
+-- >     = StdNormal
+-- >     | Normal a a
+-- 
+-- Where the two parameters of the 'Normal' data constructor are the mean and
+-- standard deviation of the random variable, respectively.  To make use of
+-- the 'Normal' type, one can convert it to an 'rvar' and manipulate it or
+-- sample it directly:
+--
+-- > x <- sample (rvar (Normal 10 2))
+-- > x <- sample (Normal 10 2)
+-- 
+-- A 'Distribution' is typically more transparent than an 'RVar'
+-- but less composable (precisely because of that transparency).  There are 
+-- several practical uses for types implementing 'Distribution':
+-- 
+-- * Typically, a 'Distribution' will expose several parameters of a standard 
+-- mathematical model of a probability distribution, such as mean and std deviation for
+-- the normal distribution.  Thus, they can be manipulated analytically using
+-- mathematical insights about the distributions they represent.  For example,
+-- a collection of bernoulli variables could be simplified into a (hopefully) smaller
+-- collection of binomial variables.
+-- 
+-- * Because they are generally just containers for parameters, they can be
+-- easily serialized to persistent storage or read from user-supplied 
+-- configurations (eg, initialization data for a simulation).
+-- 
+-- * If a type additionally implements the 'CDF' subclass, which extends
+-- 'Distribution' with a cumulative density function, an arbitrary random
+-- variable 'x' can be tested against the distribution by testing
+-- @fmap (cdf dist) x@ for uniformity.
+-- 
+-- On the other hand, most 'Distribution's will not be closed under all the
+-- same operations as 'RVar' (which, being a monad, has a fully turing-complete
+-- internal computational model).  The sum of two uniformly-distributed 
+-- variables, for example, is not uniformly distributed.  To support general 
+-- composition, the 'Distribution' class defines a function 'rvar' to 
+-- construct the more-abstract and more-composable 'RVar' representation 
+-- of a random variable.
 class Distribution d t where
     -- |Return a random variable with this distribution.
     rvar :: d t -> RVar t
+    rvar = rvarT
+    
+    -- |Return a random variable with the given distribution, pre-lifted to an arbitrary 'RVarT'.
+    -- Any arbitrary 'RVar' can also be converted to an 'RVarT m' for an arbitrary 'm', using
+    -- either 'lift' or 'sample'.
+    rvarT :: d t -> RVarT n t
+    rvarT d = lift (rvar d)
 
+
 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
@@ -31,16 +75,15 @@
     -- must be uniformly distributed over (0,1).  Inclusion of either endpoint is optional,
     -- though the preferred range is (0,1].
     -- 
-    -- Thus, 'cdf' for a product type should not be a joint CDF as commonly 
-    -- defined, as that definition violates both conditions.
+    -- Note that this definition requires that  'cdf' for a product type 
+    -- should _not_ be a joint CDF as commonly defined, as that definition 
+    -- violates both conditions.
     -- Instead, it should be a univariate CDF over the product type.  That is,
     -- it should represent the CDF with respect to the lexicographic order
-    -- of the tuple.
+    -- of the product.
+    -- 
+    -- The present specification is probably only really useful for testing
+    -- conformance of a variable to its target distribution, and I am open to
+    -- suggestions for more-useful specifications (especially with regard to
+    -- the interaction with product types).
     cdf :: d t -> t -> Double
-
--- |Return a random variable with the given distribution, pre-lifted to an arbitrary 'RVarT'.
--- Any arbitrary 'RVar' can also be converted to an 'RVarT m' for an arbitrary 'm', using
--- either 'lift' or 'sample'.
-rvarT :: Distribution d t => d t -> RVarT n t
-rvarT d = lift (rvar d)
-
diff --git a/src/Data/Random/Distribution/Bernoulli.hs b/src/Data/Random/Distribution/Bernoulli.hs
--- a/src/Data/Random/Distribution/Bernoulli.hs
+++ b/src/Data/Random/Distribution/Bernoulli.hs
@@ -16,8 +16,6 @@
 import Data.Random.Distribution
 import Data.Random.Distribution.Uniform
 
-import Data.Int
-import Data.Word
 import Data.Ratio
 import Data.Complex
 
@@ -35,7 +33,7 @@
     return (x <= p)
 
 boolBernoulliCDF :: (Real a) => a -> Bool -> Double
-boolBernoulliCDF p True  = 1
+boolBernoulliCDF _ True  = 1
 boolBernoulliCDF p False = (1 - realToFrac p)
 
 -- | @generalBernoulli t f p@ generates a random variable whose value is @t@
@@ -46,10 +44,10 @@
     return (if x then t else f)
 
 generalBernoulliCDF :: CDF (Bernoulli b) Bool => (a -> a -> Bool) -> a -> a -> b -> a -> Double
-generalBernoulliCDF (>=) f t p x
-    | f >= t    = error "generalBernoulliCDF: f >= t"
-    | x >= t    = cdf (Bernoulli p) True
-    | x >= f    = cdf (Bernoulli p) False
+generalBernoulliCDF gte f t p x
+    | f `gte` t = error "generalBernoulliCDF: f >= t"
+    | x `gte` t = cdf (Bernoulli p) True
+    | x `gte` f = cdf (Bernoulli p) False
     | otherwise = 0
 
 data Bernoulli b a = Bernoulli b
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
@@ -4,18 +4,21 @@
 {-# LANGUAGE
     MultiParamTypeClasses,
     FlexibleInstances, FlexibleContexts,
-    UndecidableInstances
+    UndecidableInstances,
+    TemplateHaskell
   #-}
 
 module Data.Random.Distribution.Beta where
 
+import Data.Random.Internal.TH
+
 import Data.Random.RVar
 import Data.Random.Distribution
 import Data.Random.Distribution.Gamma
 import Data.Random.Distribution.Uniform
 
-import Control.Monad
-
+{-# SPECIALIZE fractionalBeta :: Float  -> Float  -> RVar Float #-}
+{-# SPECIALIZE fractionalBeta :: Double -> Double -> RVar Double #-}
 fractionalBeta :: (Fractional a, Distribution Gamma a, Distribution StdUniform a) => a -> a -> RVar a
 fractionalBeta 1 1 = stdUniform
 fractionalBeta a b = do
@@ -23,16 +26,14 @@
     y <- gamma b 1
     return (x / (x + y))
 
-fractionalBetaFromIntegral :: (Fractional c, Distribution (Erlang a) c, Distribution (Erlang b) c) => a -> b -> RVar c
-fractionalBetaFromIntegral a b =  do
-    x <- erlang a
-    y <- erlang b
-    return (x / (x + y))
-
+{-# SPECIALIZE beta :: Float  -> Float  -> RVar Float #-}
+{-# SPECIALIZE beta :: Double -> Double -> RVar Double #-}
 beta :: Distribution Beta a => a -> a -> RVar a
 beta a b = rvar (Beta a b)
 
 data Beta a = Beta a a
 
-instance (Fractional a, Distribution Gamma a, Distribution StdUniform a) => Distribution Beta a where
-    rvar (Beta a b) = fractionalBeta a b
+$( replicateInstances ''Float realFloatTypes [d|
+        instance Distribution Beta Float
+              where rvar (Beta a b) = fractionalBeta a b
+    |])
diff --git a/src/Data/Random/Distribution/Binomial.hs b/src/Data/Random/Distribution/Binomial.hs
--- a/src/Data/Random/Distribution/Binomial.hs
+++ b/src/Data/Random/Distribution/Binomial.hs
@@ -4,7 +4,8 @@
 {-# LANGUAGE
     MultiParamTypeClasses,
     FlexibleInstances, FlexibleContexts,
-    UndecidableInstances, TemplateHaskell
+    UndecidableInstances, TemplateHaskell,
+    BangPatterns
   #-}
 
 module Data.Random.Distribution.Binomial where
@@ -16,29 +17,19 @@
 import Data.Random.Distribution.Beta
 import Data.Random.Distribution.Uniform
 
-import Data.Int
-import Data.Word
-import Control.Monad
-
     -- algorithm from Knuth's TAOCP, 3rd ed., p 136
     -- specific choice of cutoff size taken from gsl source
     -- note that although it's fast enough for large (eg, 2^10000) 
     -- @Integer@s, it's not accurate enough when using @Double@ as
     -- the @b@ parameter.
+{-# SPECIALIZE integralBinomial :: Int -> Float  -> RVar Int #-}
+{-# SPECIALIZE integralBinomial :: Int -> Double -> RVar Int #-}
+{-# SPECIALIZE integralBinomial :: Integer -> Float  -> RVar Integer #-}
+{-# SPECIALIZE integralBinomial :: Integer -> Double -> RVar Integer #-}
 integralBinomial :: (Integral a, Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => a -> b -> RVar a
-integralBinomial t p = bin 0 t p
+integralBinomial = bin 0
     where
-        -- GHC likes to discharge Beta to the Beta instance's context, which
-        -- @integralBinomial@'s context doesn't (directly) satisfy.
-        -- Seems like GHC could do better.  GHC could discharge it in @integralBinomial@'s
-        -- context when attempting to satisfy bin, since the Beta instance covers all @b@,
-        -- whereupon it would find that the contexts do, in fact, match.
-        -- Of course, it's a pretty obscure case, so maybe it's better to just 
-        -- leave it to the coder who's doing weird stuff to tell the compiler what
-        -- he/she wants.
-        -- Anyway, this type signature makes GHC happy.
-        bin :: (Integral a, Floating b, Ord b, Distribution Beta b, Distribution StdUniform b) => a -> a -> b -> RVar a
-        bin k t p
+        bin !k !t !p
             | t > 10    = do
                 let a = 1 + t `div` 2
                     b = 1 + t - a
@@ -50,24 +41,16 @@
         
             | otherwise = count k t
                 where
-                    count k  0    = return k
-                    count k (n+1) = do
+                    count !k'  0    = return k'
+                    count !k' (n+1) = do
                         x <- stdUniform
-                        (count $! (if x < p then k + 1 else k)) n
-        
-        -- this gives 'believable' results... but is it really any good?
-        selectA cutoff t = return a
---             | cutoff >= a   = return a
---             | otherwise     = generalBernoulli a (a `div` cutoff) p
---             
-             where
-                 a = 1 + t `div` 2
---                 p = 0.5 :: Float -- fromIntegral cutoff / fromIntegral a :: Float
-                
+                        count (if x < p then k' + 1 else k') n
+                    count _ _ = error "integralBinomial: negative number of trials specified"
+
 -- TODO: improve performance
 integralBinomialCDF :: (Integral a, Real b) => a -> b -> a -> Double
-integralBinomialCDF n p x = sum
-    [ fromIntegral (n `c` i) * p' ^^ i * (1-p') ^^ (n-i)
+integralBinomialCDF t p x = sum
+    [ fromIntegral (t `c` i) * p' ^^ i * (1-p') ^^ (t-i)
     | i <- [0 .. x]
     ]
     
@@ -77,12 +60,24 @@
 
 -- would it be valid to repeat the above computation using fractional @t@?
 -- obviously something different would have to be done with @count@ as well...
+{-# SPECIALIZE floatingBinomial :: Float  -> Float  -> RVar Float  #-}
+{-# SPECIALIZE floatingBinomial :: Float  -> Double -> RVar Float  #-}
+{-# SPECIALIZE floatingBinomial :: Double -> Float  -> RVar Double #-}
+{-# SPECIALIZE floatingBinomial :: Double -> Double -> RVar Double #-}
 floatingBinomial :: (RealFrac a, Distribution (Binomial b) Integer) => a -> b -> RVar a
 floatingBinomial t p = fmap fromInteger (rvar (Binomial (truncate t) p))
 
 floatingBinomialCDF :: (CDF (Binomial b) Integer, RealFrac a) => a -> b -> a -> Double
 floatingBinomialCDF t p x = cdf (Binomial (truncate t :: Integer) p) (floor x)
 
+{-# SPECIALIZE binomial :: Int -> Float  -> RVar Int #-}
+{-# SPECIALIZE binomial :: Int -> Double -> RVar Int #-}
+{-# SPECIALIZE binomial :: Integer -> Float  -> RVar Integer #-}
+{-# SPECIALIZE binomial :: Integer -> Double -> RVar Integer #-}
+{-# SPECIALIZE binomial :: Float  -> Float  -> RVar Float  #-}
+{-# SPECIALIZE binomial :: Float  -> Double -> RVar Float  #-}
+{-# SPECIALIZE binomial :: Double -> Float  -> RVar Double #-}
+{-# SPECIALIZE binomial :: Double -> Double -> RVar Double #-}
 binomial :: Distribution (Binomial b) a => a -> b -> RVar a
 binomial t p = rvar (Binomial t p)
 
diff --git a/src/Data/Random/Distribution/Categorical.hs b/src/Data/Random/Distribution/Categorical.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution/Categorical.hs
@@ -0,0 +1,168 @@
+{-
+ -      ``Data/Random/Distribution/Categorical''
+ -}
+{-# LANGUAGE
+    MultiParamTypeClasses,
+    FlexibleInstances, FlexibleContexts
+  #-}
+
+module Data.Random.Distribution.Categorical where
+
+import Data.Random.RVar
+import Data.Random.Distribution
+import Data.Random.Distribution.Uniform
+
+import Control.Arrow
+import Control.Monad
+import Control.Applicative
+import Data.Foldable (Foldable(foldMap))
+import Data.Traversable (Traversable(traverse, sequenceA))
+
+import Data.List
+import Data.Function
+
+-- |Construct a 'Categorical' distribution from a list of probabilities
+-- and categories, where the probabilities all sum to 1.
+categorical :: Distribution (Categorical p) a => [(p,a)] -> RVar a
+categorical ps = rvar (Categorical ps)
+
+-- | Construct a 'Categorical' distribution from a list of weighted categories,
+-- where the weights do not necessarily sum to 1.
+{-# INLINE weightedCategorical #-}
+weightedCategorical :: (Fractional p) => [(p,a)] -> Categorical p a
+weightedCategorical = normalizeCategoricalPs . Categorical
+
+-- |Construct a 'Categorical' distribution from a list of observed outcomes.
+-- Equivalent events will be grouped and counted, and the probabilities of each
+-- event in the returned distribution will be proportional to the number of 
+-- occurrences of that event.
+empirical :: (Fractional p, Ord a) => [a] -> Categorical p a
+empirical xs = normalizeCategoricalPs (Categorical bins)
+    where bins = [ (genericLength bin, x)
+                 | bin@(x:_) <- group (sort xs)
+                 ]
+
+-- |Categorical distribution; a list of events with corresponding probabilities.
+-- The sum of the probabilities must be 1, and no event should have a zero 
+-- or negative probability (at least, at time of sampling; very clever users
+-- can do what they want with the numbers before sampling, just make sure 
+-- that if you're one of those clever ones, you normalize before sampling).
+newtype Categorical p a = Categorical [(p, a)]
+    deriving (Eq, Show)
+
+instance (Fractional p, Ord p, Distribution StdUniform p) => Distribution (Categorical p) a where
+    rvar (Categorical []) = fail "categorical distribution over empty set cannot be sampled"
+    rvar (Categorical ds) = do
+        let (ps, xs) = unzip ds
+            cs = scanl1 (+) ps
+        
+        u <- stdUniform
+        getEvent u cs xs
+        
+        where
+            -- In the (hopefully) extremely rare event that, due to numerical
+            -- instability, the last 'c' is less than 1 _and_ a number greater than 
+            -- it is drawn, simply retry the sampling.  If it comes to that, also
+            -- do one last sanity check that lastC > 0, to make sure that there
+            -- is some nonzero chance of termination.
+            getEvent u cs0 xs0 = go 0 cs0 xs0
+                where
+                    go lastC [] _
+                        | lastC > 0 = do {newU <- stdUniform; getEvent newU cs0 xs0}
+                        | otherwise = fail "categorical distribution sampling error: total probablility not greater than zero"
+                    go lastC (c:cs) (x:xs)
+                        | c < lastC = fail "categorical distribution sampling error: negative probability for an event!"
+                        | u > c     = go c cs xs
+                        | c == c    = return x
+                        | otherwise = fail "categorical distribution sampling error: NaN probability"
+                    
+                    go _ _ _ = error "rvar/Categorical: programming error! this case should be impossible!"
+
+instance Functor (Categorical p) where
+    fmap f (Categorical ds) = Categorical [(p, f x) | ~(p, x) <- ds]
+
+instance Foldable (Categorical p) where
+    foldMap f (Categorical ds) = foldMap (f . snd) ds
+
+instance Traversable (Categorical p) where
+    traverse f (Categorical ds) = Categorical <$> traverse (\(p,e) -> (\e' -> (p,e')) <$> f e) ds
+    sequenceA  (Categorical ds) = Categorical <$> traverse (\(p,e) -> (\e' -> (p,e')) <$>   e) ds
+
+instance Fractional p => Monad (Categorical p) where
+    return x = Categorical [(1, x)]
+    
+    -- I'm not entirely sure whether this is a valid form of failure; see next
+    -- set of comments.
+    fail _ = Categorical []
+    
+    -- Should the normalize step be included here, or should normalization
+    -- be assumed?  It seems like there is (at least) 1 valid situation where
+    -- non-normal results would arise:  the distribution being modeled is 
+    -- "conditional" and some event arose that contradicted the assumed 
+    -- condition and thus was eliminated ('f' returned an empty or 
+    -- zero-probability consequent, possibly by 'fail'ing).
+    -- 
+    -- It seems reasonable to continue in such circumstances, but should there
+    -- be any renormalization?  If so, does it make a difference when that 
+    -- renormalization is done?  I'm pretty sure it does, actually.  So, the
+    -- normalization will be omitted here for now, as it's easier for the
+    -- user (who really better know what they mean if they're returning
+    -- non-normalized probability anyway) to normalize explicitly than to
+    -- undo any normalization that was done automatically.
+    (Categorical xs) >>= f = {- normalizeCategoricalPs . -} Categorical $ do
+        (p, x) <- xs
+        
+        let Categorical fx = f x
+        (q, y) <- fx
+        
+        return (p * q, y)
+
+instance Fractional p => Applicative (Categorical p) where
+    pure = return
+    (<*>) = ap
+
+-- |Like 'fmap', but for the probabilities of a categorical distribution.
+mapCategoricalPs :: (p -> q) -> Categorical p e -> Categorical q e
+mapCategoricalPs f (Categorical ds) = Categorical [(f p, x) | (p, x) <- ds]
+
+-- |Adjust all the weights of a categorical distribution so that they 
+-- sum to unity.
+normalizeCategoricalPs :: (Fractional p) => Categorical p e -> Categorical p e
+normalizeCategoricalPs orig@(Categorical ds) = 
+    -- For practical purposes the scale factor is strict anyway,
+    -- so check if the total probability is 1 and, if so, skip 
+    -- the actual scaling part.
+    --
+    -- Along the way, discard any zero-probability events.
+    if null ds || ps =~ 1
+        then orig
+        else Categorical
+                [ (p * scale, e)
+                | (p, e) <- ds
+                , p /= 0
+                ] 
+    where
+        ps = foldl1' (+) (map fst ds)
+        scale = recip ps
+        
+        -- Using same implicit-epsilon trick as in Distribution instance
+        -- (see comments there)
+        x =~ y  = (100 + (x-y) == 100)
+
+
+-- |Simplify a categorical distribution by combining equivalent categories (the new
+-- category will have a probability equal to the sum of all the originals).
+collectEvents :: (Ord e, Num p, Ord p) => Categorical p e -> Categorical p e
+collectEvents = collectEventsBy compare ((sum *** head) . unzip)
+        
+-- |Simplify a categorical distribution by combining equivalent events (the new
+-- event will have a weight equal to the sum of all the originals).
+-- The comparator function is used to identify events to combine.  Once chosen,
+-- the events and their weights are combined by the provided probability and
+-- event aggregation function.
+collectEventsBy :: (e -> e -> Ordering) -> ([(p,e)] -> (p,e))-> Categorical p e -> Categorical p e
+collectEventsBy compareE combine (Categorical ds) = 
+    Categorical . map combine . groupEvents . sortEvents $ ds
+    where
+        groupEvents = groupBy (\x y -> snd x `compareE` snd y == EQ)
+        sortEvents  = sortBy (compareE `on` snd)
diff --git a/src/Data/Random/Distribution/Dirichlet.hs b/src/Data/Random/Distribution/Dirichlet.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution/Dirichlet.hs
@@ -0,0 +1,30 @@
+{-# LANGUAGE
+    MultiParamTypeClasses,
+    FlexibleInstances, FlexibleContexts,
+    UndecidableInstances, GADTs
+  #-}
+
+module Data.Random.Distribution.Dirichlet where
+
+import Data.Random.RVar
+import Data.Random.Distribution
+import Data.Random.Distribution.Gamma
+
+import Data.List
+
+fractionalDirichlet :: (Fractional a, Distribution Gamma a) => [a] -> RVar [a]
+fractionalDirichlet []  = return []
+fractionalDirichlet [_] = return [1]
+fractionalDirichlet as = do
+    xs <- sequence [gamma a 1 | a <- as]
+    let total = foldl1' (+) xs
+    
+    return (map (* recip total) xs)
+
+dirichlet :: Distribution Dirichlet [a] => [a] -> RVar [a]
+dirichlet as = rvar (Dirichlet as)
+
+newtype Dirichlet a = Dirichlet a deriving (Eq, Show)
+
+instance (Fractional a, Distribution Gamma a) => Distribution Dirichlet [a] where
+    rvar (Dirichlet as) = fractionalDirichlet as
diff --git a/src/Data/Random/Distribution/Discrete.hs b/src/Data/Random/Distribution/Discrete.hs
deleted file mode 100644
--- a/src/Data/Random/Distribution/Discrete.hs
+++ /dev/null
@@ -1,136 +0,0 @@
-{-
- -      ``Data/Random/Distribution/Discrete''
- -}
-{-# LANGUAGE
-    MultiParamTypeClasses,
-    FlexibleInstances, FlexibleContexts
-  #-}
-
-module Data.Random.Distribution.Discrete where
-
-import Data.Random.RVar
-import Data.Random.Distribution
-import Data.Random.Distribution.Uniform
-import Data.Random.List (randomElement)
-
-import Control.Arrow
-import Control.Monad
-import Control.Applicative
-import Data.Foldable (Foldable(foldMap))
-import Data.Traversable (Traversable(traverse, sequenceA))
-
-import Data.List
-import Data.Function
-
-discrete :: Distribution (Discrete p) a => [(p,a)] -> RVar a
-discrete ps = rvar (Discrete ps)
-
-empirical :: (Num p, Ord a) => [a] -> Discrete p a
-empirical xs = Discrete bins
-    where bins = [ (genericLength bin, x)
-                 | bin@(x:_) <- group (sort xs)
-                 ]
-
-newtype Discrete p a = Discrete [(p, a)]
-    deriving (Eq, Show)
-
-instance (Num p, Ord p, Distribution Uniform p) => Distribution (Discrete p) a where
-    rvar (Discrete []) = fail "discrete distribution over empty set cannot be sampled"
-    rvar (Discrete ds) = do
-        let (ps, xs) = unzip ds
-            cs = scanl1 (+) ps
-        
-        when (any (<0) ps) $ fail "negative probability in discrete distribution"
-        
-        let totalWeight = last cs
-        if  totalWeight <= 0
-            -- if all events are "equally impossible", just pick an arbitrary one.
-            then randomElement xs
-            else do
-                let getU = do
-                        u <- uniform 0 totalWeight
-                        -- reject 0; this causes integral weights to behave sensibly (although it 
-                        -- potentially wastes up to 50% of all sampled integral values in the
-                        -- case where totalWeight = 1) and is still valid for fractional weights
-                        -- (it only prevents zero-probability events from ever occurring, which 
-                        -- is reasonable).
-                        if u == 0 then getU
-                                  else return u
-                u <- getU
-                return $ head
-                    [ x
-                    | (c,x) <- zip cs xs
-                    , c >= u
-                    ]
-
-instance Functor (Discrete p) where
-    fmap f (Discrete ds) = Discrete [(p, f x) | (p, x) <- ds]
-
-instance Foldable (Discrete p) where
-    foldMap f (Discrete ds) = foldMap (f . snd) ds
-
-instance Traversable (Discrete p) where
-    traverse f (Discrete ds) = Discrete <$> traverse (\(p,e) -> (\e -> (p,e)) <$> f e) ds
-    sequenceA  (Discrete ds) = Discrete <$> traverse (\(p,e) -> (\e -> (p,e)) <$>   e) ds
-
--- We want each subset of cases in fx derived from a given case 
--- in x to have the same relative weight as the set in x from whence they came.
-instance Num p => Monad (Discrete p) where
-    return x = Discrete [(1, x)]
-    (Discrete x) >>= f = Discrete $ do
-        (p, x) <- x
-        
-        let Discrete fx = f x
-        (q, x) <- fx
-        
-        return (p * q, x)
-
-instance Num p => Applicative (Discrete p) where
-    pure = return
-    (<*>) = ap
-
--- |Like 'fmap', but for the weights of a discrete distribution.
-mapDiscreteWeights :: (p -> q) -> Discrete p e -> Discrete q e
-mapDiscreteWeights f (Discrete ds) = Discrete [(f p, x) | (p, x) <- ds]
-
--- |Adjust all the weights of a discrete distribution so that they 
--- sum to unity.  If not possible, returns the original distribution 
--- unchanged.
-normalizeDiscreteWeights :: (Fractional p) => Discrete p e -> Discrete p e
-normalizeDiscreteWeights orig@(Discrete ds) = 
-    -- For practical purposes the scale factor is strict anyway,
-    -- so check if it's 0 or 1 and, if so, skip the actual scaling part.
-    if ws `elem` [0,1]
-        then orig
-        else Discrete
-                [ (w * scale, e)
-                | (w, e) <- ds
-                ] 
-    where
-        ws = sum (map fst ds)
-        scale = recip ws
-
--- |Simplify a discrete distribution by combining equivalent events (the new
--- event will have a weight equal to the sum of all the originals).
-collectDiscreteEvents :: (Ord e, Num p, Ord p) => Discrete p e -> Discrete p e
-collectDiscreteEvents = collectDiscreteEventsBy compare sum head
-        
--- |Simplify a discrete distribution by combining equivalent events (the new
--- event will have a weight equal to the sum of all the originals).
--- The comparator function is used to identify events to combine.  Once chosen,
--- the events and their weights are combined (independently) by the provided
--- weight and event aggregation functions.
-collectDiscreteEventsBy :: (e -> e -> Ordering) -> ([p] -> p) -> ([e] -> e)-> Discrete p e -> Discrete p e
-collectDiscreteEventsBy compareE sumWeights mergeEvents (Discrete ds) = 
-    Discrete . map ((sumWeights *** mergeEvents) . unzip) . groupEvents . sortEvents $ ds
-    
-    where
-        groupEvents = groupBy (\x y -> snd x `compareE` snd y == EQ)
-        sortEvents  = sortBy (compareE `on` snd)
-        
-        weight (p,x)
-            | p < 0     = error "negative probability in discrete distribution"
-            | otherwise = p
-        event ((p,x):_) = x
-        
-        combine (ps, xs) = (sumWeights ps, mergeEvents xs)
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
@@ -1,17 +1,7 @@
-{-
- -      ``Data/Random/Distribution/Gamma''
- -
- -  needs cleanup, verification, and automagic selection of appropriate
- -  algorithms, and proper citations.
- -
- -  should eliminate spurious 'border crossings' betweer RVars and sampleFrom
- -  perhaps Distribution class should have as its basis a function of type
- -  (d t -> RVar t)
- -}
 {-# LANGUAGE
     MultiParamTypeClasses,
     FlexibleInstances, FlexibleContexts,
-    UndecidableInstances
+    UndecidableInstances, BangPatterns
   #-}
 
 module Data.Random.Distribution.Gamma where
@@ -21,74 +11,44 @@
 import Data.Random.Distribution.Uniform
 import Data.Random.Distribution.Normal
 
-import Control.Monad
+import Data.Ratio
 
-    -- translated from gsl source - seems to be best I've found by far.
-    -- originally comes from Marsaglia & Tang, "A Simple Method for
-    -- generating gamma variables", ACM Transactions on Mathematical
-    -- Software, Vol 26, No 3 (2000), p363-372.
-realFloatGamma :: (Floating a, Ord a, Distribution Normal a, Distribution StdUniform a) => a -> a -> RVar a
-realFloatGamma a b
-    | a < 1 
-    = do
+-- derived from  Marsaglia & Tang, "A Simple Method for generating gamma
+-- variables", ACM Transactions on Mathematical Software, Vol 26, No 3 (2000), p363-372.
+{-# SPECIALIZE mtGamma :: Double -> Double -> RVar Double #-}
+{-# SPECIALIZE mtGamma :: Float -> Float -> RVar Float #-}
+mtGamma
+    :: (Floating a, Ord a,
+        Distribution StdUniform a, 
+        Distribution Normal a)
+    => a -> a -> RVar a
+mtGamma a b 
+    | a < 1     = do
         u <- stdUniform
-        x <- realFloatGamma (1 + a) b
-        return (x * u ** recip a)
-    | otherwise
-    = go
-        where
-            d = a - (1 / 3)
-            c = recip (3 * sqrt d) -- (1 / 3) / sqrt d
-            
-            go = do
-                x <- stdNormal
-                
-                let cx = c * x
-                    v = (1 + cx) ^ 3
-                    
-                    x_2 = x * x
-                    x_4 = x_2 * x_2
-                
-                if cx <= (-1)
-                    then go
-                    else do
-                        u <- stdUniform
-                        
-                        if         u < 1 - 0.0331 * x_4
-                            || log u < 0.5 * x_2  + d * (1 - v + log v)
-                            then return (b * d * v)
-                            else go
-
-
-realFloatErlang :: (Integral a, Floating b, Ord b, Distribution Normal b, Distribution StdUniform b) => a -> RVar b
-realFloatErlang a
-    | a < 1 
-    = fail "realFloatErlang: a < 1"
-    | otherwise
-    = go
-        where
-            d = fromIntegral a - (1 / 3)
-            c = recip (3 * sqrt d) -- (1 / 3) / sqrt d
+        mtGamma (1+a) $! (b * u ** recip a)
+    | otherwise = go
+    where
+        !d = a - fromRational (1%3)
+        !c = recip (sqrt (9*d))
+        
+        go = do
+            x <- stdNormal
+            let !v   = 1 + c*x
             
-            go = do
-                x <- stdNormal
-                
-                let cx = c * x
-                    v = (1 + cx) ^ 3
-                    
-                    x_2 = x * x
-                    x_4 = x_2 * x_2
-                
-                if cx <= (-1)
-                    then go
-                    else do
-                        u <- stdUniform
-                        
-                        if         u < 1 - 0.0331 * x_4
-                            || log u < 0.5 * x_2  + d * (1 - v + log v)
-                            then return (d * v)
-                            else go
+            if v <= 0
+                then go
+                else do
+                    u  <- stdUniform
+                    let !x_2 = x*x; !x_4 = x_2*x_2
+                        v3 = v*v*v
+                        dv = d * v3
+                    if      u < 1 - 0.0331*x_4
+                     || log u < 0.5 * x_2 + d - dv + d*log v3
+                        then return (b*dv)
+                        else go
 
+{-# SPECIALIZE gamma :: Float  -> Float  -> RVar Float  #-}
+{-# SPECIALIZE gamma :: Double -> Double -> RVar Double #-}
 gamma :: (Distribution Gamma a) => a -> a -> RVar a
 gamma a b = rvar (Gamma a b)
 
@@ -99,7 +59,9 @@
 data Erlang a b = Erlang a
 
 instance (Floating a, Ord a, Distribution Normal a, Distribution StdUniform a) => Distribution Gamma a where
-    rvar (Gamma a b) = realFloatGamma a b
+    {-# SPECIALIZE instance Distribution Gamma Double #-}
+    {-# SPECIALIZE instance Distribution Gamma Float #-}
+    rvar (Gamma a b) = mtGamma a b
 
 instance (Integral a, Floating b, Ord b, Distribution Normal b, Distribution StdUniform b) => Distribution (Erlang a) b where
-    rvar (Erlang a) = realFloatErlang a
+    rvar (Erlang a) = mtGamma (fromIntegral a) 1
diff --git a/src/Data/Random/Distribution/Multinomial.hs b/src/Data/Random/Distribution/Multinomial.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution/Multinomial.hs
@@ -0,0 +1,31 @@
+{-# LANGUAGE GADTs, MultiParamTypeClasses, FlexibleInstances, FlexibleContexts #-}
+module Data.Random.Distribution.Multinomial where
+
+import Data.Random.RVar
+import Data.Random.Distribution
+import Data.Random.Distribution.Binomial
+
+multinomial :: Distribution (Multinomial p) [a] => [p] -> a -> RVar [a]
+multinomial ps n = rvar (Multinomial ps n)
+
+data Multinomial p a where
+    Multinomial :: [p] -> a -> Multinomial p [a]
+
+instance (Num a, Fractional p, Distribution (Binomial p) a) => Distribution (Multinomial p) [a] where
+    -- TODO: implement faster version based on Categorical for small n, large (length ps)
+    rvar (Multinomial ps0 t) = go t ps0 (tailSums ps0) id
+        where
+            go _ []     _            f = return (f [])
+            go n [_]    _            f = return (f [n])
+            go 0 (_:ps) (_   :psums) f = go 0 ps psums (f . (0:))
+            go n (p:ps) (psum:psums) f = do
+                x <- binomial n (p / psum)
+                go (n-x) ps psums (f . (x:))
+            
+            go _ _ _ _ = error "rvar/Multinomial: programming error! this case should be impossible!"
+            
+            -- less wasteful version of (map sum . tails)
+            tailSums [] = [0]
+            tailSums (x:xs) = case tailSums xs of
+                (s:rest) -> (x+s):s:rest
+                _ -> error "rvar/Multinomial/tailSums: programming error! this case should be impossible!"
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
@@ -1,9 +1,6 @@
-{-
- -      ``Data/Random/Distribution/Normal''
- -}
 {-# LANGUAGE
     MultiParamTypeClasses, FlexibleInstances, FlexibleContexts,
-    UndecidableInstances, ForeignFunctionInterface
+    UndecidableInstances, ForeignFunctionInterface, BangPatterns
   #-}
 
 module Data.Random.Distribution.Normal
@@ -31,8 +28,9 @@
 import Data.Random.Distribution.Ziggurat
 import Data.Random.RVar
 
-import Control.Monad
-import Foreign.Storable
+import Data.Vector.Generic (Vector)
+import qualified Data.Vector as V
+import qualified Data.Vector.Unboxed as UV
 
 import Data.Number.Erf
 
@@ -82,12 +80,12 @@
 normalTail :: (Distribution StdUniform a, Floating a, Ord a) =>
               a -> RVar a
 normalTail r = go
-    where 
+    where
         go = do
-            u <- stdUniform
-            v <- stdUniform
-            let x = log u / r
-                y = log v
+            !u <- stdUniform
+            let !x = log u / r
+            !v <- stdUniform
+            let !y = log v
             if x*x + y+y > 0
                 then go
                 else return (r - x)
@@ -95,8 +93,8 @@
 -- |Construct a 'Ziggurat' for sampling a normal distribution, given
 -- @logBase 2 c@ and the 'zGetIU' implementation.
 normalZ ::
-  (RealFloat a, Erf a, Storable a, Distribution Uniform a, Integral b) =>
-  b -> RVar (Int, a) -> Ziggurat a
+  (RealFloat a, Erf a, Vector v a, Distribution Uniform a, Integral b) =>
+  b -> RVar (Int, a) -> Ziggurat v a
 normalZ p = mkZigguratRec True normalF normalFInv normalFInt normalFVol (2^p)
 
 -- | Ziggurat target function (upper half of a non-normalized gaussian PDF)
@@ -119,7 +117,7 @@
 -- |A random variable sampling from the standard normal distribution
 -- over any 'RealFloat' type (subject to the rest of the constraints -
 -- it builds and uses a 'Ziggurat' internally, which requires the 'Erf'
--- and 'Storable' classes).  
+-- class).  
 -- 
 -- Because it computes a 'Ziggurat', it is very expensive to use for
 -- just one evaluation, or even for multiple evaluations if not used and
@@ -132,18 +130,20 @@
 --
 -- As far as I know, this should be safe to use in any monomorphic
 -- @Distribution Normal@ instance declaration.
-realFloatStdNormal :: (RealFloat a, Erf a, Storable a, Distribution Uniform a) => RVar a
-realFloatStdNormal = runZiggurat (normalZ p getIU)
+realFloatStdNormal :: (RealFloat a, Erf a, Distribution Uniform a) => RVar a
+realFloatStdNormal = runZiggurat (normalZ p getIU `asTypeOf` (undefined :: Ziggurat V.Vector a))
     where 
+        p :: Int
         p = 6
         
         getIU = do
-            i <- getRandomByte
+            i <- getRandomPrim PrimWord8
             u <- uniform (-1) 1
             return (fromIntegral i .&. (2^p-1), u)
 
 -- |A random variable sampling from the standard normal distribution
 -- over the 'Double' type.
+{-# NOINLINE doubleStdNormal #-}
 doubleStdNormal :: RVar Double
 doubleStdNormal = runZiggurat doubleStdNormalZ
 
@@ -154,7 +154,8 @@
 doubleStdNormalR = 3.852046150368388
 doubleStdNormalV = 2.4567663515413507e-3
 
-doubleStdNormalZ :: Ziggurat Double
+{-# NOINLINE doubleStdNormalZ #-}
+doubleStdNormalZ :: Ziggurat UV.Vector Double
 doubleStdNormalZ = mkZiggurat_ True 
         normalF normalFInv 
         doubleStdNormalC doubleStdNormalR doubleStdNormalV 
@@ -162,23 +163,25 @@
         (normalTail doubleStdNormalR)
     where 
         getIU = do
-            w <- getRandomWord
+            !w <- getRandomPrim PrimWord64
             let (u,i) = wordToDoubleWithExcess w
-            return (fromIntegral i .&. (doubleStdNormalC-1), u+u-1)
+            return $! (fromIntegral i .&. (doubleStdNormalC-1), u+u-1)
 
 -- |A random variable sampling from the standard normal distribution
 -- over the 'Float' type.
+{-# NOINLINE floatStdNormal #-}
 floatStdNormal :: RVar Float
 floatStdNormal = runZiggurat floatStdNormalZ
 
--- floatStdNormalC must not be over 2^41 if using wordToFloatWithExcess
+-- floatStdNormalC must not be over 2^9 if using word32ToFloatWithExcess
 floatStdNormalC :: Int
 floatStdNormalC = 512
 floatStdNormalR, floatStdNormalV :: Float
 floatStdNormalR = 3.852046150368388
 floatStdNormalV = 2.4567663515413507e-3
 
-floatStdNormalZ :: Ziggurat Float
+{-# NOINLINE floatStdNormalZ #-}
+floatStdNormalZ :: Ziggurat UV.Vector Float
 floatStdNormalZ = mkZiggurat_ True 
         normalF normalFInv 
         floatStdNormalC floatStdNormalR floatStdNormalV 
@@ -186,13 +189,10 @@
         (normalTail floatStdNormalR)
     where
         getIU = do
-            w <- getRandomWord
-            let (u,i) = wordToFloatWithExcess w
+            !w <- getRandomPrim PrimWord32
+            let (u,i) = word32ToFloatWithExcess w
             return (fromIntegral i .&. (floatStdNormalC-1), u+u-1)
 
-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 = normcdf ((realToFrac x - realToFrac m) / realToFrac s)
 
@@ -200,7 +200,7 @@
 data Normal a
     -- |The \"standard\" normal distribution - mean 0, stddev 1
     = StdNormal
-    -- |@Normal m s@ is a normal distribution with mean @m@ and stddev @s@.
+    -- |@Normal m s@ is a normal distribution with mean @m@ and stddev @sd@.
     | Normal a a -- mean, sd
 
 instance Distribution Normal Double where
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
@@ -17,16 +17,14 @@
 import Data.Random.Distribution.Gamma
 import Data.Random.Distribution.Binomial
 
-import Data.Int
-import Data.Word
 import Control.Monad
 
 -- from Knuth, with interpretation help from gsl sources
 integralPoisson :: (Integral a, RealFloat b, Distribution StdUniform b, Distribution (Erlang a) b, Distribution (Binomial b) a) => b -> RVar a
-integralPoisson mu = psn 0 mu
+integralPoisson = psn 0
     where
         psn :: (Integral a, RealFloat b, Distribution StdUniform b, Distribution (Erlang a) b, Distribution (Binomial b) a) => a -> b -> RVar a
-        psn k mu
+        psn j mu
             | mu > 10   = do
                 let m = floor (mu * (7/8))
             
@@ -34,10 +32,10 @@
                 if x >= mu
                     then do
                         b <- binomial (m - 1) (mu / x)
-                        return (k + b)
-                    else psn (k + m) (mu - x)
+                        return (j + b)
+                    else psn (j + m) (mu - x)
             
-            | otherwise = prod 1 k
+            | otherwise = prod 1 j
                 where
                     emu = exp (-mu)
                 
diff --git a/src/Data/Random/Distribution/Triangular.hs b/src/Data/Random/Distribution/Triangular.hs
--- a/src/Data/Random/Distribution/Triangular.hs
+++ b/src/Data/Random/Distribution/Triangular.hs
@@ -28,12 +28,12 @@
     triUpper  :: a}
     deriving (Eq, Show)
 
--- |Compute a triangular distribution for a 'Fractional' type.  The name is
--- a historical accident and may change in the future.
-realFloatTriangular :: (Floating a, Ord a, Distribution StdUniform a) => a -> a -> a -> RVar a
-realFloatTriangular a b c
-    | a <= b && b <= c
-    = do
+-- |Compute a triangular distribution for a 'Floating' type.
+floatingTriangular :: (Floating a, Ord a, Distribution StdUniform a) => a -> a -> a -> RVar a
+floatingTriangular a b c
+    | a > b     = floatingTriangular b a c
+    | b > c     = floatingTriangular a c b
+    | otherwise = do
         let p = (c-b)/(c-a)
         u <- stdUniform
         let d   | u >= p    = a
@@ -44,19 +44,19 @@
 --        x <- stdUniform
         return (b - ((1 - sqrt x) * (b-d)))
 
--- |@realFloatTriangularCDF a b c@ is the CDF of @realFloatTriangular a b c@.
-realFloatTriangularCDF :: RealFrac a => a -> a -> a -> a -> Double
-realFloatTriangularCDF a b c x
+-- |@triangularCDF a b c@ is the CDF of @realFloatTriangular a b c@.
+triangularCDF :: RealFrac a => a -> a -> a -> a -> Double
+triangularCDF a b c x
     | x < a
     = 0
     | x <= b
-    = realToFrac ((x - a) ^ 2 / ((c - a) * (b - a)))
+    = realToFrac ((x - a)^(2 :: Int) / ((c - a) * (b - a)))
     | x <= c
-    = realToFrac (1 - (c - x) ^ 2 / ((c - a) * (c - b)))
+    = realToFrac (1 - (c - x)^(2 :: Int) / ((c - a) * (c - b)))
     | otherwise
     = 1
     
 instance (RealFloat a, Ord a, Distribution StdUniform a) => Distribution Triangular a where
-    rvar (Triangular a b c) = realFloatTriangular a b c
+    rvar (Triangular a b c) = floatingTriangular a b c
 instance (RealFrac a, Distribution Triangular a) => CDF Triangular a where
-    cdf  (Triangular a b c) = realFloatTriangularCDF a b c
+    cdf  (Triangular a b c) = triangularCDF a b c
diff --git a/src/Data/Random/Distribution/Uniform.hs b/src/Data/Random/Distribution/Uniform.hs
--- a/src/Data/Random/Distribution/Uniform.hs
+++ b/src/Data/Random/Distribution/Uniform.hs
@@ -5,7 +5,8 @@
     MultiParamTypeClasses, FunctionalDependencies,
     FlexibleContexts, FlexibleInstances, 
     UndecidableInstances, EmptyDataDecls,
-    TemplateHaskell
+    TemplateHaskell,
+    BangPatterns
   #-}
 
 module Data.Random.Distribution.Uniform
@@ -14,7 +15,6 @@
 	
     , StdUniform(..)
     , stdUniform
-    , stdUniformNonneg
     , stdUniformPos
     
     , integralUniform
@@ -50,30 +50,55 @@
 import Control.Monad.Loops
 
 -- |Compute a random 'Integral' value between the 2 values provided (inclusive).
+{-# INLINE integralUniform #-}
 integralUniform :: (Integral a) => a -> a -> RVar a
-integralUniform a b
-    | a > b     = compute b a
-    | otherwise = compute a b
+integralUniform !x !y = if x < y then integralUniform' x y else integralUniform' y x
+
+{-# SPECIALIZE integralUniform' :: Int   -> Int   -> RVar Int   #-}
+{-# SPECIALIZE integralUniform' :: Int8  -> Int8  -> RVar Int8  #-}
+{-# SPECIALIZE integralUniform' :: Int16 -> Int16 -> RVar Int16 #-}
+{-# SPECIALIZE integralUniform' :: Int32 -> Int32 -> RVar Int32 #-}
+{-# SPECIALIZE integralUniform' :: Int64 -> Int64 -> RVar Int64 #-}
+{-# SPECIALIZE integralUniform' :: Word   -> Word   -> RVar Word   #-}
+{-# SPECIALIZE integralUniform' :: Word8  -> Word8  -> RVar Word8  #-}
+{-# SPECIALIZE integralUniform' :: Word16 -> Word16 -> RVar Word16 #-}
+{-# SPECIALIZE integralUniform' :: Word32 -> Word32 -> RVar Word32 #-}
+{-# SPECIALIZE integralUniform' :: Word64 -> Word64 -> RVar Word64 #-}
+{-# SPECIALIZE integralUniform' :: Integer -> Integer -> RVar Integer #-}
+integralUniform' :: (Integral a) => a -> a -> RVar a
+integralUniform' !l !u
+    | nReject == 0  = fmap shift prim
+    | otherwise     = fmap shift loop
     where
-        compute a b = do
-            let m = 1 + toInteger b - toInteger a
-            
-            let bytes = bytesNeeded m
-                maxXpossible = (powersOf256 !! bytes) - 1
-            
-            x <- iterateUntil (maxXpossible - maxXpossible `mod` m >) (nByteInteger bytes)
-            return (a + fromInteger (x `mod` m))
+        m = 1 + toInteger u - toInteger l
+        (bytes, nPossible) = bytesNeeded m
+        nReject = nPossible `mod` m
+        
+        !prim = getRandomPrim (PrimNByteInteger bytes)
+        !shift = \(!z) -> l + (fromInteger $! (z `mod` m))
+        
+        loop = do
+            z <- prim
+            if z < nReject
+                then loop
+                else return z
 
+integralUniformCDF :: (Integral a, Fractional b) => a -> a -> a -> b
 integralUniformCDF a b x
     | b < a     = integralUniformCDF b a x
     | x < a     = 0
     | x > b     = 1
     | otherwise = (fromIntegral x - fromIntegral a) / (fromIntegral b - fromIntegral a)
 
-bytesNeeded x = case findIndex (> x) powersOf256 of
-    Just x -> x
-powersOf256 = iterate (256 *) 1
+-- TODO: come up with a decent, fast heuristic to decide whether to return an extra
+-- byte.  May involve moving calculation of nReject into this function, and then
+-- accepting first if 4*nReject < nPossible or something similar.
+bytesNeeded :: Integer -> (Int, Integer)
+bytesNeeded x = head (dropWhile ((<= x).snd) powersOf256)
 
+powersOf256 :: [(Int, Integer)]
+powersOf256 = zip [0..] (iterate (256 *) 1)
+
 -- |Compute a random value for a 'Bounded' type, between 'minBound' and 'maxBound'
 -- (inclusive for 'Integral' or 'Enum' types, in ['minBound', 'maxBound') for Fractional types.)
 boundedStdUniform :: (Distribution Uniform a, Bounded a) => RVar a
@@ -93,12 +118,13 @@
 -- |Compute a uniform random 'Float' value in the range [0,1)
 floatStdUniform :: RVar Float
 floatStdUniform = do
-    x <- getRandomWord
-    return (wordToFloat x)
+    x <- getRandomPrim PrimWord32
+    return (word32ToFloat x)
 
 -- |Compute a uniform random 'Double' value in the range [0,1)
+{-# INLINE doubleStdUniform #-}
 doubleStdUniform :: RVar Double
-doubleStdUniform = getRandomDouble
+doubleStdUniform = getRandomPrim PrimDouble
 
 -- |Compute a uniform random value in the range [0,1) for any 'RealFloat' type 
 realFloatStdUniform :: RealFloat a => RVar a
@@ -137,6 +163,7 @@
     return (a + x * (b - a))
 
 -- |@doubleUniform a b@ computes a uniform random 'Double' value in the range [a,b)
+{-# INLINE doubleUniform #-}
 doubleUniform :: Double -> Double -> RVar Double
 doubleUniform 0 1 = doubleStdUniform
 doubleUniform a b = do
@@ -190,22 +217,22 @@
 
 -- |Get a \"standard\" uniformly distributed value.
 -- For integral types, this means uniformly distributed over the full range
--- of the type (and hence there is no support for Integer).  For fractional
+-- of the type (there is no support for 'Integer').  For fractional
 -- types, this means uniformly distributed on the interval [0,1).
+{-# SPECIALIZE stdUniform :: RVar Double #-}
+{-# SPECIALIZE stdUniform :: RVar Float #-}
 stdUniform :: (Distribution StdUniform a) => RVar a
 stdUniform = rvar StdUniform
 
--- |Like 'stdUniform', but uses 'abs' to return only positive or zero values.
+-- |Like 'stdUniform', but returns only positive or zero values.  Not 
+-- exported because it is not truly uniform: nonzero values are twice
+-- as likely as zero on signed types.
 stdUniformNonneg :: (Distribution StdUniform a, Num a) => RVar a
-stdUniformNonneg = abs `fmap` stdUniform
+stdUniformNonneg = fmap abs stdUniform
 
 -- |Like 'stdUniform' but only returns positive values.
 stdUniformPos :: (Distribution StdUniform a, Num a) => RVar a
-stdUniformPos = do
-    x <- stdUniformNonneg
-    if x /= 0
-        then return x
-        else stdUniformPos
+stdUniformPos = iterateUntil (/= 0) stdUniformNonneg
 
 -- |A definition of a uniform distribution over the type @t@.  See also 'uniform'.
 data Uniform t = 
@@ -229,15 +256,34 @@
         instance CDF Uniform Int            where cdf  (Uniform a b) = integralUniformCDF a b
     |])
 
--- Some integral types have specialized StdUniform rvars:
-instance Distribution StdUniform Int8       where rvar ~StdUniform = fmap fromIntegral getRandomByte
-instance Distribution StdUniform Word8      where rvar ~StdUniform = getRandomByte
-instance Distribution StdUniform Word64     where rvar ~StdUniform = getRandomWord
--- and Integer has none...
-$( replicateInstances ''Int (integralTypes \\ [''Int8, ''Word8, ''Word64, ''Integer]) [d|
-        instance Distribution StdUniform Int    where rvar ~StdUniform = fmap fromIntegral getRandomWord
-    |])
+instance Distribution StdUniform Word8      where rvarT ~StdUniform = getRandomPrim PrimWord8
+instance Distribution StdUniform Word16     where rvarT ~StdUniform = getRandomPrim PrimWord16
+instance Distribution StdUniform Word32     where rvarT ~StdUniform = getRandomPrim PrimWord32
+instance Distribution StdUniform Word64     where rvarT ~StdUniform = getRandomPrim PrimWord64
 
+instance Distribution StdUniform Int8       where rvarT ~StdUniform = fromIntegral `fmap` getRandomPrim PrimWord8
+instance Distribution StdUniform Int16      where rvarT ~StdUniform = fromIntegral `fmap` getRandomPrim PrimWord16
+instance Distribution StdUniform Int32      where rvarT ~StdUniform = fromIntegral `fmap` getRandomPrim PrimWord32
+instance Distribution StdUniform Int64      where rvarT ~StdUniform = fromIntegral `fmap` getRandomPrim PrimWord64
+
+instance Distribution StdUniform Int where
+    rvar
+        | toInteger (maxBound :: Int) > toInteger (maxBound :: Int32)
+        = const (fromIntegral `fmap` getRandomPrim PrimWord64)
+        
+        | otherwise
+        = const (fromIntegral `fmap` getRandomPrim PrimWord32)
+
+instance Distribution StdUniform Word where
+    rvar
+        | toInteger (maxBound :: Word) > toInteger (maxBound :: Word32)
+        = const (fromIntegral `fmap` getRandomPrim PrimWord64)
+        
+        | otherwise
+        = const (fromIntegral `fmap` getRandomPrim PrimWord32)
+
+-- Integer has no StdUniform...
+
 $( replicateInstances ''Int (integralTypes \\ [''Integer]) [d|
         instance CDF StdUniform Int         where cdf  ~StdUniform = boundedStdUniformCDF
     |])
@@ -249,7 +295,7 @@
 instance CDF Uniform Double                 where cdf  (Uniform a b) = realUniformCDF a b
 
 instance Distribution StdUniform Float      where rvar ~StdUniform = floatStdUniform
-instance Distribution StdUniform Double     where rvar ~StdUniform = doubleStdUniform
+instance Distribution StdUniform Double     where rvar ~StdUniform = getRandomPrim PrimDouble; rvarT ~StdUniform = getRandomPrim PrimDouble
 instance CDF StdUniform Float               where cdf  ~StdUniform = realStdUniformCDF
 instance CDF StdUniform Double              where cdf  ~StdUniform = realStdUniformCDF
 
@@ -262,18 +308,18 @@
 instance HasResolution r => 
          CDF StdUniform (Fixed r)           where cdf  ~StdUniform = realStdUniformCDF
 
-instance Distribution Uniform ()            where rvar (Uniform a b) = return ()
-instance CDF Uniform ()                     where cdf  (Uniform a b) = return 1
+instance Distribution Uniform ()            where rvar (Uniform _ _) = return ()
+instance CDF Uniform ()                     where cdf  (Uniform _ _) = return 1
 $( replicateInstances ''Char [''Char, ''Bool, ''Ordering] [d|
         instance Distribution Uniform Char  where rvar (Uniform a b) = enumUniform a b
         instance CDF Uniform Char           where cdf  (Uniform a b) = enumUniformCDF a b
 
     |])
 
-instance Distribution StdUniform ()         where rvar ~StdUniform = return ()
-instance CDF StdUniform ()                  where cdf  ~StdUniform = return 1
-instance Distribution StdUniform Bool       where rvar ~StdUniform = fmap even getRandomByte
-instance CDF StdUniform Bool                where cdf  ~StdUniform = boundedEnumStdUniformCDF
+instance Distribution StdUniform ()         where rvarT ~StdUniform = return ()
+instance CDF StdUniform ()                  where cdf   ~StdUniform = return 1
+instance Distribution StdUniform Bool       where rvarT ~StdUniform = fmap even (getRandomPrim PrimWord8)
+instance CDF StdUniform Bool                where cdf   ~StdUniform = boundedEnumStdUniformCDF
 
 instance Distribution StdUniform Char       where rvar ~StdUniform = boundedEnumStdUniform
 instance CDF StdUniform Char                where cdf  ~StdUniform = boundedEnumStdUniformCDF
diff --git a/src/Data/Random/Distribution/Weibull.hs b/src/Data/Random/Distribution/Weibull.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Distribution/Weibull.hs
@@ -0,0 +1,16 @@
+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, UndecidableInstances #-}
+module Data.Random.Distribution.Weibull where
+
+import Data.Random.Distribution
+import Data.Random.Distribution.Uniform
+
+data Weibull a = Weibull { weibullLambda :: !a, weibullK :: !a }
+    deriving (Eq, Show)
+
+instance (Floating a, Distribution StdUniform a) => Distribution Weibull a where
+    rvarT (Weibull lambda k) = do
+        u <- rvarT StdUniform
+        return (lambda * (negate (log u)) ** recip k)
+
+instance (Real a, Distribution Weibull a) => CDF Weibull a where
+    cdf (Weibull lambda k) x = 1 - exp (negate ((realToFrac x / realToFrac lambda) ** realToFrac k))
diff --git a/src/Data/Random/Distribution/Ziggurat.hs b/src/Data/Random/Distribution/Ziggurat.hs
--- a/src/Data/Random/Distribution/Ziggurat.hs
+++ b/src/Data/Random/Distribution/Ziggurat.hs
@@ -1,7 +1,7 @@
 {-# LANGUAGE
         MultiParamTypeClasses,
         FlexibleInstances, FlexibleContexts,
-        RecordWildCards
+        RecordWildCards, BangPatterns
   #-}
 
 -- |A generic \"ziggurat algorithm\" implementation.  Fairly rough right
@@ -31,10 +31,10 @@
 import Data.Random.Distribution.Uniform
 import Data.Random.Distribution
 import Data.Random.RVar
-import Data.StorableVector as Vec
-import Foreign.Storable
-
-vec ! i = index vec i
+import Data.Vector.Generic as Vec
+import qualified Data.Vector as V
+import qualified Data.Vector.Unboxed as UV
+import Data.Function (fix)
 
 -- |A data structure containing all the data that is needed
 -- to implement Marsaglia & Tang's \"ziggurat\" algorithm for
@@ -45,7 +45,7 @@
 -- be.  There are several helper functions that will build 'Ziggurat's.
 -- The pathologically curious may wish to read the 'runZiggurat' source.
 -- That is the ultimate specification of the semantics of all these fields.
-data Ziggurat t = Ziggurat {
+data Ziggurat v t = Ziggurat {
         -- |The X locations of each bin in the distribution.  Bin 0 is the
         -- 'infinite' one.
         -- 
@@ -55,11 +55,11 @@
         -- faster by not needing to specially-handle bin 0 quite as often.
         -- If you really need to know why it works, see the 'runZiggurat'
         -- source or \"the literature\" - it's a fairly standard setup.
-        zTable_xs         :: Vector t,
+        zTable_xs         :: !(v t),
         -- |The ratio of each bin's Y value to the next bin's Y value
-        zTable_x_ratios   :: Vector t,
+        zTable_y_ratios   :: !(v t),
         -- |The Y value (zFunc x) of each bin
-        zTable_ys         :: Vector t,
+        zTable_ys         :: !(v t),
         -- |An RVar providing a random tuple consisting of:
         --
         --  * a bin index, uniform over [0,c) :: Int (where @c@ is the
@@ -73,58 +73,63 @@
         -- single random word (64 bits) can be efficiently converted to
         -- a double (using 52 bits) and a bin number (using up to 12 bits),
         -- for example.
-        zGetIU            :: RVar (Int, t),
+        zGetIU            :: !(RVar (Int, t)),
         
         -- |The distribution for the final \"virtual\" bin
         -- (the ziggurat algorithm does not handle distributions
         -- that wander off to infinity, so another distribution is needed
         -- to handle the last \"bin\" that stretches to infinity)
-        zTailDist         :: RVar t,
+        zTailDist         :: (RVar t),
         
         -- |A copy of the uniform RVar generator for the base type,
         -- so that @Distribution Uniform t@ is not needed when sampling
         -- from a Ziggurat (makes it a bit more self-contained).
-        zUniform          :: t -> t -> RVar t,
+        zUniform          :: !(t -> t -> RVar t),
         
         -- |The (one-sided antitone) PDF, not necessarily normalized
-        zFunc             :: t -> t,
+        zFunc             :: !(t -> t),
         
         -- |A flag indicating whether the distribution should be
         -- mirrored about the origin (the ziggurat algorithm it
         -- its native form only samples from one-sided distributions.
         -- By mirroring, we can extend it to symmetric distributions
         -- such as the normal distribution)
-        zMirror           :: Bool
+        zMirror           :: !Bool
     }
 
 -- |Sample from the distribution encoded in a 'Ziggurat' data structure.
 {-# INLINE runZiggurat #-}
-runZiggurat :: (Num a, Ord a, Storable a) =>
-               Ziggurat a -> RVar a
-runZiggurat Ziggurat{..} = go
+{-# SPECIALIZE runZiggurat :: Ziggurat UV.Vector Float  -> RVar Float #-}
+{-# SPECIALIZE runZiggurat :: Ziggurat UV.Vector Double -> RVar Double #-}
+{-# SPECIALIZE runZiggurat :: Ziggurat  V.Vector Float  -> RVar Float #-}
+{-# SPECIALIZE runZiggurat :: Ziggurat  V.Vector Double -> RVar Double #-}
+runZiggurat :: (Num a, Ord a, Vector v a) =>
+               Ziggurat v a -> RVar a
+runZiggurat !Ziggurat{..} = go
     where
+        {-# NOINLINE go #-}
         go = do
             -- Select a bin (I) and a uniform value (U) from -1 to 1
             -- (or 0 to 1 if not mirroring the distribution).
             -- Let X be U scaled to the size of the selected bin.
-            (i,u) <- zGetIU
-            let x = u * zTable_xs ! i
+            (!i,!u) <- zGetIU
             
             -- if the uniform value U falls in the area "clearly inside" the
             -- bin, accept X immediately.
             -- Otherwise, depending on the bin selected, use either the
             -- tail distribution or an accept/reject test.
-            if abs u < zTable_x_ratios ! i
-                then return $! x
+            if abs u < zTable_y_ratios ! i
+                then return $! (u * zTable_xs ! i)
                 else if i == 0
-                    then sampleTail x
-                    else sampleGreyArea i x
+                    then sampleTail u
+                    else sampleGreyArea i $! (u * zTable_xs ! i)
         
         -- when the sample falls in the "grey area" (the area between
         -- the Y values of the selected bin and the bin after that one),
         -- use an accept/reject method based on the target PDF.
+        {-# INLINE sampleGreyArea #-}
         sampleGreyArea i x = do
-            v <- zUniform (zTable_ys ! (i+1)) (zTable_ys ! i)
+            !v <- zUniform (zTable_ys ! (i+1)) (zTable_ys ! i)
             if v < zFunc (abs x)
                 then return $! x
                 else go
@@ -132,9 +137,10 @@
         -- if the selected bin is the "infinite" one, call it quits and
         -- defer to the tail distribution (mirroring if needed to ensure
         -- the result has the sign already selected by zGetIU)
+        {-# INLINE sampleTail #-}
         sampleTail x
-            | x < 0     = fmap negate zTailDist
-            | otherwise = zTailDist
+            | zMirror && x < 0  = fmap negate zTailDist
+            | otherwise         = zTailDist
 
 
 -- |Build the tables to implement the \"ziggurat algorithm\" devised by 
@@ -159,7 +165,12 @@
 -- 
 --  * an RVar sampling from the tail (the region where x > R)
 -- 
-mkZiggurat_ :: (RealFloat t, Storable t,
+{-# INLINE mkZiggurat_ #-}
+{-# SPECIALIZE mkZiggurat_ :: Bool -> (Float  ->  Float) -> (Float  ->  Float) -> Int -> Float  -> Float  -> RVar (Int,  Float) -> RVar Float  -> Ziggurat UV.Vector Float #-}
+{-# SPECIALIZE mkZiggurat_ :: Bool -> (Double -> Double) -> (Double -> Double) -> Int -> Double -> Double -> RVar (Int, Double) -> RVar Double -> Ziggurat UV.Vector Double #-}
+{-# SPECIALIZE mkZiggurat_ :: Bool -> (Float  ->  Float) -> (Float  ->  Float) -> Int -> Float  -> Float  -> RVar (Int,  Float) -> RVar Float  -> Ziggurat V.Vector Float #-}
+{-# SPECIALIZE mkZiggurat_ :: Bool -> (Double -> Double) -> (Double -> Double) -> Int -> Double -> Double -> RVar (Int, Double) -> RVar Double -> Ziggurat V.Vector Double #-}
+mkZiggurat_ :: (RealFloat t, Vector v t,
                Distribution Uniform t) =>
               Bool
               -> (t -> t)
@@ -169,18 +180,19 @@
               -> t
               -> RVar (Int, t)
               -> RVar t
-              -> Ziggurat t
-mkZiggurat_ m f fInv c r v getIU tailDist = z
-    where z = Ziggurat
-            { zTable_xs         = zigguratTable f fInv c r v
-            , zTable_x_ratios   = precomputeRatios (zTable_xs z)
-            , zTable_ys         = Vec.map f (zTable_xs z)
-            , zGetIU            = getIU
-            , zUniform          = uniform
-            , zFunc             = f
-            , zTailDist         = tailDist
-            , zMirror           = m
-            }
+              -> Ziggurat v t
+mkZiggurat_ m f fInv c r v getIU tailDist = Ziggurat
+    { zTable_xs         = xs
+    , zTable_y_ratios   = precomputeRatios xs
+    , zTable_ys         = Vec.map f xs
+    , zGetIU            = getIU
+    , zUniform          = uniform
+    , zFunc             = f
+    , zTailDist         = tailDist
+    , zMirror           = m
+    }
+    where 
+        xs = zigguratTable f fInv c r v
 
 -- |Build the tables to implement the \"ziggurat algorithm\" devised by 
 -- Marsaglia & Tang, attempting to automatically compute the R and V
@@ -189,7 +201,7 @@
 -- Arguments are the same as for 'mkZigguratRec', with an additional
 -- argument for the tail distribution as a function of the selected
 -- R value.
-mkZiggurat :: (RealFloat t, Storable t,
+mkZiggurat :: (RealFloat t, Vector v t,
                Distribution Uniform t) =>
               Bool
               -> (t -> t)
@@ -199,7 +211,7 @@
               -> Int
               -> RVar (Int, t)
               -> (t -> RVar t)
-              -> Ziggurat t
+              -> Ziggurat v t
 mkZiggurat m f fInv fInt fVol c getIU tailDist =
     mkZiggurat_ m f fInv c r v getIU (tailDist r) 
         where
@@ -226,7 +238,7 @@
 --  * an RVar providing the 'zGetIU' random tuple
 --
 mkZigguratRec ::
-  (RealFloat t, Storable t,
+  (RealFloat t, Vector v t,
    Distribution Uniform t) =>
   Bool
   -> (t -> t)
@@ -235,14 +247,23 @@
   -> t
   -> Int
   -> RVar (Int, t)
-  -> Ziggurat t
-mkZigguratRec m f fInv fInt fVol c getIU = 
-    mkZiggurat m f fInv fInt fVol c getIU (fix (mkTail m f fInv fInt fVol c getIU))
+  -> Ziggurat v t
+mkZigguratRec m f fInv fInt fVol c getIU = z
         where
-            fix f = f (fix f)
+            z = mkZiggurat m f fInv fInt fVol c getIU (fix (mkTail m f fInv fInt fVol c getIU z))
 
-mkTail m f fInv fInt fVol c getIU nextTail r = do
-     x <- rvar (mkZiggurat m f' fInv' fInt' fVol' c getIU nextTail)
+mkTail :: 
+    (RealFloat a, Vector v a, Distribution Uniform a) =>
+    Bool
+    -> (a -> a) -> (a -> a) -> (a -> a)
+    -> a
+    -> Int
+    -> RVar (Int, a)
+    -> Ziggurat v a
+    -> (a -> RVar a)
+    -> (a -> RVar a)
+mkTail m f fInv fInt fVol c getIU typeRep nextTail r = do
+     x <- rvar (mkZiggurat m f' fInv' fInt' fVol' c getIU nextTail `asTypeOf` typeRep)
      return (x + r * signum x)
         where
             fIntR = fInt r
@@ -256,14 +277,15 @@
             fVol' = fVol - fIntR
         
 
-zigguratTable :: (Fractional a, Storable a, Ord a) =>
-                 (a -> a) -> (a -> a) -> Int -> a -> a -> Vector a
+zigguratTable :: (Fractional a, Vector v a, Ord a) =>
+                 (a -> a) -> (a -> a) -> Int -> a -> a -> v a
 zigguratTable f fInv c r v = case zigguratXs f fInv c r v of
-    (xs, excess) -> pack xs
-    where epsilon = 1e-3*v
+    (xs, _excess) -> fromList xs
 
+zigguratExcess :: (Fractional a, Ord a) => (a -> a) -> (a -> a) -> Int -> a -> a -> a
 zigguratExcess f fInv c r v = snd (zigguratXs f fInv c r v)
 
+zigguratXs :: (Fractional a, Ord a) => (a -> a) -> (a -> a) -> Int -> a -> a -> ([a], a)
 zigguratXs f fInv c r v = (xs, excess)
     where
         xs = Prelude.map x [0..c] -- sample c x
@@ -273,6 +295,7 @@
         x 1 = r
         x i | i == c = 0
         x (i+1) = next i
+        x _ = error "zigguratXs: programming error! this case should be impossible!"
         
         next i = let x_i = xs!!i
                   in if x_i <= 0 then -1 else fInv (ys!!i + (v / x_i))
@@ -280,7 +303,8 @@
         excess = xs!!(c-1) * (f 0 - ys !! (c-1)) - v 
 
 
-precomputeRatios zTable_xs = sample (c-1) $ \i -> zTable_xs!(i+1) / zTable_xs!i
+precomputeRatios :: (Vector v a, Fractional a) => v a -> v a
+precomputeRatios zTable_xs = generate (c-1) $ \i -> zTable_xs!(i+1) / zTable_xs!i
     where
         c = Vec.length zTable_xs
 
@@ -302,7 +326,7 @@
 -- Result: (R,V)
 findBin0 :: (RealFloat b) => 
     Int -> (b -> b) -> (b -> b) -> (b -> b) -> b -> (b, b)
-findBin0 cInt f fInv fInt fVol = (r,v r)
+findBin0 cInt f fInv fInt fVol = (rMin,v rMin)
     where
         c = fromIntegral cInt
         v r = r * f r + fVol - fInt r
@@ -310,11 +334,11 @@
         -- initial R guess:
         r0 = findMin (\r -> v r <= fVol / c)
         -- find a better R:
-        r = findMinFrom r0 1 $ \r -> 
+        rMin = findMinFrom r0 1 $ \r -> 
             let e = exc r 
              in e >= 0 && not (isNaN e)
         
         exc x = zigguratExcess f fInv cInt x (v x)
 
-instance (Num t, Ord t, Storable t) => Distribution Ziggurat t where
+instance (Num t, Ord t, Vector v t) => Distribution (Ziggurat v) t where
     rvar = runZiggurat
diff --git a/src/Data/Random/Internal/Find.hs b/src/Data/Random/Internal/Find.hs
--- a/src/Data/Random/Internal/Find.hs
+++ b/src/Data/Random/Internal/Find.hs
@@ -9,21 +9,25 @@
 findMax :: (Fractional a, Ord a) => (a -> Bool) -> a
 findMax p = negate (findMin (p.negate))
 
-findMaxFrom :: (Fractional a, Ord a) => a -> a -> (a -> Bool) -> a
-findMaxFrom z 0 p = findMaxFrom z 1 p
-findMaxFrom z step1 p = findMinFrom z (negate step1) p
-
 -- |Given an upward-closed predicate on an ordered Fractional type,
 -- find the smallest value satisfying the predicate.
 findMin :: (Fractional a, Ord a) => (a -> Bool) -> a
 findMin = findMinFrom 0 1
 
+-- |Given an upward-closed predicate on an ordered Fractional type,
+-- find the smallest value satisfying the predicate.  Starts at the
+-- specified point with the specified stepsize, performs an exponential
+-- search out from there until it finds an interval bracketing the
+-- change-point of the predicate, and then performs a bisection search
+-- to isolate the change point.  Note that infinitely-divisible domains 
+-- such as 'Rational' cannot be searched by this function because it does
+-- not terminate until it reaches a point where further subdivision of the
+-- interval has no effect.
 findMinFrom :: (Fractional a, Ord a) => a -> a -> (a -> Bool) -> a
-findMinFrom z 0 p = findMinFrom z 1 p
-findMinFrom z step1 p
-    | p z   = descend (z-step1) z
-    | otherwise
-    = fixZero (ascend z (z+step1))
+findMinFrom z0 0 p = findMinFrom z0 1 p
+findMinFrom z0 step1 p
+    | p z0      = descend (z0-step1) z0
+    | otherwise = fixZero (ascend z0 (z0+step1))
     where
         -- eliminate negative zero, which, in many domains, is technically
         -- a feasible answer
@@ -35,13 +39,13 @@
         -- 0 <= l < x
         ascend l x 
             | p x       = bisect l x
-            | otherwise = ascend x $! 2*x-z
+            | otherwise = ascend x $! 2*x-z0
         
         -- preconditions:
         -- p h
         -- x < h <= 0
         descend x h 
-            | p x       = (descend $! 2*x-z) x
+            | p x       = (descend $! 2*x-z0) x
             | otherwise = bisect x h
         
         -- preconditions:
@@ -49,11 +53,11 @@
         -- p h
         -- l <= h
         bisect l h 
-            | l >= h    = h
-            | l >= mid || mid >= h
+            | l /< h    = h
+            | l /< mid || mid /< h
             = if p mid then mid else h
             | p mid     = bisect l mid
             | otherwise = bisect mid h
             where 
-                a >= b = not (a < b)
+                a /< b = not (a < b)
                 mid = (l+h)*0.5
diff --git a/src/Data/Random/Internal/Primitives.hs b/src/Data/Random/Internal/Primitives.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Internal/Primitives.hs
@@ -0,0 +1,235 @@
+{-# LANGUAGE GADTs, RankNTypes, DeriveDataTypeable #-}
+-- |This is an experimental interface to support an extensible set of primitives,
+-- where a RandomSource will be able to support whatever subset of them they want
+-- and have well-founded defaults generated automatically for any unsupported
+-- primitives.
+--
+-- The purpose, in case it's not clear, is to decouple the implementations of
+-- entropy sources from any particular set of primitives, so that implementors
+-- of random variates can make use of a large number of primitives, supported
+-- on all entropy sources, while the burden on entropy-source implementors
+-- is only to provide one or two basic primitives of their choice.
+-- 
+-- One challenge I foresee with this interface is optimization - different 
+-- compilers or even different versions of GHC may treat this interface 
+-- radically differently, making it very hard to achieve reliable performance
+-- on all platforms.  It may even be that no compiler optimizes sufficiently
+-- to make the flexibility this system provides worth the overhead.  I hope
+-- this is not the case, but if it turns out to be a major problem, this
+-- system may disappear or be modified in significant ways.
+module Data.Random.Internal.Primitives (Prim(..), decomposePrimWhere) where
+
+import Data.Random.Internal.Words
+import Data.Word
+import Data.Bits
+import Data.Typeable
+
+import Control.Monad.Prompt
+
+-- |A 'Prompt' GADT describing a request for a primitive random variate.
+-- Random variable definitions will request their entropy via these prompts,
+-- and entropy sources will satisfy some or all of them.  The 'decomposePrimWhere'
+-- function extends an entropy source's incomplete definition to a complete 
+-- definition, essentially defining a very flexible implementation-defaulting
+-- system.
+-- 
+-- Some possible future additions:
+--    PrimFloat :: Prim Float
+--    PrimInt :: Prim Int
+--    PrimPair :: Prim a -> Prim b -> Prim (a :*: b)
+--    PrimNormal :: Prim Double
+--    PrimChoice :: [(Double :*: a)] -> Prim a
+--
+-- Unfortunately, I cannot get Haddock to accept my comments about the 
+-- data constructors, but hopefully they should be reasonably self-explanatory.
+data Prim a where
+    -- An unsigned byte, uniformly distributed from 0 to 0xff
+    PrimWord8           :: Prim Word8
+    -- An unsigned 16-bit word, uniformly distributed from 0 to 0xffff
+    PrimWord16          :: Prim Word16
+    -- An unsigned 32-bit word, uniformly distributed from 0 to 0xffffffff
+    PrimWord32          :: Prim Word32
+    -- An unsigned 64-bit word, uniformly distributed from 0 to 0xffffffffffffffff
+    PrimWord64          :: Prim Word64
+    -- A double-precision float U, uniformly distributed 0 <= U < 1
+    PrimDouble          :: Prim Double
+    -- A uniformly distributed 'Integer' 0 <= U < 2^(8*n)
+    PrimNByteInteger    :: !Int -> Prim Integer
+    deriving (Typeable)
+
+instance Show (Prim a) where
+    showsPrec _p PrimWord8               = showString "PrimWord8"
+    showsPrec _p PrimWord16              = showString "PrimWord16"
+    showsPrec _p PrimWord32              = showString "PrimWord32"
+    showsPrec _p PrimWord64              = showString "PrimWord64"
+    showsPrec _p PrimDouble              = showString "PrimDouble"
+    showsPrec  p (PrimNByteInteger n)    = showParen (p > 10) (showString "PrimNByteInteger " . showsPrec 11 n)
+
+-- |This is essentially a suite of interrelated default implementations,
+-- each definition making use of only \"supported\" primitives.  It _really_
+-- ought to be inlined to the point where the @supported@ predicate
+-- is able to be inlined into it and eliminated.  
+-- 
+-- When inlined sufficiently, it should in theory be optimized down to the
+-- static set of "best" definitions for each required primitive in terms of 
+-- only supported primitives.
+-- 
+-- Hopefully, when not inlined, it does not impose too much overhead.
+{-# INLINE decomposePrimWhere #-}
+{-# SPECIALIZE decomposePrimWhere :: (forall t. Prim t -> Bool) -> Prim Word8   -> Prompt Prim Word8   #-}
+{-# SPECIALIZE decomposePrimWhere :: (forall t. Prim t -> Bool) -> Prim Word16  -> Prompt Prim Word16  #-}
+{-# SPECIALIZE decomposePrimWhere :: (forall t. Prim t -> Bool) -> Prim Word32  -> Prompt Prim Word32  #-}
+{-# SPECIALIZE decomposePrimWhere :: (forall t. Prim t -> Bool) -> Prim Word64  -> Prompt Prim Word64  #-}
+{-# SPECIALIZE decomposePrimWhere :: (forall t. Prim t -> Bool) -> Prim Double  -> Prompt Prim Double  #-}
+{-# SPECIALIZE decomposePrimWhere :: (forall t. Prim t -> Bool) -> Prim Integer -> Prompt Prim Integer #-}
+decomposePrimWhere :: (forall t. Prim t -> Bool) -> Prim a -> Prompt Prim a
+decomposePrimWhere supported requested = decomp requested
+    where
+        {-# INLINE decomp #-}
+
+        {-# SPECIALIZE decomp :: Prim Word8   -> Prompt Prim Word8   #-}
+        {-# SPECIALIZE decomp :: Prim Word16  -> Prompt Prim Word16  #-}
+        {-# SPECIALIZE decomp :: Prim Word32  -> Prompt Prim Word32  #-}
+        {-# SPECIALIZE decomp :: Prim Word64  -> Prompt Prim Word64  #-}
+        {-# SPECIALIZE decomp :: Prim Double  -> Prompt Prim Double  #-}
+        {-# SPECIALIZE decomp :: Prim Integer -> Prompt Prim Integer #-}
+        -- First, all supported prims should just be evaluated directly.
+        decomp :: Prim a -> Prompt Prim a
+        decomp prim
+            | supported prim = prompt prim
+        -- beyond this point, all definitions must be in terms of
+        -- 'prompt's referring to other supported primitives or 
+        -- 'decomp's referring to other primitives in a well-founded way
+        
+        decomp PrimWord8
+            | supported PrimWord16 = do
+                w <- prompt PrimWord16
+                return (fromIntegral w)
+            | supported PrimWord32 = do
+                w <- prompt PrimWord32
+                return (fromIntegral w)
+            | supported PrimWord64 = do
+                w <- prompt PrimWord64
+                return (fromIntegral w)
+            | supported PrimDouble = do
+                d <- prompt PrimDouble
+                return (truncate (d * 256))
+            | supported (PrimNByteInteger 1) = do
+                i <- prompt (PrimNByteInteger 1)
+                return (fromInteger i)
+        
+        decomp PrimWord16
+            | supported PrimWord8 = do
+                b0 <- prompt PrimWord8
+                b1 <- prompt PrimWord8
+                return (buildWord16 b0 b1)
+            | supported PrimWord32 = do
+                w <- prompt PrimWord32
+                return (fromIntegral w)
+            | supported PrimWord64 = do
+                w <- prompt PrimWord64
+                return (fromIntegral w)
+            | supported PrimDouble = do
+                d <- prompt PrimDouble
+                return (truncate (d * 65536))
+            | supported (PrimNByteInteger 2) = do
+                i <- prompt (PrimNByteInteger 2)
+                return (fromInteger i)
+        
+        decomp PrimWord32
+            | supported PrimWord16 = do
+                w0 <- prompt PrimWord16
+                w1 <- prompt PrimWord16
+                
+                return (buildWord32' w0 w1)
+            | supported PrimWord8 = do
+                b0 <- prompt PrimWord8
+                b1 <- prompt PrimWord8
+                b2 <- prompt PrimWord8
+                b3 <- prompt PrimWord8
+                
+                return (buildWord32 b0 b1 b2 b3)
+            | supported PrimWord64 = do
+                w <- prompt PrimWord64
+                return (fromIntegral w)
+            | supported PrimDouble = do
+                d <- prompt PrimDouble
+                return (truncate (d * 4294967296))
+            | supported (PrimNByteInteger 4) = do
+                i <- prompt (PrimNByteInteger 4)
+                return (fromInteger i)
+        
+        decomp PrimWord64
+            | supported PrimWord32 = do
+                w0 <- prompt PrimWord32
+                w1 <- prompt PrimWord32
+                
+                return (buildWord64'' w0 w1)
+            | supported PrimWord16 = do
+                w0 <- prompt PrimWord16
+                w1 <- prompt PrimWord16
+                w2 <- prompt PrimWord16
+                w3 <- prompt PrimWord16
+                
+                return (buildWord64' w0 w1 w2 w3)
+            | supported PrimWord8 = do
+                b0 <- prompt PrimWord8
+                b1 <- prompt PrimWord8
+                b2 <- prompt PrimWord8
+                b3 <- prompt PrimWord8
+                b4 <- prompt PrimWord8
+                b5 <- prompt PrimWord8
+                b6 <- prompt PrimWord8
+                b7 <- prompt PrimWord8
+                
+                return (buildWord64 b0 b1 b2 b3 b4 b5 b6 b7)
+            | supported PrimDouble = do
+                -- Need 2 doubles, because a uniform [0,1) double only has
+                -- about 52 bits of reliable entropy
+                d0 <- prompt PrimDouble
+                d1 <- prompt PrimDouble
+                
+                let w0 = truncate (d0 * 4294967296)
+                    w1 = truncate (d1 * 4294967296)
+                
+                return (w0 .|. (w1 `shiftL` 32))
+            | supported (PrimNByteInteger 8) = do
+                i <- prompt (PrimNByteInteger 8)
+                return (fromInteger i)
+        
+        decomp PrimDouble = do
+            word <- decomp PrimWord64
+            return (wordToDouble word)
+        
+        decomp (PrimNByteInteger 1) = do
+            x <- decomp PrimWord8
+            return $! toInteger x
+        decomp (PrimNByteInteger 2) = do
+            x <- decomp PrimWord16
+            return $! toInteger x
+        decomp (PrimNByteInteger 4) = do
+            x <- decomp PrimWord32
+            return $! toInteger x
+        decomp (PrimNByteInteger 8) = do
+            x <- decomp PrimWord64
+            return $! toInteger x
+        decomp (PrimNByteInteger (n+8))  = do
+            x <- decomp PrimWord64
+            y <- decomp (PrimNByteInteger n)
+            return $! (toInteger x `shiftL` (n `shiftL` 3)) .|. y
+        decomp (PrimNByteInteger (n+4))  = do
+            x <- decomp PrimWord32
+            y <- decomp (PrimNByteInteger n)
+            return $! (toInteger x `shiftL` (n `shiftL` 3)) .|. y
+        decomp (PrimNByteInteger (n+2))  = do
+            x <- decomp PrimWord16
+            y <- decomp (PrimNByteInteger n)
+            return $! (toInteger x `shiftL` (n `shiftL` 3)) .|. y
+-- REDUNDANT CASE
+--        decomp (PrimNByteInteger (n+1))  = do
+--            x <- decomp PrimWord8
+--            y <- decomp (PrimNByteInteger n)
+--            return $! (toInteger x `shiftL` (n `shiftL` 3)) .|. y
+        decomp (PrimNByteInteger _) = return 0
+        
+        decomp _ = error ("decomposePrimWhere: no supported primitive to satisfy " ++ show requested)
diff --git a/src/Data/Random/Internal/TH.hs b/src/Data/Random/Internal/TH.hs
--- a/src/Data/Random/Internal/TH.hs
+++ b/src/Data/Random/Internal/TH.hs
@@ -26,13 +26,14 @@
 
 import Data.Word
 import Data.Int
+import Control.Monad
 
 -- |Names of standard 'Integral' types
 integralTypes :: [Name]
 integralTypes = 
-    [ ''Int,   ''Integer
-    , ''Int8,  ''Int16,  ''Int32,  ''Int64
-    , ''Word8, ''Word16, ''Word32, ''Word64
+    [ ''Integer
+    , ''Int,  ''Int8,  ''Int16,  ''Int32,  ''Int64
+    , ''Word, ''Word8, ''Word16, ''Word32, ''Word64
     ]
 
 -- |Names of standard 'RealFloat' types
@@ -69,11 +70,13 @@
 -- This code takes those 2 instance declarations and creates identical ones for
 -- every type named in 'integralTypes'.
 replicateInstances :: (Monad m, Data t) => Name -> [Name] -> m [t] -> m [t]
-replicateInstances standin types decls = do
-    decls <- decls
-    sequence
-        [ everywhereM (mkM (return . replaceName standin t)) dec
-        | t <- types
-        , dec <- decls
-        ]
+replicateInstances standin types getDecls = liftM concat $ sequence
+    [ do
+        decls <- getDecls
+        sequence
+            [ everywhereM (mkM (return . replaceName standin t)) dec
+            | dec <- decls
+            ]
+    | t <- types
+    ]
 
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
@@ -6,19 +6,55 @@
 module Data.Random.Internal.Words where
 
 import Foreign
-import GHC.IOBase
 
-import Data.Bits
-import Data.Word
+-- TODO: add a build flag for endianness-invariance, or just find a way
+-- to make sure these operations all do the right thing without costing 
+-- anything extra at runtime
 
-{-# INLINE buildWord #-}
+{-# INLINE buildWord16 #-}
 -- |Build a word out of 8 bytes.  No promises are made regarding the order
 -- in which the bytes are stuffed.  Note that this means that a 'RandomSource'
 -- or 'MonadRandom' making use of the default definition of 'getRandomWord', etc.,
 -- may return different random values on different platforms when started 
 -- with the same seed, depending on the platform's endianness.
-buildWord :: Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word64
-buildWord b0 b1 b2 b3 b4 b5 b6 b7
+buildWord16 :: Word8 -> Word8 -> Word16
+buildWord16 b0 b1
+    = unsafePerformIO . allocaBytes 2 $ \p -> do
+        pokeByteOff p 0 b0
+        pokeByteOff p 1 b1
+        peek (castPtr p)
+
+{-# INLINE buildWord32 #-}
+-- |Build a word out of 8 bytes.  No promises are made regarding the order
+-- in which the bytes are stuffed.  Note that this means that a 'RandomSource'
+-- or 'MonadRandom' making use of the default definition of 'getRandomWord', etc.,
+-- may return different random values on different platforms when started 
+-- with the same seed, depending on the platform's endianness.
+buildWord32 :: Word8 -> Word8 -> Word8 -> Word8 -> Word32
+buildWord32 b0 b1 b2 b3
+    = unsafePerformIO . allocaBytes 4 $ \p -> do
+        pokeByteOff p 0 b0
+        pokeByteOff p 1 b1
+        pokeByteOff p 2 b2
+        pokeByteOff p 3 b3
+        peek (castPtr p)
+
+{-# INLINE buildWord32' #-}
+buildWord32' :: Word16 -> Word16 -> Word32
+buildWord32' w0 w1
+    = unsafePerformIO . allocaBytes 4 $ \p -> do
+        pokeByteOff p 0 w0
+        pokeByteOff p 2 w1
+        peek (castPtr p)
+
+{-# INLINE buildWord64 #-}
+-- |Build a word out of 8 bytes.  No promises are made regarding the order
+-- in which the bytes are stuffed.  Note that this means that a 'RandomSource'
+-- or 'MonadRandom' making use of the default definition of 'getRandomWord', etc.,
+-- may return different random values on different platforms when started 
+-- with the same seed, depending on the platform's endianness.
+buildWord64 :: Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word8 -> Word64
+buildWord64 b0 b1 b2 b3 b4 b5 b6 b7
     = unsafePerformIO . allocaBytes 8 $ \p -> do
         pokeByteOff p 0 b0
         pokeByteOff p 1 b1
@@ -30,6 +66,36 @@
         pokeByteOff p 7 b7
         peek (castPtr p)
 
+{-# INLINE buildWord64' #-}
+buildWord64' :: Word16 -> Word16 -> Word16 -> Word16 -> Word64
+buildWord64' w0 w1 w2 w3
+    = unsafePerformIO . allocaBytes 8 $ \p -> do
+        pokeByteOff p 0 w0
+        pokeByteOff p 2 w1
+        pokeByteOff p 4 w2
+        pokeByteOff p 6 w3
+        peek (castPtr p)
+
+{-# INLINE buildWord64'' #-}
+buildWord64'' :: Word32 -> Word32 -> Word64
+buildWord64'' w0 w1
+    = unsafePerformIO . allocaBytes 8 $ \p -> do
+        pokeByteOff p 0 w0
+        pokeByteOff p 4 w1
+        peek (castPtr p)
+
+{-# INLINE word32ToFloat #-}
+-- |Pack the low 23 bits from a 'Word32' into a 'Float' in the range [0,1).
+-- Used to convert a 'stdUniform' 'Word32' to a 'stdUniform' 'Double'.
+word32ToFloat :: Word32 -> Float
+word32ToFloat x = (encodeFloat $! toInteger (x .&. 0x007fffff {- 2^23-1 -} )) $ (-23)
+
+{-# INLINE word32ToFloatWithExcess #-}
+-- |Same as word32ToFloat, but also return the unused bits (as the 9
+-- least significant bits of a 'Word32')
+word32ToFloatWithExcess :: Word32 -> (Float, Word32)
+word32ToFloatWithExcess x = (word32ToFloat x, x `shiftR` 23)
+
 {-# INLINE wordToFloat #-}
 -- |Pack the low 23 bits from a 'Word64' into a 'Float' in the range [0,1).
 -- Used to convert a 'stdUniform' 'Word64' to a 'stdUniform' 'Double'.
@@ -47,6 +113,12 @@
 -- Used to convert a 'stdUniform' 'Word64' to a 'stdUniform' 'Double'.
 wordToDouble :: Word64 -> Double
 wordToDouble x = (encodeFloat $! toInteger (x .&. 0x000fffffffffffff {- 2^52-1 -})) $ (-52)
+
+{-# INLINE word32ToDouble #-}
+-- |Pack a 'Word32' into a 'Double' in the range [0,1).  Note that a Double's 
+-- mantissa is 52 bits, so this does not fill all of them.
+word32ToDouble :: Word32 -> Double
+word32ToDouble x = (encodeFloat $! toInteger x) $ (-32)
 
 {-# INLINE wordToDoubleWithExcess #-}
 -- |Same as wordToDouble, but also return the unused bits (as the 12
diff --git a/src/Data/Random/Lift.hs b/src/Data/Random/Lift.hs
--- a/src/Data/Random/Lift.hs
+++ b/src/Data/Random/Lift.hs
@@ -13,6 +13,12 @@
 -- For instances where 'm' and 'n' have 'return'/'pure' defined,
 -- these instances must satisfy
 -- @lift (return x) == return x@.
+-- 
+-- This form of 'lift' has an extremely general type and is used primarily to
+-- support 'sample'.  Its excessive generality is the main reason it's not
+-- exported from "Data.Random".  'RVarT' is, however, an instance of 
+-- 'T.MonadTrans', which in most cases is the preferred way
+-- to do the lifting.
 class Lift m n where
     lift :: m a -> n a
 
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
@@ -24,16 +24,20 @@
     
     return (SRS.shuffle xs (reverse is))
 
--- | A random variable that shuffles a list of a known length. Useful for 
--- shuffling large lists when the length is known in advance.
--- Avoids needing to traverse the list to discover its length.  Each ordering
--- has equal probability.
---
--- Throws an error the list is not exactly as long as claimed.
+-- | A random variable that shuffles a list of a known length (or a list
+-- prefix of the specified length). Useful for shuffling large lists when 
+-- the length is known in advance.  Avoids needing to traverse the list to
+-- discover its length.  Each ordering has equal probability.
 shuffleN :: Int -> [a] -> RVar [a]
-shuffleN 0 xs = return []
-shuffleN m@(n+1) xs = do
-    is <- sequence [uniform 0 i | i <- [n,n-1..1]]
-    return (SRS.shuffle xs is)
+shuffleN n xs = shuffleNofM n n xs
 
-    
+-- | A random variable that selects N arbitrary elements of a list of known length M.
+shuffleNofM :: Int -> Int -> [a] -> RVar [a]
+shuffleNofM 0 _ _ = return []
+shuffleNofM n m xs
+    | n > m     = error "shuffleNofM: n > m"
+    | otherwise = do
+        is <- sequence [uniform 0 i | i <- take n [m-1, m-2 ..1]]
+        return (take n $ SRS.shuffle (take m xs) is)
+shuffleNofM _ _ _ = error "shuffleNofM: negative length specified"
+
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
@@ -18,68 +18,183 @@
     , runRVar
     , RVarT
     , runRVarT
-    , nByteInteger
-    , nBitInteger
+    , runRVarTWith
     ) where
 
 
+import Data.Random.Internal.Primitives
 import Data.Random.Source
 import Data.Random.Lift as L
 
-import Data.Bits
-
 import qualified Control.Monad.Trans as T
 import Control.Applicative
-import Control.Monad.Reader
 import Control.Monad.Identity
+import Control.Monad.Prompt (PromptT, runPromptT, prompt)
 
--- |An opaque type containing a \"random variable\" - a value 
--- which depends on the outcome of some random process.
+-- |An opaque type modeling a \"random variable\" - a value 
+-- which depends on the outcome of some random event.  'RVar's 
+-- can be conveniently defined by an imperative-looking style:
+-- 
+-- > normalPair =  do
+-- >     u <- stdUniform
+-- >     t <- stdUniform
+-- >     let r = sqrt (-2 * log u)
+-- >         theta = (2 * pi) * t
+-- >         
+-- >         x = r * cos theta
+-- >         y = r * sin theta
+-- >     return (x,y)
+-- 
+-- OR by a more applicative style:
+-- 
+-- > logNormal = exp <$> stdNormal
+--
+-- Once defined (in any style), there are a couple ways to sample 'RVar's:
+-- 
+-- * In a monad, using a 'RandomSource':
+-- 
+-- > sampleFrom DevRandom (uniform 1 100) :: IO Int
+-- 
+-- * As a pure function transforming a functional RNG:
+-- 
+-- > sampleState (uniform 1 100) :: StdGen -> (Int, StdGen)
 type RVar = RVarT Identity
 
--- | single combined container allowing all the relevant 
--- dictionaries (plus the RandomSource item itself) to be passed
--- with one pointer.
-data RVarDict n m where
-    RVarDict :: (Lift n m, Monad m, RandomSource m s) => s -> RVarDict n m
-
+-- |\"Run\" an 'RVar' - samples the random variable from the provided
+-- source of entropy.
 runRVar :: RandomSource m s => RVar a -> s -> m a
 runRVar = runRVarT
 
 -- |A random variable with access to operations in an underlying monad.  Useful
 -- examples include any form of state for implementing random processes with hysteresis,
 -- or writer monads for implementing tracing of complicated algorithms.
-newtype RVarT n a = RVarT { unRVarT :: forall m r. (a -> m r) -> RVarDict n m -> m r }
+-- 
+-- For example, a simple random walk can be implemented as an 'RVarT' 'IO' value:
+--
+-- > rwalkIO :: IO (RVarT IO Double)
+-- > rwalkIO d = do
+-- >     lastVal <- newIORef 0
+-- >     
+-- >     let x = do
+-- >             prev    <- lift (readIORef lastVal)
+-- >             change  <- rvarT StdNormal
+-- >             
+-- >             let new = prev + change
+-- >             lift (writeIORef lastVal new)
+-- >             return new
+-- >         
+-- >     return x
+--
+-- To run the random walk, it must first be initialized, and then it can be sampled as usual:
+--
+-- > do
+-- >     rw <- rwalkIO
+-- >     x <- sampleFrom DevURandom rw
+-- >     y <- sampleFrom DevURandom rw
+-- >     ...
+--
+-- The same random-walk process as above can be implemented using MTL types
+-- as follows (using @import Control.Monad.Trans as MTL@):
+-- 
+-- > rwalkState :: RVarT (State Double) Double
+-- > rwalkState = do
+-- >     prev <- MTL.lift get
+-- >     change  <- rvarT StdNormal
+-- >     
+-- >     let new = prev + change
+-- >     MTL.lift (put new)
+-- >     return new
+-- 
+-- Invocation is straightforward (although a bit noisy) if you're used 
+-- to MTL, but there is a gotcha lurking here: @sample@ and 'runRVarT' 
+-- inherit the extreme generality of 'lift', so there will almost always
+-- need to be an explicit type signature lurking somewhere in any client 
+-- code making use of 'RVarT' with MTL types.  In this example, the 
+-- inferred type of @start@ would be too general to be practical, so the
+-- signature for @rwalk@  explicitly fixes it to 'Double'.  Alternatively, 
+-- in this case @sample@ could be replaced with
+-- @\\x -> runRVarTWith MTL.lift x StdRandom@.
+-- 
+-- > rwalk :: Int -> Double -> StdGen -> ([Double], StdGen)
+-- > rwalk count start gen = evalState (runStateT (sample (replicateM count rwalkState)) gen) start
+newtype RVarT m a = RVarT { unRVarT :: PromptT Prim m a }
 
--- | \"Runs\" the monad.
+-- | \"Runs\" an 'RVarT', sampling the random variable it defines.
+-- 
+-- The 'Lift' context allows random variables to be defined using a minimal
+-- underlying functor ('Identity' is sufficient for \"conventional\" random
+-- variables) and then sampled in any monad into which the underlying functor 
+-- can be embedded (which, for 'Identity', is all monads).
+-- 
+-- The lifting is very important - without it, every 'RVar' would have
+-- to either be given access to the full capability of the monad in which it
+-- will eventually be sampled (which, incidentally, would also have to be 
+-- monomorphic so you couldn't sample one 'RVar' in more than one monad)
+-- or functions manipulating 'RVar's would have to use higher-ranked 
+-- types to enforce the same kind of isolation and polymorphism.
+-- 
+-- For non-standard liftings or those where you would rather not introduce a
+-- 'Lift' instance, see 'runRVarTWith'.
+{-# INLINE runRVarT #-}
 runRVarT :: (Lift n m, RandomSource m s) => RVarT n a -> s -> m a
-runRVarT (RVarT m) (src) = m return (RVarDict src)
+runRVarT (RVarT m) src = runPromptT return bindP bindN m
+    where
+        bindP prim cont = getRandomPrimFrom src prim >>= cont
+        bindN nExp cont = lift nExp >>= cont
 
+-- |Like 'runRVarT' but allowing a user-specified lift operation.  This 
+-- operation must obey the \"monad transformer\" laws:
+--
+-- > lift . return = return
+-- > lift (x >>= f) = (lift x) >>= (lift . f)
+--
+-- One example of a useful non-standard lifting would be one that takes @State s@ to
+-- another monad with a different state representation (such as @IO@ with the
+-- state mapped to an @IORef@):
+--
+-- > embedState :: (Monad m) => m s -> (s -> m ()) -> State s a -> m a
+-- > embedState get put = \m -> do
+-- >     s <- get
+-- >     (res,s) <- return (runState m s)
+-- >     put s
+-- >     return res
+{-# INLINE runRVarTWith #-}
+runRVarTWith :: (RandomSource m s) => (forall t. n t -> m t) -> RVarT n a -> s -> m a
+runRVarTWith liftN (RVarT m) src = runPromptT return bindP bindN m
+    where
+        bindP prim cont = getRandomPrimFrom src prim >>= cont
+        bindN nExp cont = liftN nExp >>= cont
+
 instance Functor (RVarT n) where
     fmap = liftM
 
 instance Monad (RVarT n) where
-    return x = RVarT $ \k _ -> k x
-    fail s   = RVarT $ \_ (RVarDict _) -> fail s
-    (RVarT m) >>= k = RVarT $ \c s -> m (\a -> unRVarT (k a) c s) s
+    return x = RVarT (return $! x)
+    fail s   = RVarT (fail s)
+    (RVarT m) >>= k = RVarT (m >>= \x -> x `seq` unRVarT (k x))
 
 instance Applicative (RVarT n) where
     pure  = return
     (<*>) = ap
 
 instance T.MonadTrans RVarT where
-    lift m = RVarT $ \k r@(RVarDict _) -> L.lift m >>= \a -> k a
+    lift m = RVarT (T.lift m)
 
 instance Lift (RVarT Identity) (RVarT m) where
-    lift (RVarT m) = RVarT $ \k (RVarDict src) -> m k (RVarDict src)
+    lift (RVarT m) = RVarT (runPromptT return bindP bindN m)
+        where
+            bindP prim  cont = prompt prim >>= cont
+            bindN idExp cont = cont (runIdentity idExp)
 
-instance MonadIO m => MonadIO (RVarT m) where
-    liftIO = T.lift . liftIO
+instance T.MonadIO m => T.MonadIO (RVarT m) where
+    liftIO = T.lift . T.liftIO
 
 instance MonadRandom (RVarT n) where
-    getRandomByte   = RVarT $ \k (RVarDict s) -> getRandomByteFrom   s >>= \a -> k a
-    getRandomWord   = RVarT $ \k (RVarDict s) -> getRandomWordFrom   s >>= \a -> k a
-    getRandomDouble = RVarT $ \k (RVarDict s) -> getRandomDoubleFrom s >>= \a -> k a
+    supportedPrims _ _ = True
+    {-# INLINE getSupportedRandomPrim #-}
+    getSupportedRandomPrim p    = RVarT (prompt p)
+    {-# INLINE getRandomPrim #-}
+    getRandomPrim p = RVarT (prompt p)
 
 -- I would really like to be able to do this, but I can't because of the
 -- blasted Eq and Show in Num's class context...
@@ -91,32 +206,3 @@
 --     signum = liftA signum
 --     abs = liftA abs
 --     fromInteger = pure . fromInteger
-
--- some 'fundamental' RVarTs
--- this maybe ought to even be a part of the RandomSource class...
-{-# INLINE nByteInteger #-}
--- |A random variable evenly distributed over all unsigned integers from
--- 0 to 2^(8*n)-1, inclusive.
-nByteInteger :: Int -> RVarT m Integer
-nByteInteger 1 = do
-    x <- getRandomByte
-    return $! toInteger x
-nByteInteger 8 = do
-    x <- getRandomWord
-    return $! toInteger x
-nByteInteger n = nBitInteger (n `shiftL` 3)
-
-{-# INLINE nBitInteger #-}
--- |A random variable evenly distributed over all unsigned integers from
--- 0 to 2^n-1, inclusive.
-nBitInteger :: Int -> RVarT m Integer
-nBitInteger 8  = do
-    x <- getRandomByte
-    return $! toInteger x
-nBitInteger (n+64) = do
-    x <- getRandomWord
-    y <- nBitInteger n
-    return $! (toInteger x `shiftL` n) .|. y
-nBitInteger n = do
-        x <- getRandomWord
-        return $! toInteger (x `shiftR` (64-n))
diff --git a/src/Data/Random/Sample.hs b/src/Data/Random/Sample.hs
--- a/src/Data/Random/Sample.hs
+++ b/src/Data/Random/Sample.hs
@@ -9,6 +9,7 @@
 
 module Data.Random.Sample where
 
+import Control.Monad.State 
 import Data.Random.Distribution
 import Data.Random.Lift
 import Data.Random.RVar
@@ -17,7 +18,8 @@
 
 -- |A typeclass allowing 'Distribution's and 'RVar's to be sampled.  Both may
 -- also be sampled via 'runRVar' or 'runRVarT', but I find it psychologically
--- pleasing to be able to sample both using this function.
+-- pleasing to be able to sample both using this function, as they are two
+-- separate abstractions for one base concept: a random variable.
 class Sampleable d m t where
     -- |Directly sample from a distribution or random variable, using the given source of entropy.
     sampleFrom :: RandomSource m s => s -> d t -> m t
@@ -29,7 +31,17 @@
 instance Lift m n => Sampleable (RVarT m) n t where
     sampleFrom src x = runRVarT x src
 
--- |Sample a distribution using the default source of entropy for the
+-- |Sample a random variable using the default source of entropy for the
 -- monad in which the sampling occurs.
 sample :: (Sampleable d m t, MonadRandom m) => d t -> m t
 sample = sampleFrom StdRandom
+
+-- |Sample a random variable in a \"functional\" style.  Typical instantiations
+-- of @s@ are @System.Random.StdGen@ or @System.Random.Mersenne.Pure64.PureMT@.
+sampleState :: (Sampleable d (State s) t, MonadRandom (State s)) => d t -> s -> (t, s)
+sampleState thing = runState (sample thing)
+
+-- |Sample a random variable in a \"semi-functional\" style.  Typical instantiations
+-- of @s@ are @System.Random.StdGen@ or @System.Random.Mersenne.Pure64.PureMT@.
+sampleStateT :: (Sampleable d (StateT s m) t, MonadRandom (StateT s m)) => d t -> s -> m (t, s)
+sampleStateT thing = runStateT (sample thing)
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
@@ -2,18 +2,20 @@
  -      ``Data/Random/Source''
  -}
 {-# LANGUAGE
-    MultiParamTypeClasses, FlexibleInstances
+    MultiParamTypeClasses, FlexibleInstances, GADTs
   #-}
 
 module Data.Random.Source
     ( MonadRandom(..)
     , RandomSource(..)
+    , Prim(..)
     ) where
 
 import Data.Word
-import Control.Monad
+import Control.Monad.Prompt
+import Data.Tagged
 
-import Data.Random.Internal.Words
+import Data.Random.Internal.Primitives
 
 -- |A typeclass for monads with a chosen source of entropy.  For example,
 -- 'RVar' is such a monad - the source from which it is (eventually) sampled
@@ -21,68 +23,119 @@
 -- when directly requesting entropy for a random variable these functions
 -- are used.
 -- 
--- The minimal definition is either 'getRandomByte' or 'getRandomWord'.
--- 'getRandomDouble' is defaulted in terms of 'getRandomWord'.
+-- The minimal definition is 'supportedPrims' and 'getSupportedRandomPrim'
+-- with cases for those primitives where 'supportedPrims' returns 'True'.
+--
+-- It is recommended (despite the warnings it generates) that, even when
+-- all primitives are supported, a final wildcard case of 'supportedPrims' is
+-- specified, as:
+-- 
+-- > supportedPrims _ _ = False
+--
+-- The overlapping pattern warnings can be suppressed (without suppressing 
+-- other, genuine, overlapping-pattern warnings) by the GHC flag
+-- @-fno-warn-simple-patterns@.  This is not actually the documented behavior
+-- of that flag as far as I can find in 3 google-minutes, but it works with
+-- GHC 6.12.1 anyway, and that's good enough for me.
+--
+-- Note that it is very important that at least 'supportedPrims' (and preferably
+-- 'getSupportedRandomPrim' as well) gets inlined into the default implementation
+-- of 'getRandomPrim'.  If your 'supportedPrims' is more than about 2 or 3
+-- cases, add an INLINE pragma so that it can be optimized out of 'getRandomPrim'.
 class Monad m => MonadRandom m where
-    -- |Get a random uniformly-distributed byte.
-    getRandomByte :: m Word8
-    getRandomByte = do
-        word <- getRandomWord
-        return (fromIntegral word)
-    
-    -- |Get a random 'Word64' uniformly-distributed over the full range of the type.
-    getRandomWord :: m Word64
-    getRandomWord = do
-        b0 <- getRandomByte
-        b1 <- getRandomByte
-        b2 <- getRandomByte
-        b3 <- getRandomByte
-        b4 <- getRandomByte
-        b5 <- getRandomByte
-        b6 <- getRandomByte
-        b7 <- getRandomByte
-        
-        return (buildWord b0 b1 b2 b3 b4 b5 b6 b7)
+    -- |Predicate indicating whether a given primitive is supported by the
+    -- instance.  The first parameter is a phantom used to select the instance.
+    supportedPrims :: m () -> Prim t -> Bool
     
-    -- |Get a random 'Double' uniformly-distributed over the interval [0,1)
-    getRandomDouble :: m Double
-    getRandomDouble = do
-        word <- getRandomWord
-        return (wordToDouble word)
+    -- |Generate a random value corresponding to the specified primitive.  Will
+    -- not be called unless supportedPrims returns true for that primitive.
+    getSupportedRandomPrim :: Prim t -> m t
 
+    -- This could just be a function, but placing it in a dictionary gives
+    -- GHC a place to optimize it separately for each instance, which is 
+    -- kinda the whole point of the 'Prim' machinery:
+    -- 
+    -- |Generate a random value corresponding to the specified primitive.  The
+    -- default implementation makes use of 'supportedPrims' and 'getSupportedRandomPrim'
+    -- to construct any required Prim out of the supported ones.
+    {-# NOINLINE getRandomPrim #-}
+    getRandomPrim :: Prim t -> m t
+    getRandomPrim prim = val
+        where
+            val = runPromptM getSupportedRandomPrim (decomposePrimWhere (supportedPrims mPhantom) prim)
+            mPhantom = error "supportedPrims tried to evaluate a phantom parameter" `asTypeOf` (val >> return ())
+
 -- |A source of entropy which can be used in the given monad.
 --
--- The minimal definition is either 'getRandomByteFrom' or 'getRandomWordFrom'.
--- 'getRandomDoubleFrom' is defaulted in terms of 'getRandomWordFrom'
+-- The minimal definition is 'supportedPrimsFrom' and 'getSupportedRandomPrimFrom'
+-- with cases for those primitives where 'supportedPrimsFrom' returns 'True'.
+-- 
+-- Note that it is very important that at least 'supportedPrimsFrom' (and preferably
+-- 'getSupportedRandomPrimFrom' as well) gets inlined into the default implementation
+-- of 'getRandomPrimFrom'.  If your 'supportedPrimsFrom' is more than about 2 or 3
+-- cases, add an INLINE pragma so that it can be optimized out of 'getRandomPrimFrom'.
+-- 
+-- See also 'MonadRandom'.
 class Monad m => RandomSource m s where
-    -- |Get a random uniformly-distributed byte.
-    getRandomByteFrom :: s -> m Word8
-    getRandomByteFrom src = do
-        word <- getRandomWordFrom src
-        return (fromIntegral word)
+    -- |Predicate indicating whether a given primitive is supported by the
+    -- instance.  The tag on the first parameter is a phantom used only to
+    -- select the instance, but the value itself may be inspected.
+    supportedPrimsFrom :: Tagged (m ()) s -> Prim t -> Bool
     
-    -- |Get a random 'Word64' uniformly-distributed over the full range of the type.
-    getRandomWordFrom :: s -> m Word64
-    getRandomWordFrom src = do
-        b0 <- getRandomByteFrom src
-        b1 <- getRandomByteFrom src
-        b2 <- getRandomByteFrom src
-        b3 <- getRandomByteFrom src
-        b4 <- getRandomByteFrom src
-        b5 <- getRandomByteFrom src
-        b6 <- getRandomByteFrom src
-        b7 <- getRandomByteFrom src
-        
-        return (buildWord b0 b1 b2 b3 b4 b5 b6 b7)
+    -- |Generate a random value corresponding to the specified primitive
+    getSupportedRandomPrimFrom :: s -> Prim t -> m t
     
-    -- |Get a random 'Double' uniformly-distributed over the interval [0,1)
-    getRandomDoubleFrom :: s -> m Double
-    getRandomDoubleFrom src = do
-        word <- getRandomWordFrom src
-        return (wordToDouble word)
+    
+    -- This could just be a function, but placing it in a dictionary gives
+    -- GHC a place to optimize it separately for each instance, which is 
+    -- kinda the whole point of the 'Prim' machinery:
+    -- 
+    -- |Generate a random value corresponding to the specified primitive.  The
+    -- default implementation makes use of 'supportedPrimsFrom' and
+    -- 'getSupportedRandomPrimFrom' to construct any required Prim out of 
+    -- the supported ones.
+    {-# NOINLINE getRandomPrimFrom #-}
+    getRandomPrimFrom :: s -> Prim t -> m t
+    getRandomPrimFrom src prim = val
+        where
+            val = runPromptM (getSupportedRandomPrimFrom src) (decomposePrimWhere supported prim)
+            supported :: Prim t -> Bool
+            supported = supportedPrimsFrom (tagIt (val >> return ()) src)
+            
+            tagIt :: a -> b -> Tagged a b
+            tagIt _ it = Tagged it
 
 instance Monad m => RandomSource m (m Word8) where
-    getRandomByteFrom = id
+    supportedPrimsFrom _ PrimWord8 = True
+    supportedPrimsFrom _ _ = False
+    
+    getSupportedRandomPrimFrom f PrimWord8 = f
+    getSupportedRandomPrimFrom _ p = error ("getSupportedRandomPrimFrom/RandomSource m (m Word8): unsupported prim requested: " ++ show p)
 
+instance Monad m => RandomSource m (m Word16) where
+    supportedPrimsFrom _ PrimWord16 = True
+    supportedPrimsFrom _ _ = False
+    
+    getSupportedRandomPrimFrom f PrimWord16 = f
+    getSupportedRandomPrimFrom _ p = error ("getSupportedRandomPrimFrom/RandomSource m (m Word16): unsupported prim requested: " ++ show p)
+
+instance Monad m => RandomSource m (m Word32) where
+    supportedPrimsFrom _ PrimWord32 = True
+    supportedPrimsFrom _ _ = False
+    
+    getSupportedRandomPrimFrom f PrimWord32 = f
+    getSupportedRandomPrimFrom _ p = error ("getSupportedRandomPrimFrom/RandomSource m (m Word32): unsupported prim requested: " ++ show p)
+
 instance Monad m => RandomSource m (m Word64) where
-    getRandomWordFrom = id
+    supportedPrimsFrom _ PrimWord64 = True
+    supportedPrimsFrom _ _ = False
+    
+    getSupportedRandomPrimFrom f PrimWord64 = f
+    getSupportedRandomPrimFrom _ p = error ("getSupportedRandomPrimFrom/RandomSource m (m Word64): unsupported prim requested: " ++ show p)
+
+instance Monad m => RandomSource m (m Double) where
+    supportedPrimsFrom _ PrimDouble = True
+    supportedPrimsFrom _ _ = False
+    
+    getSupportedRandomPrimFrom f PrimDouble = f
+    getSupportedRandomPrimFrom _ p = error ("getSupportedRandomPrimFrom/RandomSource m (m Double): unsupported prim requested: " ++ show p)
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
@@ -2,7 +2,7 @@
  -      ``Data/Random/Source/DevRandom''
  -}
 {-# LANGUAGE
-    MultiParamTypeClasses
+    MultiParamTypeClasses, GADTs
   #-}
 
 module Data.Random.Source.DevRandom 
@@ -11,7 +11,7 @@
 
 import Data.Random.Source
 
-import System.IO (openBinaryFile, hGetBuf, IOMode(..))
+import System.IO (openBinaryFile, hGetBuf, Handle, IOMode(..))
 import Foreign
 
 -- |On systems that have it, \/dev\/random is a handy-dandy ready-to-use source
@@ -22,19 +22,32 @@
 -- purposes other than cryptography, \/dev\/urandom is preferable because when it
 -- runs out of real entropy it'll still churn out pseudorandom data.
 data DevRandom = DevRandom | DevURandom
+    deriving (Eq, Show)
 
 {-# NOINLINE devRandom  #-}
+devRandom :: Handle
 devRandom  = unsafePerformIO (openBinaryFile "/dev/random"  ReadMode)
 {-# NOINLINE devURandom #-}
+devURandom :: Handle
 devURandom = unsafePerformIO (openBinaryFile "/dev/urandom" ReadMode)
 
+dev :: DevRandom -> Handle
 dev DevRandom  = devRandom
 dev DevURandom = devURandom
 
 instance RandomSource IO DevRandom where
-    getRandomByteFrom src  = allocaBytes 1 $ \buf -> do
+    supportedPrimsFrom _ PrimWord8          = True
+    supportedPrimsFrom _ PrimWord32         = True
+    supportedPrimsFrom _ PrimWord64         = True
+    supportedPrimsFrom _ _ = False
+    
+    getSupportedRandomPrimFrom src PrimWord8  = allocaBytes 1 $ \buf -> do
         1 <- hGetBuf (dev src) buf  1
         peek buf
-    getRandomWordFrom src  = allocaBytes 8 $ \buf -> do
+    getSupportedRandomPrimFrom src PrimWord32  = allocaBytes 1 $ \buf -> do
+        4 <- hGetBuf (dev src) buf  4
+        peek (castPtr buf)
+    getSupportedRandomPrimFrom src PrimWord64  = allocaBytes 8 $ \buf -> do
         8 <- hGetBuf (dev src) buf  8
         peek (castPtr buf)
+    getSupportedRandomPrimFrom src prim = error ("getSupportedRandomPrimFrom/" ++ show src ++ ": unsupported prim requested: " ++ show prim)
diff --git a/src/Data/Random/Source/MWC.hs b/src/Data/Random/Source/MWC.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Random/Source/MWC.hs
@@ -0,0 +1,46 @@
+{-# LANGUAGE
+        MultiParamTypeClasses,
+        FlexibleInstances,
+        GADTs
+  #-}
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+-- |This module defines the following instances:
+-- 
+-- > instance RandomSource (ST s) (Gen s)
+-- > instance RandomSource IO (Gen RealWorld)
+module Data.Random.Source.MWC where
+
+import Data.Random.Internal.Primitives
+import Data.Random.Internal.Words
+import Data.Random.Source
+import System.Random.MWC
+import Control.Monad.ST
+
+instance RandomSource (ST s) (Gen s) where
+    {-# INLINE supportedPrimsFrom #-}
+    supportedPrimsFrom _ PrimWord8  = True
+    supportedPrimsFrom _ PrimWord16 = True
+    supportedPrimsFrom _ PrimWord32 = True
+    supportedPrimsFrom _ PrimWord64 = True
+    supportedPrimsFrom _ PrimDouble = True
+    supportedPrimsFrom _ _ = False
+    
+    {-# INLINE getSupportedRandomPrimFrom #-}
+    getSupportedRandomPrimFrom gen PrimWord8    = uniform gen
+    getSupportedRandomPrimFrom gen PrimWord16   = uniform gen
+    getSupportedRandomPrimFrom gen PrimWord32   = uniform gen
+    getSupportedRandomPrimFrom gen PrimWord64   = uniform gen
+    getSupportedRandomPrimFrom gen PrimDouble   = fmap wordToDouble (uniform gen)
+    getSupportedRandomPrimFrom _ p = error ("getSupportedRandomPrimFrom/Gen s: unsupported prim requested: " ++ show p)
+
+instance RandomSource IO (Gen RealWorld) where
+    {-# INLINE supportedPrimsFrom #-}
+    supportedPrimsFrom _ PrimWord8  = True
+    supportedPrimsFrom _ PrimWord16 = True
+    supportedPrimsFrom _ PrimWord32 = True
+    supportedPrimsFrom _ PrimWord64 = True
+    supportedPrimsFrom _ PrimDouble = True
+    supportedPrimsFrom _ _ = False
+    
+    {-# INLINE getSupportedRandomPrimFrom #-}
+    getSupportedRandomPrimFrom gen prim = stToIO (getSupportedRandomPrimFrom gen prim)
diff --git a/src/Data/Random/Source/PureMT.hs b/src/Data/Random/Source/PureMT.hs
--- a/src/Data/Random/Source/PureMT.hs
+++ b/src/Data/Random/Source/PureMT.hs
@@ -1,8 +1,11 @@
 {-# LANGUAGE
+    BangPatterns,
     MultiParamTypeClasses,
     FlexibleContexts, FlexibleInstances,
-    UndecidableInstances
+    UndecidableInstances,
+    GADTs
   #-}
+{-# OPTIONS_GHC -fno-warn-orphans #-}
 
 -- |This module provides functions useful for implementing new 'MonadRandom'
 -- and 'RandomSource' instances for state-abstractions containing 'PureMT'
@@ -11,52 +14,38 @@
 -- cases.
 module Data.Random.Source.PureMT where
 
-import Data.Random.Internal.Words
+import Data.Random.Internal.Primitives
 import Data.Random.Source
 import System.Random.Mersenne.Pure64
 
 import Data.StateRef
-import Data.Word
 
+import Control.Monad.Prompt
 import Control.Monad.State
 import qualified Control.Monad.ST.Strict as S
 import qualified Control.Monad.State.Strict as S
 
--- |Given a mutable reference to a 'PureMT' generator, we can make a
--- 'RandomSource' usable in any monad in which the reference can be modified.
---
--- For example, if @x :: TVar PureMT@, @getRandomWordFromMTRef x@ can be
--- used as a 'RandomSource' in 'IO', 'STM', or any monad which is an instance
--- of 'MonadIO'.  These functions can also be used to implement additional
--- 'RandomSource' instances for mutable references to 'PureMT' states.
-getRandomWordFromMTRef :: ModifyRef sr m PureMT => sr -> m Word64
-getRandomWordFromMTRef ref = do
-    atomicModifyReference ref (swap . randomWord64)
-    
-    where
-        swap (a,b) = (b,a)
-
-getRandomByteFromMTRef :: (Monad m, ModifyRef sr m PureMT) => sr -> m Word8
-getRandomByteFromMTRef ref = do
-    x <- atomicModifyReference ref (swap . randomInt)
-    return (fromIntegral x)
-    
-    where
-        swap (a,b) = (b,a)
-
--- for whatever reason, my simple wordToDouble is faster than whatever
--- the mersenne random library is using, at least in the version I have.
--- if this changes, switch to the commented version.
--- Same thing below, in getRandomDoubleFromMTState.
-getRandomDoubleFromMTRef :: (Monad m, ModifyRef sr m PureMT) => sr -> m Double
-getRandomDoubleFromMTRef src = liftM wordToDouble (getRandomWordFromMTRef src)
--- getRandomDoubleFromMTRef ref = do
---     atomicModifyReference ref (swap . randomDouble)
---     
---     where
---         swap (a,b) = (b,a)
+-- |Given a mutable reference to a 'PureMT' generator, we can implement
+-- 'RandomSource' for in any monad in which the reference can be modified.
+getRandomPrimFromMTRef :: (Monad m, ModifyRef sr m PureMT) => sr -> Prim a -> m a
+getRandomPrimFromMTRef ref prim
+    | supported prim = getThing (genPrim prim)
+    | otherwise = runPromptM (getRandomPrimFromMTRef ref) (decomposePrimWhere supported prim)
+    where 
+        supported :: Prim a -> Bool
+        supported PrimWord64 = True
+        supported PrimDouble = True
+        supported _          = False
+        
+        genPrim :: Prim a -> (PureMT -> (a, PureMT))
+        genPrim PrimWord64 = randomWord64
+        genPrim PrimDouble = randomDouble
+        genPrim p = error ("getRandomPrimFromMTRef: genPrim called for unsupported prim " ++ show p)
+        
+        getThing thing = atomicModifyReference ref $ \(!oldMT) -> case thing oldMT of (!w, !newMT) -> (newMT, w)
+            
 
--- |Similarly, @getRandomWordFromMTState x@ can be used in any \"state\"
+-- |Similarly, @getRandomPrimFromMTState x@ can be used in any \"state\"
 -- monad in the mtl sense whose state is a 'PureMT' generator.
 -- Additionally, the standard mtl state monads have 'MonadRandom' instances
 -- which do precisely that, allowing an easy conversion of 'RVar's and
@@ -64,71 +53,65 @@
 -- @runState . sample :: Distribution d t => d t -> PureMT -> (t, PureMT)@.
 -- 'PureMT' in the type there can be replaced by 'StdGen' or anything else 
 -- satisfying @MonadRandom (State s) => s@).
-getRandomWordFromMTState :: MonadState PureMT m => m Word64
-getRandomWordFromMTState = do
-    mt <- get
-    let (ws, newMt) = randomWord64 mt
-    put newMt
-    return ws
-
-getRandomByteFromMTState :: MonadState PureMT m => m Word8
-getRandomByteFromMTState = do
-    mt <- get
-    let (ws, newMt) = randomInt mt
-    put newMt
-    return (fromIntegral ws)
-
-
-getRandomDoubleFromMTState :: MonadState PureMT m => m Double
-getRandomDoubleFromMTState = liftM wordToDouble getRandomWordFromMTState
--- getRandomDoubleFromMTState = do
---     mt <- get
---     let (x, newMt) = randomDouble mt
---     put newMt
---     return x
+getRandomPrimFromMTState :: MonadState PureMT m => Prim a -> m a
+getRandomPrimFromMTState prim
+    | supported prim = getThing (genPrim prim)
+    | otherwise = runPromptM getRandomPrimFromMTState (decomposePrimWhere supported prim)
+    where
+        supported :: Prim a -> Bool
+        supported PrimWord64 = True
+        supported PrimDouble = True
+        supported _          = False
+        
+        genPrim :: Prim a -> (PureMT -> (a, PureMT))
+        genPrim PrimWord64 = randomWord64
+        genPrim PrimDouble = randomDouble
+        genPrim p = error ("getRandomPrimFromMTRef: genPrim called for unsupported prim " ++ show p)
+        
+        getThing thing = do
+            !mt <- get
+            let (!ws, !newMt) = thing mt
+            put newMt
+            return ws
 
 instance MonadRandom (State PureMT) where
-    getRandomByte   = getRandomByteFromMTState
-    getRandomWord   = getRandomWordFromMTState
-    getRandomDouble = getRandomDoubleFromMTState
+    supportedPrims _ _ = True
+    getSupportedRandomPrim = getRandomPrimFromMTState
 
 instance MonadRandom (S.State PureMT) where
-    getRandomByte   = getRandomByteFromMTState
-    getRandomWord   = getRandomWordFromMTState
-    getRandomDouble = getRandomDoubleFromMTState
+    supportedPrims _ _ = True
+    getSupportedRandomPrim = getRandomPrimFromMTState
 
 instance (Monad m1, ModifyRef (Ref m2 PureMT) m1 PureMT) => RandomSource m1 (Ref m2 PureMT) where
-    getRandomByteFrom   = getRandomByteFromMTRef
-    getRandomWordFrom   = getRandomWordFromMTRef
-    getRandomDoubleFrom = getRandomDoubleFromMTRef
-
+    supportedPrimsFrom _ _ = True
+    getSupportedRandomPrimFrom = getRandomPrimFromMTRef
+    
 instance Monad m => MonadRandom (StateT PureMT m) where
-    getRandomByte   = getRandomByteFromMTState
-    getRandomWord   = getRandomWordFromMTState
-    getRandomDouble = getRandomDoubleFromMTState
+    supportedPrims _ _ = True
+    getSupportedRandomPrim = getRandomPrimFromMTState
 
 instance Monad m => MonadRandom (S.StateT PureMT m) where
-    getRandomByte   = getRandomByteFromMTState
-    getRandomWord   = getRandomWordFromMTState
-    getRandomDouble = getRandomDoubleFromMTState
+    supportedPrims _ _ = True
+    getSupportedRandomPrim = getRandomPrimFromMTState
 
 instance (Monad m, ModifyRef (IORef PureMT) m PureMT) => RandomSource m (IORef PureMT) where
     {-# SPECIALIZE instance RandomSource IO (IORef PureMT)#-}
-    getRandomByteFrom   = getRandomByteFromMTRef
-    getRandomWordFrom   = getRandomWordFromMTRef
-    getRandomDoubleFrom = getRandomDoubleFromMTRef
-
+    supportedPrimsFrom _ _ = True
+    getSupportedRandomPrimFrom = getRandomPrimFromMTRef
+    
 instance (Monad m, ModifyRef (STRef s PureMT) m PureMT) => RandomSource m (STRef s PureMT) where
     {-# SPECIALIZE instance RandomSource (ST s) (STRef s PureMT) #-}
     {-# SPECIALIZE instance RandomSource (S.ST s) (STRef s PureMT) #-}
-    getRandomByteFrom   = getRandomByteFromMTRef
-    getRandomWordFrom   = getRandomWordFromMTRef
-    getRandomDoubleFrom = getRandomDoubleFromMTRef
-
-instance (Monad m, ModifyRef (TVar PureMT) m PureMT) => RandomSource m (TVar PureMT) where
-    {-# SPECIALIZE instance RandomSource IO  (TVar PureMT) #-}
-    {-# SPECIALIZE instance RandomSource STM (TVar PureMT) #-}
-    getRandomByteFrom   = getRandomByteFromMTRef
-    getRandomWordFrom   = getRandomWordFromMTRef
-    getRandomDoubleFrom = getRandomDoubleFromMTRef
+    supportedPrimsFrom _ _ = True
+    getSupportedRandomPrimFrom = getRandomPrimFromMTRef
 
+-- Note that this instance is probably a Bad Idea.  STM allows random variables
+-- to interact in spooky quantum-esque ways - One transaction can 'retry' until
+-- it gets a \"random\" answer it likes, which causes it to selectively consume 
+-- entropy, biasing the supply from which other random variables will draw.
+-- instance (Monad m, ModifyRef (TVar PureMT) m PureMT) => RandomSource m (TVar PureMT) where
+--     {-# SPECIALIZE instance RandomSource IO  (TVar PureMT) #-}
+--     {-# SPECIALIZE instance RandomSource STM (TVar PureMT) #-}
+--     supportedPrimsFrom _ _ = True
+--     getSupportedRandomPrimFrom = getRandomPrimFromMTRef
+    
diff --git a/src/Data/Random/Source/Std.hs b/src/Data/Random/Source/Std.hs
--- a/src/Data/Random/Source/Std.hs
+++ b/src/Data/Random/Source/Std.hs
@@ -8,6 +8,7 @@
 module Data.Random.Source.Std where
 
 import Data.Random.Source
+import Data.Tagged
 
 -- |A token representing the \"standard\" entropy source in a 'MonadRandom'
 -- monad.  Its sole purpose is to make the following true (when the types check):
@@ -16,6 +17,11 @@
 data StdRandom = StdRandom
 
 instance MonadRandom m => RandomSource m StdRandom where
-    getRandomByteFrom   StdRandom = getRandomByte
-    getRandomWordFrom   StdRandom = getRandomWord
-    getRandomDoubleFrom StdRandom = getRandomDouble
+    {-SPECIALIZE instance MonadRandom m => RandomSource m StdRandom -}
+    supportedPrimsFrom w = supportedPrims (mkWit w)
+        where
+            mkWit :: Tagged a b -> a
+            mkWit = error "supportedPrims tried to evaluate its phantom parameter"
+    getSupportedRandomPrimFrom   StdRandom = getSupportedRandomPrim
+    
+    getRandomPrimFrom StdRandom = getRandomPrim
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
@@ -1,9 +1,8 @@
-{-
- -      ``Data/Random/Source/StdGen''
- -}
 {-# LANGUAGE
-    MultiParamTypeClasses, FlexibleInstances, UndecidableInstances
+    MultiParamTypeClasses, FlexibleInstances, UndecidableInstances, GADTs,
+    BangPatterns, RankNTypes
   #-}
+{-# OPTIONS_GHC -fno-warn-orphans #-}
 
 -- |This module provides functions useful for implementing new 'MonadRandom'
 -- and 'RandomSource' instances for state-abstractions containing 'StdGen'
@@ -13,131 +12,165 @@
 module Data.Random.Source.StdGen where
 
 import Data.Random.Internal.Words
+import Data.Random.Internal.Primitives
 import Data.Random.Source
 import System.Random
+import Control.Monad.Prompt
 import Control.Monad.State
 import qualified Control.Monad.ST.Strict as S
 import qualified Control.Monad.State.Strict as S
 import Data.StateRef
 import Data.Word
 
+
 instance (Monad m1, ModifyRef (Ref m2 StdGen) m1 StdGen) => RandomSource m1 (Ref m2 StdGen) where
-    getRandomByteFrom   = getRandomByteFromRandomGenRef
-    getRandomWordFrom   = getRandomWordFromRandomGenRef
-    getRandomDoubleFrom = getRandomDoubleFromRandomGenRef
+    supportedPrimsFrom _ _ = True
+    getSupportedRandomPrimFrom = getRandomPrimFromRandomGenRef
 
 instance (Monad m, ModifyRef (IORef   StdGen) m StdGen) => RandomSource m (IORef   StdGen) where
     {-# SPECIALIZE instance RandomSource IO (IORef StdGen) #-}
-    getRandomByteFrom   = getRandomByteFromRandomGenRef
-    getRandomWordFrom   = getRandomWordFromRandomGenRef
-    getRandomDoubleFrom = getRandomDoubleFromRandomGenRef
-instance (Monad m, ModifyRef (TVar    StdGen) m StdGen) => RandomSource m (TVar    StdGen) where
-    {-# SPECIALIZE instance RandomSource IO  (TVar StdGen) #-}
-    {-# SPECIALIZE instance RandomSource STM (TVar StdGen) #-}
-    getRandomByteFrom   = getRandomByteFromRandomGenRef
-    getRandomWordFrom   = getRandomWordFromRandomGenRef
-    getRandomDoubleFrom = getRandomDoubleFromRandomGenRef
+    supportedPrimsFrom _ _ = True
+    getSupportedRandomPrimFrom = getRandomPrimFromRandomGenRef
+
+-- Note that this instance is probably a Bad Idea.  STM allows random variables
+-- to interact in spooky quantum-esque ways - One transaction can 'retry' until
+-- it gets a \"random\" answer it likes, which causes it to selectively consume 
+-- entropy, biasing the supply from which other random variables will draw.
+-- instance (Monad m, ModifyRef (TVar    StdGen) m StdGen) => RandomSource m (TVar    StdGen) where
+--     {-# SPECIALIZE instance RandomSource IO  (TVar StdGen) #-}
+--     {-# SPECIALIZE instance RandomSource STM (TVar StdGen) #-}
+--     supportedPrimsFrom _ _ = True
+--     getSupportedRandomPrimFrom = getRandomPrimFromRandomGenRef
+
 instance (Monad m, ModifyRef (STRef s StdGen) m StdGen) => RandomSource m (STRef s StdGen) where
     {-# SPECIALIZE instance RandomSource (ST s) (STRef s StdGen) #-}
     {-# SPECIALIZE instance RandomSource (S.ST s) (STRef s StdGen) #-}
-    getRandomByteFrom   = getRandomByteFromRandomGenRef
-    getRandomWordFrom   = getRandomWordFromRandomGenRef
-    getRandomDoubleFrom = getRandomDoubleFromRandomGenRef
-
-getRandomByteFromStdGenIO :: IO Word8
-getRandomByteFromStdGenIO = do
-    int <- randomRIO (0, 255) :: IO Int
-    return (fromIntegral int)
+    supportedPrimsFrom _ _ = True
+    getSupportedRandomPrimFrom = getRandomPrimFromRandomGenRef
 
-getRandomWordFromStdGenIO :: IO Word64
-getRandomWordFromStdGenIO = do
-    int <- randomRIO (0, 0xffffffffffffffff)
-    return (fromInteger int)
+getRandomPrimFromStdGenIO :: Prim a -> IO a
+getRandomPrimFromStdGenIO prim
+    | supported prim = genPrim prim
+    | otherwise = runPromptM getRandomPrimFromStdGenIO (decomposePrimWhere supported prim)
+    where 
+        {-# INLINE supported #-}
+        supported :: Prim a -> Bool
+        supported PrimWord8             = True
+        supported PrimWord16            = True
+        supported PrimWord32            = True
+        supported PrimWord64            = True
+        supported PrimDouble            = True
+        supported (PrimNByteInteger _)  = True
+        supported _                     = False
+        
+        -- based on reading the source of the "random" library's implementation, I do
+        -- not believe that the randomRIO (0,1) implementation for Double is capable of producing
+        -- the value 0.  Therefore, I'm not using it.  If this is an incorrect reading on
+        -- my part, or if this changes, then feel free to change the implementation.
+        -- Same goes for the other getRandomDouble... functions here.
 
--- based on reading the source of the "random" library's implementation, I do
--- not believe that the randomRIO (0,1) implementation for Double is capable of producing
--- the value 0.  Therefore, I'm not using it.  If this is an incorrect reading on
--- my part, or if this changes, then feel free to use the commented version.
--- Same goes for the other getRandomDouble... functions here.
-getRandomDoubleFromStdGenIO :: IO Double
-getRandomDoubleFromStdGenIO = liftM wordToDouble getRandomWordFromStdGenIO
--- getRandomDoubleFromStdGenIO = randomRIO (0, 1)
+        {-# INLINE genPrim #-}
+        genPrim :: Prim a -> IO a
+        genPrim PrimWord8            = fmap fromIntegral                  (randomRIO (0, 0xff) :: IO Int)
+        genPrim PrimWord16           = fmap fromIntegral                  (randomRIO (0, 0xffff) :: IO Int)
+        genPrim PrimWord32           = fmap fromInteger                   (randomRIO (0, 0xffffffff))
+        genPrim PrimWord64           = fmap fromInteger                   (randomRIO (0, 0xffffffffffffffff))
+        genPrim PrimDouble           = fmap (wordToDouble . fromInteger)  (randomRIO (0, 0xffffffffffffffff))
+        genPrim (PrimNByteInteger n) = randomRIO (0, iterate (*256) 1 !! n)
+        genPrim p = error ("getRandomPrimFromStdGenIO: genPrim called for unsupported prim " ++ show p)
 
 -- |Given a mutable reference to a 'RandomGen' generator, we can make a
 -- 'RandomSource' usable in any monad in which the reference can be modified.
---
--- For example, if @x :: TVar StdGen@, @getRandomByteFromRandomGenRef x@ can be
--- used as a 'RandomSource' in 'IO', 'STM', or any monad which is an instance
--- of 'MonadIO'.  It's generally probably better to use
--- 'getRandomWordFromRandomGenRef' though, as this one is likely to throw
--- away a lot of perfectly good entropy.  Better still is to use these 3 functions
--- together to create a 'RandomSource' instance for the reference you're using,
--- if one does not already exist.
-getRandomByteFromRandomGenRef :: (Monad m, ModifyRef sr m g, RandomGen g) =>
-                                  sr -> m Word8
-getRandomByteFromRandomGenRef g = atomicModifyReference g (swap . randomR (0,255))
+getRandomPrimFromRandomGenRef :: (Monad m, ModifyRef sr m g, RandomGen g) =>
+                                  sr -> Prim a -> m a
+getRandomPrimFromRandomGenRef ref prim
+    | supported prim = genPrim prim getThing
+    | otherwise = runPromptM (getRandomPrimFromRandomGenRef ref) (decomposePrimWhere supported prim)
     where 
-        swap :: (Int, a) -> (a, Word8)
-        swap (a,b) = (b,fromIntegral a)
+        {-# INLINE supported #-}
+        supported :: Prim a -> Bool
+        supported PrimWord8             = True
+        supported PrimWord16            = True
+        supported PrimWord32            = True
+        supported PrimWord64            = True
+        supported PrimDouble            = True
+        supported (PrimNByteInteger _)  = True
+        supported _                     = False
+        
+        {-# INLINE genPrim #-}
+        genPrim :: (RandomGen g) => Prim a -> (forall b. (g -> (b, g)) -> (b -> a) -> c) -> c
+        genPrim PrimWord8            f = f (randomR (0, 0xff))                (fromIntegral :: Int -> Word8)
+        genPrim PrimWord16           f = f (randomR (0, 0xffff))              (fromIntegral :: Int -> Word16)
+        genPrim PrimWord32           f = f (randomR (0, 0xffffffff))          (fromInteger)
+        genPrim PrimWord64           f = f (randomR (0, 0xffffffffffffffff))  (fromInteger)
+        genPrim PrimDouble           f = f (randomR (0, 0x000fffffffffffff))  (flip encodeFloat (-52))
+        genPrim (PrimNByteInteger n) f = f (randomR (0, iterate (*256) 1 !! n)) (id :: Integer -> Integer)
+        genPrim p _ = error ("getRandomPrimFromRandomGenRef: genPrim called for unsupported prim " ++ show p)
+        
+        {-# INLINE getThing #-}
+        getThing thing f = atomicModifyReference ref $ \(!oldMT) -> case thing oldMT of (!w, !newMT) -> (newMT, f w)
 
-getRandomWordFromRandomGenRef :: (Monad m, ModifyRef sr m g, RandomGen g) =>
-                                  sr -> m Word64
-getRandomWordFromRandomGenRef g = atomicModifyReference g (swap . randomR (0,0xffffffffffffffff))
-    where swap (a,b) = (b,fromInteger a)
 
-getRandomDoubleFromRandomGenRef :: (Monad m, ModifyRef sr m g, RandomGen g) =>
-                                  sr -> m Double
-getRandomDoubleFromRandomGenRef g = liftM wordToDouble (getRandomWordFromRandomGenRef g)
--- getRandomDoubleFromRandomGenRef g = atomicModifyRef g (swap . randomR (0,1))
---     where swap (a,b) = (b,a)
-
 -- |Similarly, @getRandomWordFromRandomGenState x@ can be used in any \"state\"
 -- monad in the mtl sense whose state is a 'RandomGen' generator.
 -- Additionally, the standard mtl state monads have 'MonadRandom' instances
 -- which do precisely that, allowing an easy conversion of 'RVar's and
 -- other 'Distribution' instances to \"pure\" random variables.
-getRandomByteFromRandomGenState :: (RandomGen g, MonadState g m) => m Word8
-getRandomByteFromRandomGenState = do
-    g <- get
-    case randomR (0, 255 :: Int) g of
-        (i,g) -> do
-            put g
-            return (fromIntegral i)
-
-getRandomWordFromRandomGenState :: (RandomGen g, MonadState g m) => m Word64
-getRandomWordFromRandomGenState = do
-    g <- get
-    case randomR (0, 0xffffffffffffffff) g of
-        (i,g) -> do
-            put g
-            return (fromInteger i)
-
-getRandomDoubleFromRandomGenState :: (RandomGen g, MonadState g m) => m Double
-getRandomDoubleFromRandomGenState = liftM wordToDouble getRandomWordFromRandomGenState
--- getRandomDoubleFromRandomGenState = do
---     g <- get
---     case randomR (0, 1) g of
---         (x,g) -> do
---             put g
---             return x
-
+{-# SPECIALIZE getRandomPrimFromRandomGenState :: Prim a -> State StdGen a #-}
+{-# SPECIALIZE getRandomPrimFromRandomGenState :: Monad m => Prim a -> StateT StdGen m a #-}
+getRandomPrimFromRandomGenState :: (RandomGen g, MonadState g m) => Prim a -> m a
+getRandomPrimFromRandomGenState prim
+    | supported prim = genSupported prim
+    | otherwise = runPromptM genSupported (decomposePrimWhere supported prim)
+    where 
+        {-# INLINE genSupported #-}
+        genSupported prim = genPrim prim getThing
+        
+        {-# INLINE supported #-}
+        supported :: Prim a -> Bool
+        supported PrimWord8             = True
+        supported PrimWord16            = True
+        supported PrimWord32            = True
+        supported PrimWord64            = True
+        supported PrimDouble            = True
+        supported (PrimNByteInteger _)  = True
+        supported _                     = False
+        
+        {-# INLINE genPrim #-}
+        genPrim :: (RandomGen g) => Prim a -> (forall b. (g -> (b, g)) -> (b -> a) -> c) -> c
+        genPrim PrimWord8            f = f (randomR (0, 0xff))                (fromIntegral :: Int -> Word8)
+        genPrim PrimWord16           f = f (randomR (0, 0xffff))              (fromIntegral :: Int -> Word16)
+        genPrim PrimWord32           f = f (randomR (0, 0xffffffff))          (fromInteger)
+        genPrim PrimWord64           f = f (randomR (0, 0xffffffffffffffff))  (fromInteger)
+        genPrim PrimDouble           f = f (randomR (0, 0x000fffffffffffff))  (flip encodeFloat (-52))
+          {- not using the Random Double instance for 2 reasons.  1st, it only generates 32 bits of entropy, when 
+             a [0,1) Double has room for 52.  Second, it appears there's a bug where it can actually generate a 
+             negative number in the case where randomIvalInteger returns minBound::Int32. -}
+--        genPrim PrimDouble f = f (randomR (0, 1.0))  (id)
+        genPrim (PrimNByteInteger n) f = f (randomR (0, iterate (*256) 1 !! n)) id
+        genPrim p _ = error ("getRandomPrimFromRandomGenState: genPrim called for unsupported prim " ++ show p)
+        
+        {-# INLINE getThing #-}
+        getThing thing f = do
+            !oldGen <- get
+            case thing oldGen of
+                (!i,!newGen) -> do
+                    put newGen
+                    return (f $! i)
 
 instance MonadRandom (State StdGen) where
-    getRandomByte   = getRandomByteFromRandomGenState
-    getRandomWord   = getRandomWordFromRandomGenState
-    getRandomDouble = getRandomDoubleFromRandomGenState
+    supportedPrims _ _ = True
+    getSupportedRandomPrim = getRandomPrimFromRandomGenState
 
 instance Monad m => MonadRandom (StateT StdGen m) where
-    getRandomByte   = getRandomByteFromRandomGenState
-    getRandomWord   = getRandomWordFromRandomGenState
-    getRandomDouble = getRandomDoubleFromRandomGenState
+    supportedPrims _ _ = True
+    getSupportedRandomPrim = getRandomPrimFromRandomGenState
 
 instance MonadRandom (S.State StdGen) where
-    getRandomByte   = getRandomByteFromRandomGenState
-    getRandomWord   = getRandomWordFromRandomGenState
-    getRandomDouble = getRandomDoubleFromRandomGenState
+    supportedPrims _ _ = True
+    getSupportedRandomPrim = getRandomPrimFromRandomGenState
 
 instance Monad m => MonadRandom (S.StateT StdGen m) where
-    getRandomByte   = getRandomByteFromRandomGenState
-    getRandomWord   = getRandomWordFromRandomGenState
-    getRandomDouble = getRandomDoubleFromRandomGenState
+    supportedPrims _ _ = True
+    getSupportedRandomPrim = getRandomPrimFromRandomGenState
+