diff --git a/random-fu.cabal b/random-fu.cabal
--- a/random-fu.cabal
+++ b/random-fu.cabal
@@ -1,5 +1,5 @@
 name:                   random-fu
-version:                0.0.3
+version:                0.0.3.2
 stability:              experimental
 
 cabal-version:          >= 1.2
@@ -19,6 +19,8 @@
                         as well as reasonably fast.
 
 Flag base4
+Flag base4_2
+    Description:        base-4.2 has an incompatible change in Data.Fixed (HasResolution)
 
 Library
   hs-source-dirs:       src
@@ -50,8 +52,13 @@
                         Data.Random.Source.Std
                         Data.Random.Source.StdGen
   if flag(base4)
-    build-depends:      base >= 4 && <5,
-                        syb
+    build-depends:      syb
+    
+    if flag(base4_2)
+      build-depends:    base >= 4 && <4.2
+    else
+      cpp-options:      -Dbase_4_2
+      build-depends:    base >= 4.2 && <5
   else
     build-depends:      base >= 3 && < 4
     
diff --git a/src/Data/Random.hs b/src/Data/Random.hs
--- a/src/Data/Random.hs
+++ b/src/Data/Random.hs
@@ -7,14 +7,14 @@
 
 -- |Random numbers and stuff...
 -- 
--- Data.Random.Source exports the typeclasses for entropy sources, and
+-- "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.
 -- 
--- Data.Random.Distribution exports the typeclasses for sampling distributions,
+-- "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
+-- "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.
 
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
@@ -16,7 +16,6 @@
 class Distribution d t where
     -- |Return a random variable with this distribution.
     rvar :: d t -> RVar t
-    rvar = rvarT
 
 class Distribution d t => CDF d t where
     -- |Return the cumulative distribution function of this distribution.
@@ -34,7 +33,9 @@
     -- 
     -- Thus, '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.
+    -- 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.
     cdf :: d t -> t -> Double
 
 -- |Return a random variable with the given distribution, pre-lifted to an arbitrary 'RVarT'.
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
@@ -36,9 +36,17 @@
 
 import Data.Number.Erf
 
+-- |A random variable that produces a pair of independent
+-- normally-distributed values.
 normalPair :: (Floating a, Distribution StdUniform a) => RVar (a,a)
 normalPair = boxMullerNormalPair
 
+-- |A random variable that produces a pair of independent
+-- normally-distributed values, computed using the Box-Muller method.
+-- This algorithm is slightly slower than Knuth's method but using a 
+-- constant amount of entropy (Knuth's method is a rejection method).
+-- It is also slightly more general (Knuth's method require an 'Ord'
+-- instance).
 {-# INLINE boxMullerNormalPair #-}
 boxMullerNormalPair :: (Floating a, Distribution StdUniform a) => RVar (a,a)
 boxMullerNormalPair = do
@@ -51,6 +59,10 @@
         y = r * sin theta
     return (x,y)
 
+-- |A random variable that produces a pair of independent
+-- normally-distributed values, computed using Knuth's polar method.
+-- Slightly faster than 'boxMullerNormalPair' when it accepts on the 
+-- first try, but does not always do so.
 {-# INLINE knuthPolarNormalPair #-}
 knuthPolarNormalPair :: (Floating a, Ord a, Distribution Uniform a) => RVar (a,a)
 knuthPolarNormalPair = do
@@ -65,8 +77,7 @@
             else let scale = sqrt (-2 * log s / s) 
                   in (v1 * scale, v2 * scale)
 
--- |Draw from the tail of a normal distribution (the region beyond the provided value), 
--- returning a negative value if the Bool parameter is True.
+-- |Draw from the tail of a normal distribution (the region beyond the provided value)
 {-# INLINE normalTail #-}
 normalTail :: (Distribution StdUniform a, Floating a, Ord a) =>
               a -> RVar a
@@ -82,13 +93,13 @@
                 else return (r - x)
 
 -- |Construct a 'Ziggurat' for sampling a normal distribution, given
--- logBase 2 c, and the 'zGetIU' implementation.
+-- @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
 normalZ p = mkZigguratRec True normalF normalFInv normalFInt normalFVol (2^p)
 
--- | Ziggurat target function
+-- | Ziggurat target function (upper half of a non-normalized gaussian PDF)
 normalF :: (Floating a, Ord a) => a -> a
 normalF x
     | x <= 0    = 1
@@ -105,6 +116,22 @@
 normalFVol :: Floating a => a
 normalFVol = sqrt (0.5 * pi)
 
+-- |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).  
+-- 
+-- Because it computes a 'Ziggurat', it is very expensive to use for
+-- just one evaluation, or even for multiple evaluations if not used and
+-- reused monomorphically (to enable the ziggurat table to be let-floated
+-- out).  If you don't know whether your use case fits this description
+-- then you're probably better off using a different algorithm, such as
+-- 'boxMullerNormalPair' or 'knuthPolarNormalPair'.  And of course if
+-- you don't need the full generality of this definition then you're much
+-- better off using 'doubleStdNormal' or 'floatStdNormal'.
+--
+-- 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)
     where 
@@ -115,6 +142,8 @@
             u <- uniform (-1) 1
             return (fromIntegral i .&. (2^p-1), u)
 
+-- |A random variable sampling from the standard normal distribution
+-- over the 'Double' type.
 doubleStdNormal :: RVar Double
 doubleStdNormal = runZiggurat doubleStdNormalZ
 
@@ -137,6 +166,8 @@
             let (u,i) = wordToDoubleWithExcess w
             return (fromIntegral i .&. (doubleStdNormalC-1), u+u-1)
 
+-- |A random variable sampling from the standard normal distribution
+-- over the 'Float' type.
 floatStdNormal :: RVar Float
 floatStdNormal = runZiggurat floatStdNormalZ
 
@@ -165,8 +196,11 @@
 normalCdf :: (Real a) => a -> a -> a -> Double
 normalCdf m s x = normcdf ((realToFrac x - realToFrac m) / realToFrac s)
 
+-- |A specification of a normal distribution over the type 'a'.
 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 a a -- mean, sd
 
 instance Distribution Normal Double where
@@ -189,8 +223,10 @@
 
 {-# SPECIALIZE stdNormal :: RVar Double #-}
 {-# SPECIALIZE stdNormal :: RVar Float #-}
+-- |'stdNormal' is a normal variable with distribution 'StdNormal'.
 stdNormal :: Distribution Normal a => RVar a
 stdNormal = rvar StdNormal
 
+-- |@normal m s@ is a random variable with distribution @'Normal' m s@.
 normal :: Distribution Normal a => a -> a -> RVar a
 normal m s = rvar (Normal m s)
diff --git a/src/Data/Random/Distribution/Rayleigh.hs b/src/Data/Random/Distribution/Rayleigh.hs
--- a/src/Data/Random/Distribution/Rayleigh.hs
+++ b/src/Data/Random/Distribution/Rayleigh.hs
@@ -17,14 +17,11 @@
 floatingRayleigh s = do
     u <- stdUniformPos
     return (s * sqrt (-2 * log u))
-    
-    where 
-        stdUniformPos = do
-            u <- stdUniform
-            if u == 0
-                then stdUniformPos
-                else return u
 
+-- |The rayleigh distribution with a specified mode (\"sigma\") parameter.
+-- Its mean will be @sigma*sqrt(pi/2)@ and its variance will be @sigma^2*(4-pi)/2@
+-- 
+-- (therefore if you want one with a particular mean @m@, @sigma@ should be @m*sqrt(2/pi)@)
 newtype Rayleigh a = Rayleigh a
 
 rayleigh :: Distribution Rayleigh a => a -> RVar a
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
@@ -13,12 +13,23 @@
 import Data.Random.Distribution
 import Data.Random.Distribution.Uniform
 
-data Triangular a = Triangular
-    { triLower  :: a
-    , triMid    :: a
-    , triUpper  :: a
-    } deriving (Eq, Show)
+-- |A description of a triangular distribution - a distribution whose PDF
+-- is a triangle ramping up from a lower bound to a specified midpoint
+-- and back down to the upper bound.  This is a very simple distribution
+-- that does not generally occur naturally but is used sometimes as an
+-- estimate of a true distribution when only the range of the values and
+-- an approximate mode of the true distribution are known.
+data Triangular a = Triangular {
+    -- |The lower bound of the triangle in the PDF (the smallest number the distribution can generate)
+    triLower  :: a,
+    -- |The midpoint of the triangle (also the mode of the distribution)
+    triMid    :: a,
+    -- |The upper bound of the triangle (and the largest number the distribution can generate)
+    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
@@ -33,6 +44,7 @@
 --        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
     | x < a
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
@@ -14,6 +14,8 @@
 	
     , StdUniform(..)
     , stdUniform
+    , stdUniformNonneg
+    , stdUniformPos
     
     , integralUniform
     , realFloatUniform
@@ -47,6 +49,7 @@
 
 import Control.Monad.Loops
 
+-- |Compute a random 'Integral' value between the 2 values provided (inclusive).
 integralUniform :: (Integral a) => a -> a -> RVar a
 integralUniform a b
     | a > b     = compute b a
@@ -71,26 +74,33 @@
     Just x -> x
 powersOf256 = 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
 boundedStdUniform = uniform minBound maxBound
 
 boundedStdUniformCDF :: (CDF Uniform a, Bounded a) => a -> Double
 boundedStdUniformCDF = cdf (Uniform minBound maxBound)
 
+-- |Compute a random value for a 'Bounded' 'Enum' type, between 'minBound' and
+-- 'maxBound' (inclusive)
 boundedEnumStdUniform :: (Enum a, Bounded a) => RVar a
 boundedEnumStdUniform = enumUniform minBound maxBound
 
 boundedEnumStdUniformCDF :: (Enum a, Bounded a, Ord a) => a -> Double
 boundedEnumStdUniformCDF = enumUniformCDF minBound maxBound
 
+-- |Compute a uniform random 'Float' value in the range [0,1)
 floatStdUniform :: RVar Float
 floatStdUniform = do
     x <- getRandomWord
     return (wordToFloat x)
 
+-- |Compute a uniform random 'Double' value in the range [0,1)
 doubleStdUniform :: RVar Double
 doubleStdUniform = getRandomDouble
 
+-- |Compute a uniform random value in the range [0,1) for any 'RealFloat' type 
 realFloatStdUniform :: RealFloat a => RVar a
 realFloatStdUniform = do
     let (b, e) = decodeFloat one
@@ -102,6 +112,8 @@
     
     where one = 1
 
+-- |Compute a uniform random 'Fixed' value in the range [0,1), with any
+-- desired precision.
 fixedStdUniform :: HasResolution r => RVar (Fixed r)
 fixedStdUniform = x
     where
@@ -110,42 +122,52 @@
             u <- uniform 0 (res)
             return (mkFixed u)
 
+-- |The CDF of the random variable 'realFloatStdUniform'.
 realStdUniformCDF :: Real a => a -> Double
 realStdUniformCDF x
     | x <= 0    = 0
     | x >= 1    = 1
     | otherwise = realToFrac x
 
+-- |@floatUniform a b@ computes a uniform random 'Float' value in the range [a,b)
 floatUniform :: Float -> Float -> RVar Float
 floatUniform 0 1 = floatStdUniform
 floatUniform a b = do
     x <- floatStdUniform
     return (a + x * (b - a))
 
+-- |@doubleUniform a b@ computes a uniform random 'Double' value in the range [a,b)
 doubleUniform :: Double -> Double -> RVar Double
 doubleUniform 0 1 = doubleStdUniform
 doubleUniform a b = do
     x <- doubleStdUniform
     return (a + x * (b - a))
 
+-- |@realFloatUniform a b@ computes a uniform random value in the range [a,b) for
+-- any 'RealFloat' type
 realFloatUniform :: RealFloat a => a -> a -> RVar a
 realFloatUniform 0 1 = realFloatStdUniform
 realFloatUniform a b = do
     x <- realFloatStdUniform
     return (a + x * (b - a))
 
+-- |@fixedUniform a b@ computes a uniform random 'Fixed' value in the range 
+-- [a,b), with any desired precision.
 fixedUniform :: HasResolution r => Fixed r -> Fixed r -> RVar (Fixed r)
 fixedUniform a b = do
     u <- integralUniform (unMkFixed a) (unMkFixed b)
     return (mkFixed u)
 
-realUniformCDF :: Real a => a -> a -> a -> Double
+-- |@realUniformCDF a b@ is the CDF of the random variable @realFloatUniform a b@.
+realUniformCDF :: RealFrac a => a -> a -> a -> Double
 realUniformCDF a b x
     | b < a     = realUniformCDF b a x
     | x <= a    = 0
     | x >= b    = 1
-    | otherwise = realToFrac (x-a) / realToFrac (b-a)
+    | otherwise = realToFrac ((x-a) / (b-a))
 
+-- |@realFloatUniform a b@ computes a uniform random value in the range [a,b) for
+-- any 'Enum' type
 enumUniform :: Enum a => a -> a -> RVar a
 enumUniform a b = do
     x <- integralUniform (fromEnum a) (fromEnum b)
@@ -160,6 +182,9 @@
     
     where e2f = fromIntegral . fromEnum
 
+-- @uniform a b@ computes a uniformly distributed random value in the range
+-- [a,b] for 'Integral' or 'Enum' types and in the range [a,b) for 'Fractional'
+-- types.  Requires a @Distribution Uniform@ instance for the type.
 uniform :: Distribution Uniform a => a -> a -> RVar a
 uniform a b = rvar (Uniform a b)
 
@@ -170,14 +195,33 @@
 stdUniform :: (Distribution StdUniform a) => RVar a
 stdUniform = rvar StdUniform
 
-stdUniformPos :: (Distribution StdUniform a, Ord a, Num a) => RVar a
+-- |Like 'stdUniform', but uses 'abs' to return only positive or zero values.
+stdUniformNonneg :: (Distribution StdUniform a, Num a) => RVar a
+stdUniformNonneg = abs `fmap` stdUniform
+
+-- |Like 'stdUniform' but only returns positive values.
+stdUniformPos :: (Distribution StdUniform a, Num a) => RVar a
 stdUniformPos = do
-    x <- stdUniform
-    if x > 0
+    x <- stdUniformNonneg
+    if x /= 0
         then return x
         else stdUniformPos
 
-data Uniform t = Uniform !t !t
+-- |A definition of a uniform distribution over the type @t@.  See also 'uniform'.
+data Uniform t = 
+    -- |A uniform distribution defined by a lower and upper range bound.
+    -- For 'Integral' and 'Enum' types, the range is inclusive.  For 'Fractional'
+    -- types the range includes the lower bound but not the upper.
+    Uniform !t !t
+
+-- |A name for the \"standard\" uniform distribution over the type @t@,
+-- if one exists.  See also 'stdUniform'.
+--
+-- For 'Integral' and 'Enum' types that are also 'Bounded', this is
+-- the uniform distribution over the full range of the type.
+-- For un-'Bounded' 'Integral' types this is not defined.
+-- For 'Fractional' types this is a random variable in the range [0,1)
+-- (that is, 0 to 1 including 0 but not including 1).
 data StdUniform t = StdUniform
 
 $( replicateInstances ''Int integralTypes [d|
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,24 +1,22 @@
-{-
- -      ``Data/Random/Distribution/Ziggurat''
- -      A generic "ziggurat algorithm" implementation.  Fairly rough right
- -      now.
- -      
- -      There is a lot of room for improvement in 'findBin0' especially.
- -      It needs a fair amount of cleanup and elimination of redundant
- -      calculation, as well as either a justification for using the simple
- -      'findMinFrom' or a proper root-finding algorithm. 
- -      
- -      It would also be nice to add (preferably via its own library)
- -      support for numerical integration and differentiation, so that
- -      tables can be derived from only a PDF (if the end user is
- -      willing to take the performance hit for the convenience).
- -}
 {-# LANGUAGE
         MultiParamTypeClasses,
         FlexibleInstances, FlexibleContexts,
         RecordWildCards
   #-}
 
+-- |A generic \"ziggurat algorithm\" implementation.  Fairly rough right
+--  now.
+--  
+--  There is a lot of room for improvement in 'findBin0' especially.
+--  It needs a fair amount of cleanup and elimination of redundant
+--  calculation, as well as either a justification for using the simple
+--  'findMinFrom' or a proper root-finding algorithm. 
+--  
+--  It would also be nice to add (preferably by pulling in an 
+--  external package) support for numerical integration and 
+--  differentiation, so that tables can be derived from only a 
+--  PDF (if the end user is willing to take the performance and 
+--  accuracy hit for the convenience).
 module Data.Random.Distribution.Ziggurat
     ( Ziggurat(..)
     , mkZigguratRec
@@ -38,54 +36,129 @@
 
 vec ! i = index vec i
 
-data Ziggurat t = Ziggurat
-    { zTable_xs         :: Vector t
-    , zTable_x_ratios   :: Vector t
-    , zTable_ys         :: Vector t
-    , zGetIU            :: RVar (Int, t)
-    , zTailDist         :: RVar t
-    , zUniform          :: t -> t -> RVar t
-    , zFunc             :: t -> t
-    , zMirror           :: Bool
+-- |A data structure containing all the data that is needed
+-- to implement Marsaglia & Tang's \"ziggurat\" algorithm for
+-- sampling certain kinds of random distributions.
+--
+-- The documentation here is probably not sufficient to tell a user exactly
+-- how to build one of these from scratch, but it is not really intended to
+-- 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 {
+        -- |The X locations of each bin in the distribution.  Bin 0 is the
+        -- 'infinite' one.
+        -- 
+        -- In the case of bin 0, the value given is sort of magical - x[0] is
+        -- defined to be V/f(R).  It's not actually the location of any bin, 
+        -- but a value computed to make the algorithm more concise and slightly 
+        -- 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,
+        -- |The ratio of each bin's Y value to the next bin's Y value
+        zTable_x_ratios   :: Vector t,
+        -- |The Y value (zFunc x) of each bin
+        zTable_ys         :: Vector t,
+        -- |An RVar providing a random tuple consisting of:
+        --
+        --  * a bin index, uniform over [0,c) :: Int (where @c@ is the
+        --    number of bins in the tables)
+        -- 
+        --  * a uniformly distributed fractional value, from -1 to 1 
+        --    if not mirrored, from 0 to 1 otherwise.
+        --
+        -- This is provided as a single 'RVar' because it can be implemented
+        -- more efficiently than naively sampling 2 separate values - a
+        -- 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),
+        
+        -- |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,
+        
+        -- |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,
+        
+        -- |The (one-sided antitone) PDF, not necessarily normalized
+        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
     }
 
+-- |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
     where
         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
             
+            -- 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
                 else if i == 0
-                    then if x < 0 then fmap negate zTailDist else zTailDist
-                    else do
-                        v <- zUniform (zTable_ys ! (i+1)) (zTable_ys ! i)
-                        if v < zFunc (abs x)
-                            then return $! x
-                            else go
+                    then sampleTail x
+                    else sampleGreyArea i x
+        
+        -- 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.
+        sampleGreyArea i x = do
+            v <- zUniform (zTable_ys ! (i+1)) (zTable_ys ! i)
+            if v < zFunc (abs x)
+                then return $! x
+                else go
+        
+        -- 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)
+        sampleTail x
+            | x < 0     = fmap negate zTailDist
+            | otherwise = zTailDist
 
 
--- |Build the tables to implement the "ziggurat algorithm" devised by 
+-- |Build the tables to implement the \"ziggurat algorithm\" devised by 
 -- Marsaglia & Tang, attempting to automatically compute the R and V
 -- values.
 -- 
 -- Arguments:
 -- 
 --  * flag indicating whether to mirror the distribution
---  * the (one-sided antitone) CDF
---  * the inverse of the CDF
+-- 
+--  * the (one-sided antitone) PDF, not necessarily normalized
+-- 
+--  * the inverse of the PDF
+-- 
 --  * the number of bins
+-- 
 --  * R, the x value of the first bin
+-- 
 --  * V, the volume of each bin
---  * an RVar providing a random tuple consisting of:
---          - a bin index, uniform over [0,c) :: Int
---          - a uniformly distributed fractional value, from -1 to 1 if not mirrored, from 0 to 1 otherwise.
+-- 
+--  * an RVar providing the 'zGetIU' random tuple
+-- 
 --  * an RVar sampling from the tail (the region where x > R)
-
+-- 
 mkZiggurat_ :: (RealFloat t, Storable t,
                Distribution Uniform t) =>
               Bool
@@ -99,21 +172,21 @@
               -> 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)
+            { 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
+            , zUniform          = uniform
+            , zFunc             = f
+            , zTailDist         = tailDist
+            , zMirror           = m
             }
 
--- |Build the tables to implement the "ziggurat algorithm" devised by 
+-- |Build the tables to implement the \"ziggurat algorithm\" devised by 
 -- Marsaglia & Tang, attempting to automatically compute the R and V
 -- values.
 -- 
--- Arguments are the same as for |mkZigguratRec|, with an additional
+-- 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,
@@ -138,15 +211,31 @@
 -- Arguments:
 -- 
 --  * flag indicating whether to mirror the distribution
---  * the (one-sided antitone) CDF
---  * the inverse of the CDF
---  * the integral of the CDF (definite, from 0)
---  * the estimated volume under the CDF (from 0 to +infinity)
---  * the chunk size (number of bins).  64 seems to perform well in practice.
---  * an RVar providing a random tuple consisting of:
---          - a bin index, uniform over [0,c) :: Int
---          - a uniformly distributed fractional value, from -1 to 1 if not mirrored, from 0 to 1 otherwise.
-
+--
+--  * the (one-sided antitone) PDF, not necessarily normalized
+--
+--  * the inverse of the PDF
+--
+--  * the integral of the PDF (definite, from 0)
+--
+--  * the estimated volume under the PDF (from 0 to +infinity)
+--
+--  * the chunk size (number of bins in each layer).  64 seems to
+--    perform well in practice.
+--
+--  * an RVar providing the 'zGetIU' random tuple
+--
+mkZigguratRec ::
+  (RealFloat t, Storable t,
+   Distribution Uniform t) =>
+  Bool
+  -> (t -> t)
+  -> (t -> t)
+  -> (t -> t)
+  -> 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))
         where
@@ -201,12 +290,18 @@
 -- Arguments:
 -- 
 --  * Number of bins
---  * function (one-sided antitone CDF, not necessarily normalized)
+--
+--  * target function (one-sided antitone PDF, not necessarily normalized)
+--
 --  * function inverse
+--
 --  * function definite integral (from 0 to _)
+--
 --  * estimate of total volume under function (integral from 0 to infinity)
 --
 -- 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)
     where
         c = fromIntegral cInt
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
@@ -1,7 +1,7 @@
 {-
  -      ``Data/Random/Internal/Find''
  -  Utilities for searching fractional domains.  Needs cleanup, testing,
- -  and such.
+ -  and such.  Used for constructing generic ziggurats.
  -}
 
 module Data.Random.Internal.Find where
diff --git a/src/Data/Random/Internal/Fixed.hs b/src/Data/Random/Internal/Fixed.hs
--- a/src/Data/Random/Internal/Fixed.hs
+++ b/src/Data/Random/Internal/Fixed.hs
@@ -1,9 +1,24 @@
+{-# LANGUAGE CPP #-}
 module Data.Random.Internal.Fixed where
 
 import Data.Fixed
 import Unsafe.Coerce
 
+#ifdef base_4_2
+-- So much for backward compatibility through base (>=5) ...
+
 resolutionOf :: HasResolution r => f r -> Integer
+resolutionOf = resolution
+
+resolutionOf2 :: HasResolution r => f (g r) -> Integer
+resolutionOf2 x = resolution (res x)
+    where
+        res :: HasResolution r => f (g r) -> g r
+        res = undefined
+
+#else
+
+resolutionOf :: HasResolution r => f r -> Integer
 resolutionOf x = resolution (res x)
     where
         res :: HasResolution r => f r -> r
@@ -15,6 +30,13 @@
         res :: HasResolution r => f (g r) -> r
         res = undefined
 
+#endif
+
+-- |The 'Fixed' type doesn't expose its constructors, but I need a way to
+-- convert them to and from their raw representation in order to sample
+-- them.  As long as 'Fixed' is a newtype wrapping 'Integer', 'mkFixed' and
+-- 'unMkFixed' as defined here will work.  Both are implemented using 
+-- 'unsafeCoerce'.
 mkFixed :: Integer -> Fixed r
 mkFixed = unsafeCoerce
 
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
@@ -5,7 +5,21 @@
         TemplateHaskell
   #-}
 
-module Data.Random.Internal.TH where
+-- |Template Haskell utility code to replicate instance declarations
+-- to cover large numbers of types.  I'm doing that rather than using
+-- class contexts because most Distribution instances need to cover
+-- multiple classes (such as Enum, Integral and Fractional) and that
+-- can't be done easily because of overlap.  
+-- 
+-- I experimented a bit with a convoluted type-level classification 
+-- scheme, but I think this is simpler and easier to understand.  It 
+-- makes the haddock docs more cluttered because of the combinatorial 
+-- explosion of instances, but overall I think it's just more sane than 
+-- anything else I've come up with yet.
+module Data.Random.Internal.TH
+    ( replicateInstances
+    , integralTypes, realFloatTypes
+    ) where
 
 import Data.Generics
 import Language.Haskell.TH
@@ -13,20 +27,48 @@
 import Data.Word
 import Data.Int
 
+-- |Names of standard 'Integral' types
+integralTypes :: [Name]
 integralTypes = 
     [ ''Int,   ''Integer
     , ''Int8,  ''Int16,  ''Int32,  ''Int64
     , ''Word8, ''Word16, ''Word32, ''Word64
     ]
 
+-- |Names of standard 'RealFloat' types
+realFloatTypes :: [Name]
 realFloatTypes =
     [ ''Float, ''Double ]
 
+-- @replaceName x y@ is a function that will
+-- replace @x@ with @y@ whenever it sees it.  That is:
+--
+-- > replaceName x y x  ==>  y
+-- > replaceName x y z  ==>  z
+--  (@z /= x@)
 replaceName :: Name -> Name -> Name -> Name
 replaceName x y z
     | x == z    = y
     | otherwise = z
 
+-- | @replicateInstances standin types decls@ will take the template-haskell
+-- 'Dec's in @decls@ and substitute every instance of the 'Name' @standin@ with
+-- each 'Name' in @types@, producing one copy of the 'Dec's in @decls@ for every
+-- 'Name' in @types@.
+-- 
+-- For example, 'Data.Random.Distribution.Uniform' has the following bit of TH code:
+-- 
+-- @ $( replicateInstances ''Int integralTypes [d|                                                  @
+-- 
+-- @       instance Distribution Uniform Int   where rvar (Uniform a b) = integralUniform a b       @
+-- 
+-- @       instance CDF Uniform Int            where cdf  (Uniform a b) = integralUniformCDF a b    @
+-- 
+-- @   |])                                                                                          @
+-- 
+-- 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
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
@@ -12,6 +12,11 @@
 import Data.Word
 
 {-# INLINE buildWord #-}
+-- |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
     = unsafePerformIO . allocaBytes 8 $ \p -> do
@@ -26,7 +31,7 @@
         peek (castPtr p)
 
 {-# INLINE wordToFloat #-}
--- |Pack 23 unspecified bits from a 'Word64' into a 'Float' in the range [0,1).
+-- |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'.
 wordToFloat :: Word64 -> Float
 wordToFloat x = (encodeFloat $! toInteger (x .&. 0x007fffff {- 2^23-1 -} )) $ (-23)
@@ -38,7 +43,7 @@
 wordToFloatWithExcess x = (wordToFloat x, x `shiftR` 23)
 
 {-# INLINE wordToDouble #-}
--- |Pack 52 unspecified bits from a 'Word64' into a 'Double' in the range [0,1).
+-- |Pack the low 52 bits from a 'Word64' into a 'Double' in the range [0,1).
 -- Used to convert a 'stdUniform' 'Word64' to a 'stdUniform' 'Double'.
 wordToDouble :: Word64 -> Double
 wordToDouble x = (encodeFloat $! toInteger (x .&. 0x000fffffffffffff {- 2^52-1 -})) $ (-52)
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
@@ -6,12 +6,17 @@
 import qualified System.Random.Shuffle as SRS
 import Control.Monad
 
+-- | A random variable returning an arbitrary element of the given list.
+-- Every element has equal probability of being chosen.  Because it is a
+-- pure 'RVar' it has no memory - that is, it \"draws with replacement.\"
 randomElement :: [a] -> RVar a
 randomElement [] = error "randomElement: empty list!"
 randomElement xs = do
     n <- uniform 0 (length xs - 1)
     return (xs !! n)
 
+-- | A random variable that returns the given list in an arbitrary shuffled
+-- order.  Every ordering of the list has equal probability.
 shuffle :: [a] -> RVar [a]
 shuffle [] = return []
 shuffle xs = do
@@ -19,9 +24,16 @@
     
     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.
 shuffleN :: Int -> [a] -> RVar [a]
 shuffleN 0 xs = return []
-shuffleN (n+1) xs = do
+shuffleN m@(n+1) xs = do
     is <- sequence [uniform 0 i | i <- [n,n-1..1]]
     return (SRS.shuffle xs is)
+
     
diff --git a/src/Data/Random/RVar.hs b/src/Data/Random/RVar.hs
--- a/src/Data/Random/RVar.hs
+++ b/src/Data/Random/RVar.hs
@@ -81,6 +81,17 @@
     getRandomWord   = RVarT $ \k (RVarDict s) -> getRandomWordFrom   s >>= \a -> k a
     getRandomDouble = RVarT $ \k (RVarDict s) -> getRandomDoubleFrom s >>= \a -> k a
 
+-- 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...
+-- instance (Applicative m, Num a) => Num (RVarT m a) where
+--     (+) = liftA2 (+)
+--     (-) = liftA2 (-)
+--     (*) = liftA2 (*)
+--     negate = liftA negate
+--     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 #-}
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
@@ -25,7 +25,7 @@
 instance Distribution d t => Sampleable d m t where
     sampleFrom src d = runRVarT (rvar d) src
 
--- This instance conflicts with the other, but because RVarT is not a Distribution there is no conflict.
+-- This instance overlaps with the other, but because RVarT is not a Distribution there is no conflict.
 instance Lift m n => Sampleable (RVarT m) n t where
     sampleFrom src x = runRVarT x src
 
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,12 +1,14 @@
-{-
- -      ``Data/Random/Source/PureMT''
- -}
 {-# LANGUAGE
     MultiParamTypeClasses,
     FlexibleContexts, FlexibleInstances,
     UndecidableInstances
   #-}
 
+-- |This module provides functions useful for implementing new 'MonadRandom'
+-- and 'RandomSource' instances for state-abstractions containing 'PureMT'
+-- values (the pure pseudorandom generator provided by the
+-- mersenne-random-pure64 package), as well as instances for some common
+-- cases.
 module Data.Random.Source.PureMT where
 
 import Data.Random.Internal.Words
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
@@ -5,6 +5,11 @@
     MultiParamTypeClasses, FlexibleInstances, UndecidableInstances
   #-}
 
+-- |This module provides functions useful for implementing new 'MonadRandom'
+-- and 'RandomSource' instances for state-abstractions containing 'StdGen'
+-- values (the pure pseudorandom generator provided by the System.Random
+-- module in the \"random\" package), as well as instances for some common
+-- cases.
 module Data.Random.Source.StdGen where
 
 import Data.Random.Internal.Words