diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown
new file mode 100644
--- /dev/null
+++ b/CHANGELOG.markdown
@@ -0,0 +1,16 @@
+3.1
+---
+* Added `Chain` mode, which is `Reverse` using a linear tape that doesn't need to be sorted.
+* Added a suite of doctests.
+* Bug fix in `Forward` mode. It was previously yielding incorrect results for anything that used `bind` or `bind'` internally.
+
+3.0
+---
+* Moved the contents of `Numeric.AD.Mode.Mixed` into `Numeric.AD`
+* Split off `Numeric.AD.Variadic` for the variadic combinators
+* Removed the `UU`, `FU`, `UF`, and `FF` type aliases.
+* Stopped exporting the types for `Mode` and `AD` from almost every module. Import `Numeric.AD.Types` if necessary.
+* Renamed `Tensors` to `Jet`
+* Dependency bump to be compatible with ghc 7.4.1 and mtl 2.1
+* More aggressive zero tracking.
+* `diff (**n) 0` for constant n and `diff (0**)` both now yield the correct answer for all modes.
diff --git a/README.markdown b/README.markdown
new file mode 100644
--- /dev/null
+++ b/README.markdown
@@ -0,0 +1,124 @@
+ad
+==
+
+[![Build Status](https://secure.travis-ci.org/ekmett/ad.png?branch=master)](http://travis-ci.org/ekmett/ad)
+
+A package that provides an intuitive API for [Automatic Differentiation](http://en.wikipedia.org/wiki/Automatic_differentiation) (AD) in Haskell. Automatic differentiation provides a means to calculate the derivatives of a function while evaluating it. Unlike numerical methods based on running the program with multiple inputs or symbolic approaches, automatic differentiation typically only decreases performance by a small multiplier.
+
+AD employs the fact that any program `y = F(x)` that computes one or more value does so by composing multiple primitive operations. If the (partial) derivatives of each of those operations is known, then they can be composed to derive the answer for the derivative of the entire program at a point.
+
+This library contains at its core a single implementation that describes how to compute the partial derivatives of a wide array of primitive operations. It then exposes an API that enables a user to safely combine them using standard higher-order functions, just as you would with any other Haskell numerical type.
+
+There are several ways to compose these individual [Jacobian matrices](http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant). We hide the choice used by the API behind an explicit "Mode" type-class and universal quantification. This prevents the end user from exploiting the properties of an individual mode, and thereby potentially violating invariants or [confusing infinitesimals](http://conway.rutgers.edu/~ccshan/wiki/blog/posts/Differentiation/).
+
+Features
+--------
+
+ * Provides forward- and reverse- mode AD combinators with a common API.
+ * Type-level "branding" is used to both prevent the end user from confusing infinitesimals and to limit unsafe access to the implementation details of each mode.
+ * Each mode has a separate module full of combinators, with a consistent look and feel.
+
+Examples
+---------
+
+You can compute derivatives of functions
+
+    Prelude Numeric.AD> diff sin 0 {-# cos 0 #-}
+    1.0
+
+Or both the answer and the derivative of a function:
+
+    Prelude Numeric.AD> diff' (exp . log) 2
+    (2.0,1.0)
+
+You can use a symbolic numeric type, like the one from `simple-reflect` to obtain symbolic derivatives:
+
+    Prelude Debug.SimpleReflect Numeric.AD> diff atanh x
+    recip (1 - x * x) * 1
+
+You can compute gradients for functions that take non-scalar values in the form of a Traversable functor full of AD variables.
+
+    Prelude Numeric.AD Debug.SimpleReflect> grad (\[x,y,z] -> x * sin (x + log y)) [x,y,z]
+    [ 0 + (0 + sin (x + log y) * 1 + 1 * (0 + cos (x + log y) * (0 + x * 1)))
+    , 0 + (0 + recip y * (0 + 1 * (0 + cos (x + log y) * (0 + x * 1))))
+    , 0
+    ]
+
+which one can simplify to:
+
+    [ sin (x + log y) + cos (x + log y) * x, recip y * cos (x + log y) * x, 0 ]
+
+If you need multiple derivatives you can calculate them with `diffs`:
+
+    Prelude Numeric.AD> take 10 $ diffs sin 1
+    [0.8414709848078965,0.5403023058681398,-0.8414709848078965,-0.5403023058681398,0.8414709848078965,0.5403023058681398,-0.8414709848078965,-0.5403023058681398,0.8414709848078965,0.5403023058681398]
+
+or if your function takes multiple inputs, you can use grads, which returns an 'f-branching stream' of derivatives. Somewhat more intuitive answers can be obtained by converting the stream into the
+polymorphically recursive `Tensors` data type. With that we can look at a single 'layer' of the answer at a time:
+
+The answer:
+
+    Prelude Numeric.AD> headJet $ tensors $  grads (\[x,y] -> exp (x * y)) [1,2]
+    7.38905609893065
+
+The gradient:
+
+    Prelude Numeric.AD> headJet $ tailJet $ tensors $  grads (\[x,y] -> exp (x * y)) [1,2]
+    [14.7781121978613,7.38905609893065]
+
+The hessian (n * n matrix of 2nd derivatives)
+
+    Prelude Numeric.AD> headJet $ tailJet $ tailJet $ tensors $  grads (\[x,y] -> exp (x * y)) [1,2]
+    [[29.5562243957226,22.16716829679195],[22.16716829679195,7.38905609893065]]
+
+Or even higher order tensors of derivatives.
+
+    Prelude Numeric.AD> headJet $ tailJet $ tailJet $ tailJet $ tensors $  grads (\[x,y] -> exp (x * y)) [1,2]
+    [[[59.1124487914452,44.3343365935839],[44.3343365935839,14.7781121978613]],[[44.3343365935839,14.7781121978613],[14.7781121978613,7.38905609893065]]]
+
+Note the redundant values caused by the various symmetries in the tensors. The 'ad' library is careful to compute each distinct derivative only once and to share the resulting thunks.
+
+Overview
+--------
+
+### Modules
+
+ * `Numeric.AD` computes using whichever mode or combination thereof is suitable to each individual combinator. This mode is the default, re-exported by `Numeric.AD`
+ * `Numeric.AD.Mode.Forward` provides basic forward-mode AD. It is good for computing simple derivatives.
+ * `Numeric.AD.Mode.Sparse` computes a sparse forward-mode AD tower. It is good for higher derivatives or large numbers of outputs.
+ * `Numeric.AD.Mode.Reverse` computes with reverse-mode AD. It is good for computing a few outputs given many inputs.
+ * `Numeric.AD.Mode.Chain` computes with reverse-mode AD. It is good for computing a few outputs given many inputs, when not using sparks.
+ * `Numeric.AD.Mode.Tower` computes a dense forward-mode AD tower useful for higher derivatives of single input functions.
+
+ * `Numeric.AD.Newton` provides a number of combinators for root finding using Newton's method with quadratic convergence.
+ * `Numeric.AD.Halley` provides a number of combinators for root finding using Halley's method with cubic convergence.
+
+### Combinators
+
+While not every mode can provide all operations, the following basic operations are supported, modified as appropriate by the suffixes below:
+
+ * `grad` computes the gradient (vector of partial derivatives at a given point) of a function.
+ * `jacobian` computes the Jacobian matrix of a function at a point.
+ * `diff` computes the derivative of a function at a point.
+ * `du` computes a directional derivative of a function at a point.
+ * `hessian` computes the Hessian matrix (matrix of second partial derivatives) of a function at a point.
+
+### Combinator Suffixes
+
+The following suffixes alter the meanings of the functions above as follows:
+
+ * `'` also return the answer
+ * `With` lets the user supply a function to blend the input with the output
+ * `F` is a version of the base function lifted to return a `Traversable` (or `Functor`) result
+ * `s` means the function returns all higher derivatives in a list or f-branching `Stream`
+ * `T` means the result is transposed with respect to the traditional formulation (usually to avoid paying for transposing back)
+ * `0` means that the resulting derivative list is padded with 0s at the end.
+
+Contact Information
+-------------------
+
+Contributions and bug reports are welcome!
+
+Please feel free to contact me through github or on the #haskell IRC channel on irc.freenode.net.
+
+-Edward Kmett
diff --git a/ad.cabal b/ad.cabal
--- a/ad.cabal
+++ b/ad.cabal
@@ -1,5 +1,5 @@
 name:         ad
-version:      3.0.1
+version:      3.1.1
 license:      BSD3
 license-File: LICENSE
 copyright:    (c) Edward Kmett 2010-2012,
@@ -11,8 +11,8 @@
 homepage:     http://github.com/ekmett/ad
 bug-reports:  http://github.com/ekmett/ad/issues
 build-type:   Simple
-cabal-version: >= 1.6
-extra-source-files: TODO .travis.yml
+cabal-version: >= 1.8
+extra-source-files: TODO .travis.yml CHANGELOG.markdown README.markdown
 synopsis:     Automatic Differentiation
 description:
     Forward-, reverse- and mixed- mode automatic differentiation combinators with a common API.
@@ -24,8 +24,10 @@
     .
     * @Numeric.AD.Mode.Forward@ provides basic forward-mode AD. It is good for computing simple derivatives.
     .
-    * @Numeric.AD.Mode.Reverse@ uses benign side-effects to compute reverse-mode AD. It is good for computing gradients in one pass.
+    * @Numeric.AD.Mode.Reverse@ uses benign side-effects to compute reverse-mode AD. It is good for computing gradients in one pass. It generates a tree-like tape that needs to be topologically sorted in the end.
     .
+    * @Numeric.AD.Mode.Chain@ uses benign side-effects to compute reverse-mode AD. It is good for computing gradients in one pass. It generates a linear tape using @Data.Reflection@.
+    .
     * @Numeric.AD.Mode.Sparse@ computes a sparse forward-mode AD tower. It is good for higher derivatives or large numbers of outputs.
     .
     * @Numeric.AD.Mode.Tower@ computes a dense forward-mode AD tower useful for higher derivatives of single input functions.
@@ -58,24 +60,6 @@
     * @T@ means the result is transposed with respect to the traditional formulation.
     .
     * @0@ means that the resulting derivative list is padded with 0s at the end.
-    .
-    /Changes since 1.3/:
-    .
-    * Moved the contents of @Numeric.AD.Mode.Mixed@ into @Numeric.AD@
-    .
-    * Split off @Numeric.AD.Variadic@ for the variadic combinators
-    .
-    * Removed the @UU@, @FU@, @UF@, and @FF@ type aliases.
-    .
-    * Stopped exporting the types for @Mode@ and @AD@ from almost every module. Import @Numeric.AD.Types@ if necessary.
-    .
-    * Renamed @Tensors@ to @Jet@
-    .
-    * Dependency bump to be compatible with ghc 7.4.1 and mtl 2.1
-    .
-    * More aggressive zero tracking.
-    .
-    * @diff (**n) 0@ for constant n and @diff (0**)@ both now yield the correct answer for all modes.
 
 source-repository head
   type: git
@@ -107,7 +91,9 @@
     comonad          == 3.0.*,
     containers       >= 0.2 && < 0.6,
     data-reify       >= 0.6 && < 0.7,
-    free             == 3.0.*,
+    free             >= 3.0 && <= 3.2,
+    reflection       >= 1.1.6 && < 1.2,
+    tagged           >= 0.4.2.1 && < 0.5,
     template-haskell >= 2.5 && < 2.9
 
   exposed-modules:
@@ -117,6 +103,7 @@
     Numeric.AD.Newton
     Numeric.AD.Halley
 
+    Numeric.AD.Mode.Chain
     Numeric.AD.Mode.Directed
     Numeric.AD.Mode.Forward
     Numeric.AD.Mode.Reverse
@@ -132,6 +119,8 @@
     Numeric.AD.Internal.Forward
     Numeric.AD.Internal.Tower
     Numeric.AD.Internal.Reverse
+    Numeric.AD.Internal.Var
+    Numeric.AD.Internal.Chain
     Numeric.AD.Internal.Sparse
     Numeric.AD.Internal.Dense
     Numeric.AD.Internal.Composition
@@ -142,3 +131,15 @@
     Numeric.AD.Internal.Identity
 
   ghc-options: -Wall -fspec-constr -fdicts-cheap -O2
+
+-- Verify the results of the examples
+test-suite doctests
+  type:    exitcode-stdio-1.0
+  main-is: doctests.hs
+  build-depends:
+    base == 4.*,
+    directory >= 1.0 && < 1.2,
+    doctest >= 0.8 && <= 0.9,
+    filepath >= 1.3 && < 1.4
+  ghc-options: -Wall -Werror -threaded
+  hs-source-dirs: tests
diff --git a/src/Numeric/AD.hs b/src/Numeric/AD.hs
--- a/src/Numeric/AD.hs
+++ b/src/Numeric/AD.hs
@@ -208,7 +208,7 @@
 hessian :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f (f a)
 hessian f = Sparse.jacobian (grad (decomposeMode . f . fmap composeMode))
 
--- | Compute the order 3 Hessian tensor on a non-scalar-to-non-scalar function using Sparse or Sparse-on-Reverse
+-- | Compute the order 3 Hessian tensor on a non-scalar-to-non-scalar function using 'Sparse' or 'Sparse'-on-'Reverse'
 hessianF :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f (f a))
 hessianF f as
     | big (size as) = decomposeFunctor $ Sparse.jacobian (ComposeFunctor . Reverse.jacobian (fmap decomposeMode . f . fmap composeMode)) as
diff --git a/src/Numeric/AD/Halley.hs b/src/Numeric/AD/Halley.hs
--- a/src/Numeric/AD/Halley.hs
+++ b/src/Numeric/AD/Halley.hs
@@ -36,11 +36,12 @@
 --
 -- Examples:
 --
---  > take 10 $ findZero (\\x->x^2-4) 1  -- converge to 2.0
---
---  > module Data.Complex
---  > take 10 $ findZero ((+1).(^2)) (1 :+ 1)  -- converge to (0 :+ 1)@
+-- >>> take 10 $ findZero (\x->x^2-4) 1
+-- [1.0,1.8571428571428572,1.9997967892704736,1.9999999999994755,2.0]
 --
+-- >>> import Data.Complex
+-- >>> last $ take 10 $ findZero ((+1).(^2)) (1 :+ 1)
+-- 0.0 :+ 1.0
 findZero :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
 findZero f = go
     where
@@ -54,8 +55,7 @@
 -- results.  (Modulo the usual caveats.)
 --
 -- Note: the @take 10 $ inverse sqrt 1 (sqrt 10)@ example that works for Newton's method
--- fails with Halley's method because the preconditions do not hold.
-
+-- fails with Halley's method because the preconditions do not hold!
 inverse :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> a -> [a]
 inverse f x0 y = findZero (\x -> f x - lift y) x0
 {-# INLINE inverse  #-}
@@ -64,7 +64,8 @@
 -- function using Halley's method; its output is a stream of
 -- increasingly accurate results.  (Modulo the usual caveats.)
 --
--- > take 10 $ fixedPoint cos 1 -- converges to 0.7390851332151607
+-- >>> last $ take 10 $ fixedPoint cos 1
+-- 0.7390851332151607
 fixedPoint :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
 fixedPoint f = findZero (\x -> f x - x)
 {-# INLINE fixedPoint #-}
@@ -73,7 +74,8 @@
 -- function using Halley's method; produces a stream of increasingly
 -- accurate results.  (Modulo the usual caveats.)
 --
--- > take 10 $ extremum cos 1 -- convert to 0
+-- >>> take 10 $ extremum cos 1
+-- [1.0,0.29616942658570555,4.59979519460002e-3,1.6220740159042513e-8,0.0]
 extremum :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
 extremum f = findZero (diff (decomposeMode . f . composeMode))
 {-# INLINE extremum #-}
diff --git a/src/Numeric/AD/Internal/Chain.hs b/src/Numeric/AD/Internal/Chain.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/AD/Internal/Chain.hs
@@ -0,0 +1,226 @@
+{-# LANGUAGE Rank2Types, TypeFamilies, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, FlexibleContexts, TemplateHaskell, UndecidableInstances, DeriveDataTypeable, GADTs, ScopedTypeVariables #-}
+-- {-# OPTIONS_HADDOCK hide, prune #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD.Internal.Chain
+-- Copyright   :  (c) Edward Kmett 2012
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only
+--
+-- Reverse-Mode Automatic Differentiation using a single tape.
+--
+-- This version uses @Data.Reflection@ to find and update the tape
+--
+-- This is asymptotically faster than using @Reverse@, which
+-- is forced to reify and topologically sort the graph, but it requires
+-- a fairly expensive rendezvous during construction.
+--
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Internal.Chain
+    ( Chain(..)
+    , Tape(..)
+    , Head(..)
+    , Cells(..)
+    , reifyTape
+    , partials
+    , partialArrayOf
+    , partialMapOf
+    , derivativeOf
+    , derivativeOf'
+    ) where
+
+import Control.Monad.ST
+import Data.Array.ST
+import Data.Array
+import Data.Array.Unsafe as Unsafe
+import Data.IORef
+import Data.IntMap (IntMap, fromDistinctAscList)
+import Data.Proxy
+import Data.Reflection
+import Data.Typeable
+import Language.Haskell.TH hiding (reify)
+import Numeric.AD.Internal.Types
+import Numeric.AD.Internal.Classes
+import Numeric.AD.Internal.Identity
+import Numeric.AD.Internal.Var
+import Prelude hiding (mapM)
+import System.IO.Unsafe (unsafePerformIO)
+import Unsafe.Coerce
+
+-- evil untyped tape
+data Cells where
+  Nil    :: Cells
+  Unary  :: {-# UNPACK #-} !Int -> a -> Cells -> Cells
+  Binary :: {-# UNPACK #-} !Int -> {-# UNPACK #-} !Int -> a -> a -> Cells -> Cells
+
+dropCells :: Int -> Cells -> Cells
+dropCells 0 xs = xs
+dropCells _ Nil = Nil
+dropCells n (Unary _ _ xs)      = (dropCells $! n - 1) xs
+dropCells n (Binary _ _ _ _ xs) = (dropCells $! n - 1) xs
+
+data Head = Head {-# UNPACK #-} !Int Cells
+
+newtype Tape = Tape { getTape :: IORef Head }
+
+un :: Int -> a -> Head -> (Head, Int)
+un i di (Head r t) = h `seq` r' `seq` (h, r') where
+  r' = r + 1
+  h = Head r' (Unary i di t)
+{-# INLINE un #-}
+
+bin :: Int -> Int -> a -> a -> Head -> (Head, Int)
+bin i j di dj (Head r t) = h `seq` r' `seq` (h, r') where
+  r' = r + 1
+  h = Head r' (Binary i j di dj t)
+{-# INLINE bin #-}
+
+modifyTape :: Reifies s Tape => p s -> (Head -> (Head, r)) -> IO r
+modifyTape p = atomicModifyIORef (getTape (reflect p))
+{-# INLINE modifyTape #-}
+
+-- | This is used to create a new entry on the chain given a unary function, its derivative with respect to its input,
+-- the variable ID of its input, and the value of its input. Used by 'unary' and 'binary' internally.
+unarily :: forall s a. Reifies s Tape => (a -> a) -> a -> Int -> a -> Chain s a
+unarily f di i b = Chain (unsafePerformIO (modifyTape (Proxy :: Proxy s) (un i di))) $! f b
+{-# INLINE unarily #-}
+
+-- | This is used to create a new entry on the chain given a binary function, its derivatives with respect to its inputs,
+-- their variable IDs and values. Used by 'binary' internally.
+binarily :: forall s a. Reifies s Tape => (a -> a -> a) -> a -> a -> Int -> a -> Int -> a -> Chain s a
+binarily f di dj i b j c = Chain (unsafePerformIO (modifyTape (Proxy :: Proxy s) (bin i j di dj))) $! f b c
+{-# INLINE binarily #-}
+
+data Chain s a where
+  Zero :: Chain s a
+  Lift :: a -> Chain s a
+  Chain :: {-# UNPACK #-} !Int -> a -> Chain s a
+  deriving (Show, Typeable)
+
+instance (Reifies s Tape, Lifted (Chain s)) => Mode (Chain s) where
+  isKnownZero Zero = True
+  isKnownZero _    = False
+
+  isKnownConstant Chain{} = False
+  isKnownConstant _ = True
+
+  lift = Lift
+  zero = Zero
+  (<+>)  = binary (+) one one
+  a *^ b = lift1 (a *) (\_ -> lift a) b
+  a ^* b = lift1 (* b) (\_ -> lift b) a
+  a ^/ b = lift1 (/ b) (\_ -> lift (recip b)) a
+
+  Zero <**> y      = lift (0 ** primal y)
+  _    <**> Zero   = lift 1
+  x    <**> Lift y = lift1 (**y) (\z -> (y *^ z ** Id (y-1))) x
+  x    <**> y      = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y
+
+instance Primal (Chain s) where
+    primal Zero = 0
+    primal (Lift a) = a
+    primal (Chain _ a) = a
+
+instance (Reifies s Tape, Lifted (Chain s)) => Jacobian (Chain s) where
+    type D (Chain s) = Id
+
+    unary f _         (Zero)   = Lift (f 0)
+    unary f _         (Lift a) = Lift (f a)
+    unary f (Id dadi) (Chain i b) = unarily f dadi i b
+
+    lift1 f df b = unary f (df (Id pb)) b
+        where pb = primal b
+
+    lift1_ f df b = unary (const a) (df (Id a) (Id pb)) b
+        where pb = primal b
+              a = f pb
+
+    binary f _         _         Zero     Zero     = Lift (f 0 0)
+    binary f _         _         Zero     (Lift c) = Lift (f 0 c)
+    binary f _         _         (Lift b) Zero     = Lift (f b 0)
+    binary f _         _         (Lift b) (Lift c) = Lift (f b c)
+
+    binary f _         (Id dadc) Zero        (Chain i c) = unarily (f 0) dadc i c
+    binary f _         (Id dadc) (Lift b)    (Chain i c) = unarily (f b) dadc i c
+    binary f (Id dadb) _         (Chain i b) Zero        = unarily (`f` 0) dadb i b
+    binary f (Id dadb) _         (Chain i b) (Lift c)    = unarily (`f` c) dadb i b
+    binary f (Id dadb) (Id dadc) (Chain i b) (Chain j c) = binarily f dadb dadc i b j c
+
+    lift2 f df b c = binary f dadb dadc b c
+        where (dadb, dadc) = df (Id (primal b)) (Id (primal c))
+
+    lift2_ f df b c = binary (\_ _ -> a) dadb dadc b c
+        where
+            pb = primal b
+            pc = primal c
+            a = f pb pc
+            (dadb, dadc) = df (Id a) (Id pb) (Id pc)
+
+let s = varT (mkName "s") in
+  deriveLifted (classP ''Reifies [s, conT ''Tape] :) (conT ''Chain `appT` s)
+
+-- | Helper that extracts the derivative of a chain when the chain was constructed with one variable.
+derivativeOf :: (Reifies s Tape, Num a) => Proxy s -> AD (Chain s) a -> a
+derivativeOf _ = sum . partials
+{-# INLINE derivativeOf #-}
+
+-- | Helper that extracts both the primal and derivative of a chain when the chain was constructed with one variable.
+derivativeOf' :: (Reifies s Tape, Num a) => Proxy s -> AD (Chain s) a -> (a, a)
+derivativeOf' p r = (primal r, derivativeOf p r)
+{-# INLINE derivativeOf' #-}
+
+-- | Used internally to push sensitivities down the chain.
+backPropagate :: Num a => Int -> Cells -> STArray s Int a -> ST s Int
+backPropagate k Nil _ = return k
+backPropagate k (Unary i g xs) ss = do
+  da <- readArray ss k
+  db <- readArray ss i
+  writeArray ss i $! db + unsafeCoerce g*da
+  (backPropagate $! k - 1) xs ss
+backPropagate k (Binary i j g h xs) ss = do
+  da <- readArray ss k
+  db <- readArray ss i
+  writeArray ss i $! db + unsafeCoerce g*da
+  dc <- readArray ss j
+  writeArray ss j $! dc + unsafeCoerce h*da
+  (backPropagate $! k - 1) xs ss
+
+-- | Extract the partials from the current chain for a given AD variable.
+{-# SPECIALIZE partials :: Reifies s Tape => AD (Chain s) Double -> [Double] #-}
+partials :: forall s a. (Reifies s Tape, Num a) => AD (Chain s) a -> [a]
+partials (AD Zero)        = []
+partials (AD (Lift _))    = []
+partials (AD (Chain k _)) = map (sensitivities !) [0..vs] where
+   Head n t = unsafePerformIO $ readIORef (getTape (reflect (Proxy :: Proxy s)))
+   tk = dropCells (n - k) t
+   (vs,sensitivities) = runST $ do
+     ss <- newArray (0, k) 0
+     writeArray ss k 1
+     v <- backPropagate k tk ss
+     as <- Unsafe.unsafeFreeze ss
+     return (v, as)
+
+-- | Return an 'Array' of 'partials' given bounds for the variable IDs.
+partialArrayOf :: (Reifies s Tape, Num a) => Proxy s -> (Int, Int) -> AD (Chain s) a -> Array Int a
+partialArrayOf _ vbounds = accumArray (+) 0 vbounds . zip [0..] . partials
+{-# INLINE partialArrayOf #-}
+
+-- | Return an 'IntMap' of sparse partials
+partialMapOf :: (Reifies s Tape, Num a) => Proxy s -> AD (Chain s) a -> IntMap a
+partialMapOf _ = fromDistinctAscList . zip [0..] . partials
+{-# INLINE partialMapOf #-}
+
+-- | Construct a tape that starts with @n@ variables.
+reifyTape :: Int -> (forall s. Reifies s Tape => Proxy s -> r) -> r
+reifyTape vs k = unsafePerformIO $ do
+  h <- newIORef (Head vs Nil)
+  return (reify (Tape h) k)
+{-# NOINLINE reifyTape #-}
+
+instance Var (Chain s) where
+    var a v = Chain v a
+    varId (Chain v _) = v
+    varId _ = error "varId: not a Var"
diff --git a/src/Numeric/AD/Internal/Classes.hs b/src/Numeric/AD/Internal/Classes.hs
--- a/src/Numeric/AD/Internal/Classes.hs
+++ b/src/Numeric/AD/Internal/Classes.hs
@@ -29,7 +29,7 @@
 import Control.Applicative hiding ((<**>))
 import Data.Char
 import Language.Haskell.TH
-import Numeric.AD.Internal.Combinators (on)
+import Data.Function (on)
 
 infixr 8 **!, <**>
 infixl 7 *!, /!, ^*, *^, ^/
diff --git a/src/Numeric/AD/Internal/Combinators.hs b/src/Numeric/AD/Internal/Combinators.hs
--- a/src/Numeric/AD/Internal/Combinators.hs
+++ b/src/Numeric/AD/Internal/Combinators.hs
@@ -8,21 +8,20 @@
 -- Stability   :  experimental
 -- Portability :  GHC only
 --
+-- Combinators used internally by @Numeric.AD@
 -----------------------------------------------------------------------------
 module Numeric.AD.Internal.Combinators
     ( zipWithT
     , zipWithDefaultT
-    , on
     ) where
 
 import Data.Traversable (Traversable, mapAccumL)
 import Data.Foldable (Foldable, toList)
 
-on :: (a -> a -> b) -> (c -> a) -> c -> c -> b
-on f g a b = f (g a) (g b)
-
+-- | Zip a @'Foldable' f@ with a @'Traversable' g@ assuming @f@ has at least as many entries as @g@.
 zipWithT :: (Foldable f, Traversable g) => (a -> b -> c) -> f a -> g b -> g c
 zipWithT f as = snd . mapAccumL (\(a:as') b -> (as', f a b)) (toList as)
 
+-- | Zip a @'Foldable' f@ with a @'Traversable' g@ assuming @f@, using a default value after @f@ is exhausted.
 zipWithDefaultT :: (Foldable f, Traversable g) => a -> (a -> b -> c) -> f a -> g b -> g c
 zipWithDefaultT z f as = zipWithT f (toList as ++ repeat z)
diff --git a/src/Numeric/AD/Internal/Forward.hs b/src/Numeric/AD/Internal/Forward.hs
--- a/src/Numeric/AD/Internal/Forward.hs
+++ b/src/Numeric/AD/Internal/Forward.hs
@@ -37,12 +37,14 @@
 import Numeric.AD.Internal.Classes
 import Numeric.AD.Internal.Identity
 
+-- | 'Forward' mode AD.
 data Forward a
   = Forward !a a
   | Lift !a
   | Zero
   deriving (Show, Data, Typeable)
 
+-- | Calculate the 'tangent' using forward mode AD.
 tangent :: Num a => AD Forward a -> a
 tangent (AD (Forward _ da)) = da
 tangent _ = 0
@@ -166,13 +168,13 @@
 bind f as = snd $ mapAccumL outer (0 :: Int) as
     where
         outer !i _ = (i + 1, f $ snd $ mapAccumL (inner i) 0 as)
-        inner !i !j a = (j + 1, if i == j then bundle a 1 else AD Zero)
+        inner !i !j a = (j + 1, if i == j then bundle a 1 else lift a)
 
 bind' :: (Traversable f, Num a) => (f (AD Forward a) -> b) -> f a -> (b, f b)
 bind' f as = dropIx $ mapAccumL outer (0 :: Int, b0) as
     where
         outer (!i, _) _ = let b = f $ snd $ mapAccumL (inner i) (0 :: Int) as in ((i + 1, b), b)
-        inner !i !j a = (j + 1, if i == j then bundle a 1 else AD Zero)
+        inner !i !j a = (j + 1, if i == j then bundle a 1 else lift a)
         b0 = f (lift <$> as)
         dropIx ((_,b),bs) = (b,bs)
 
diff --git a/src/Numeric/AD/Internal/Jet.hs b/src/Numeric/AD/Internal/Jet.hs
--- a/src/Numeric/AD/Internal/Jet.hs
+++ b/src/Numeric/AD/Internal/Jet.hs
@@ -1,4 +1,4 @@
-{-# LANGUAGE TypeOperators, TemplateHaskell, ScopedTypeVariables, FlexibleContexts #-}
+{-# LANGUAGE CPP, TypeOperators, TemplateHaskell, ScopedTypeVariables, FlexibleContexts #-}
 {-# OPTIONS_HADDOCK hide #-}
 -----------------------------------------------------------------------------
 -- |
@@ -10,7 +10,6 @@
 -- Portability :  GHC only
 --
 -----------------------------------------------------------------------------
-
 module Numeric.AD.Internal.Jet
     ( Jet(..)
     , headJet
@@ -35,9 +34,14 @@
 
 infixl 3 :-
 
--- | A jet is a tower of all (higher order) partial derivatives of a function
+-- | A 'Jet' is a tower of all (higher order) partial derivatives of a function
+--
+-- At each step, a @'Jet' f@ is wrapped in another layer worth of @f@.
+--
+-- > a :- f a :- f (f a) :- f (f (f a)) :- ...
 data Jet f a = a :- Jet f (f a)
 
+-- | Used to sidestep the need for UndecidableInstances.
 newtype Showable = Showable (Int -> String -> String)
 
 instance Show Showable where
@@ -62,14 +66,17 @@
 instance Traversable f => Traversable (Jet f) where
     traverse f (a :- as) = (:-) <$> f a <*> traverse (traverse f) as
 
+-- | Take the tail of a 'Jet'.
 tailJet :: Jet f a -> Jet f (f a)
 tailJet (_ :- as) = as
 {-# INLINE tailJet #-}
 
+-- | Take the head of a 'Jet'.
 headJet :: Jet f a -> a
 headJet (a :- _) = a
 {-# INLINE headJet #-}
 
+-- | Construct a 'Jet' by unzipping the layers of a 'Cofree' 'Comonad'.
 jet :: Functor f => Cofree f a -> Jet f a
 jet (a :< as) = a :- dist (jet <$> as)
     where
diff --git a/src/Numeric/AD/Internal/Reverse.hs b/src/Numeric/AD/Internal/Reverse.hs
--- a/src/Numeric/AD/Internal/Reverse.hs
+++ b/src/Numeric/AD/Internal/Reverse.hs
@@ -26,12 +26,6 @@
     , partialMap
     , derivative
     , derivative'
-    , Var(..)
-    , bind
-    , unbind
-    , unbindMap
-    , unbindWith
-    , unbindMapWithDefault
     , vgrad, vgrad'
     , Grad(..)
     ) where
@@ -40,16 +34,13 @@
 import Control.Applicative (Applicative(..),(<$>))
 import Control.Monad.ST
 import Control.Monad (forM_)
-import Data.List (foldl', delete)
+import Data.List (foldl')
 import Data.Array.ST
 import Data.Array
-import Data.IntMap (IntMap, fromListWith, findWithDefault, fromAscList, 
-                    updateLookupWithKey)
-import qualified Data.IntSet as IS
-import Data.Graph (graphFromEdges', Vertex, vertices, edges, transposeG, Graph)
+import Data.IntMap (IntMap, fromListWith)
+import Data.Graph (Vertex, transposeG, Graph)
 import Data.Reify (reifyGraph, MuRef(..))
 import qualified Data.Reify.Graph as Reified
-import Data.Traversable (Traversable, mapM)
 import System.IO.Unsafe (unsafePerformIO)
 import Language.Haskell.TH
 import Data.Data (Data)
@@ -57,6 +48,7 @@
 import Numeric.AD.Internal.Types
 import Numeric.AD.Internal.Classes
 import Numeric.AD.Internal.Identity
+import Numeric.AD.Internal.Var
 
 -- | A @Tape@ records the information needed back propagate from the output to each input during 'Reverse' 'Mode' AD.
 data Tape a t
@@ -82,6 +74,13 @@
     mapDeRef f (Reverse (Unary a dadb b)) = Unary a dadb <$> f b
 
 instance Lifted Reverse => Mode Reverse where
+    isKnownZero (Reverse Zero) = True
+    isKnownZero _    = False
+
+    isKnownConstant (Reverse Zero) = True
+    isKnownConstant (Reverse (Lift a)) = True
+    isKnownConstant _ = False
+
     lift a = Reverse (Lift a)
     zero   = Reverse Zero
     (<+>)  = binary (+) one one
@@ -162,40 +161,50 @@
             _ -> return ()
     where
         (node, i, _) = vmap v
-
         -- this isn't _quite_ right, as it should allow negative zeros to multiply through
 
 topSortAcyclic :: Graph -> [Vertex]
-topSortAcyclic g = go (fromAscList . assocs $ transposeG g) starters
-  where starters = IS.toList $ foldl' (flip IS.delete)
-                                      (IS.fromList $ vertices g)
-                                      (map snd $ edges g)
-        go _ [] = []
-        go g' (n:ns) = let (g'',ns') = foldl' (uncurry (prune n)) (g',[]) (g!n)
-                       in n : go g'' (ns'++ns)
-        prune n g' acc m = let f _ = Just . delete n
-                               (Just ns, g'') = updateLookupWithKey f m g'
-                           in g'' `seq` (g'', if null (tail ns) then m:acc else acc)
-
+topSortAcyclic g = reverse $ runST $ do
+    del <- newArray (bounds g) False :: ST s (STUArray s Int Bool)
+    let tg = transposeG g
+        starters = [ n | (n, []) <- assocs tg ]
+        loop [] rs = return rs
+        loop (n:ns) rs = do
+            writeArray del n True
+            let add [] = return ns
+                add (m:ms) = do
+                    b <- ok (tg!m)
+                    ms' <- add ms
+                    if b then return (m:ms') else return ms'
+                ok [] = return True
+                ok (x:xs) = do b <- readArray del x; if b then ok xs else return False
+            ns' <- add (g!n)
+            loop ns' (n : rs)
+    loop starters []
 
 -- | This returns a list of contributions to the partials.
 -- The variable ids returned in the list are likely /not/ unique!
-partials :: Num a => AD Reverse a -> [(Int, a)]
-partials (AD tape) = [ (ident, sensitivities ! ix) | (ix, Var _ ident) <- xs ]
+{-# SPECIALIZE partials :: AD Reverse Double -> [(Int, Double)] #-}
+partials :: forall a . Num a => AD Reverse a -> [(Int, a)]
+partials (AD tape) = [ let v = sensitivities ! ix in seq v (ident, v) | (ix, Var _ ident) <- xs ]
     where
         Reified.Graph xs start = unsafePerformIO $ reifyGraph tape
-        (g, vmap) = graphFromEdges' (edgeSet <$> filter nonConst xs)
+        g = array xsBounds [ (i, successors t) | (i, t) <- xs ]
+        vertexMap = array xsBounds xs
+        vmap i = (vertexMap ! i, i, [])
+        xsBounds = sbounds xs
+
         sensitivities = runSTArray $ do
-            ss <- newArray (sbounds xs) 0
+            ss <- newArray xsBounds 0
             writeArray ss start 1
             forM_ (topSortAcyclic g) $
                 backPropagate vmap ss
             return ss
+
         sbounds ((a,_):as) = foldl' (\(lo,hi) (b,_) -> let lo' = min lo b; hi' = max hi b in lo' `seq` hi' `seq` (lo', hi')) (a,a) as
         sbounds _ = undefined -- the graph can't be empty, it contains the output node!
-        edgeSet (i, t) = (t, i, successors t)
-        nonConst (_, Lift{}) = False
-        nonConst _ = True
+
+        successors :: Tape a t -> [t]
         successors (Unary _ _ b) = [b]
         successors (Binary _ _ _ b c) = [b,c]
         successors _ = []
@@ -216,37 +225,10 @@
     return a = S (\s -> (a,s))
     S g >>= f = S (\s -> let (a,s') = g s in runS (f a) s')
 
--- | Used to mark variables for inspection during the reverse pass
-class Primal v => Var v where
-    var   :: a -> Int -> v a
-    varId :: v a -> Int
-
 instance Var Reverse where
     var a v = Reverse (Var a v)
     varId (Reverse (Var _ v)) = v
     varId _ = error "varId: not a Var"
-
-instance Var (AD Reverse) where
-    var a v = AD (var a v)
-    varId (AD v) = varId v
-
-bind :: (Traversable f, Var v) => f a -> (f (v a), (Int,Int))
-bind xs = (r,(0,hi))
-    where
-        (r,hi) = runS (mapM freshVar xs) 0
-        freshVar a = S (\s -> let s' = s + 1 in s' `seq` (var a s, s'))
-
-unbind :: (Functor f, Var v)  => f (v a) -> Array Int a -> f a
-unbind xs ys = fmap (\v -> ys ! varId v) xs
-
-unbindWith :: (Functor f, Var v, Num a) => (a -> b -> c) -> f (v a) -> Array Int b -> f c
-unbindWith f xs ys = fmap (\v -> f (primal v) (ys ! varId v)) xs
-
-unbindMap :: (Functor f, Var v, Num a) => f (v a) -> IntMap a -> f a
-unbindMap xs ys = fmap (\v -> findWithDefault 0 (varId v) ys) xs
-
-unbindMapWithDefault :: (Functor f, Var v, Num a) => b -> (a -> b -> c) -> f (v a) -> IntMap b -> f c
-unbindMapWithDefault z f xs ys = fmap (\v -> f (primal v) $ findWithDefault z (varId v) ys) xs
 
 class Num a => Grad i o o' a | i -> a o o', o -> a i o', o' -> a i o where
     pack :: i -> [AD Reverse a] -> AD Reverse a
diff --git a/src/Numeric/AD/Internal/Var.hs b/src/Numeric/AD/Internal/Var.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/AD/Internal/Var.hs
@@ -0,0 +1,80 @@
+-- {-# OPTIONS_HADDOCK hide, prune #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD.Internal.Var
+-- Copyright   :  (c) Edward Kmett 2012
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only
+--
+-- Reverse-Mode Automatic Differentiation using a single tape.
+--
+-- This version uses @Data.Reflection@ to update a single tape.
+--
+-- This is asymptotically faster than using @Reverse@, which
+-- is forced to reify and topologically sort the graph, but it is
+-- less friendly to the use of sparks.
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Internal.Var
+    ( Var(..)
+    , bind
+    , unbind
+    , unbindMap
+    , unbindWith
+    , unbindMapWithDefault
+    , Variable(..)
+    , vary
+    ) where
+
+import Prelude hiding (mapM)
+import Data.Array
+import Data.IntMap (IntMap, findWithDefault)
+import Data.Traversable (Traversable, mapM)
+import Numeric.AD.Internal.Types
+import Numeric.AD.Internal.Classes
+
+-- | Used to mark variables for inspection during the reverse pass
+class Primal v => Var v where
+    var   :: a -> Int -> v a
+    varId :: v a -> Int
+
+instance Var f => Var (AD f) where
+    var a v = AD (var a v)
+    varId (AD v) = varId v
+
+-- A simple fresh variable supply monad
+newtype S a = S { runS :: Int -> (a,Int) }
+instance Monad S where
+    return a = S (\s -> (a,s))
+    S g >>= f = S (\s -> let (a,s') = g s in runS (f a) s')
+
+bind :: (Traversable f, Var v) => f a -> (f (v a), (Int,Int))
+bind xs = (r,(0,hi)) where
+  (r,hi) = runS (mapM freshVar xs) 0
+  freshVar a = S (\s -> let s' = s + 1 in s' `seq` (var a s, s'))
+
+unbind :: (Functor f, Var v)  => f (v a) -> Array Int a -> f a
+unbind xs ys = fmap (\v -> ys ! varId v) xs
+
+unbindWith :: (Functor f, Var v, Num a) => (a -> b -> c) -> f (v a) -> Array Int b -> f c
+unbindWith f xs ys = fmap (\v -> f (primal v) (ys ! varId v)) xs
+
+unbindMap :: (Functor f, Var v, Num a) => f (v a) -> IntMap a -> f a
+unbindMap xs ys = fmap (\v -> findWithDefault 0 (varId v) ys) xs
+
+unbindMapWithDefault :: (Functor f, Var v, Num a) => b -> (a -> b -> c) -> f (v a) -> IntMap b -> f c
+unbindMapWithDefault z f xs ys = fmap (\v -> f (primal v) $ findWithDefault z (varId v) ys) xs
+
+data Variable a = Variable a {-# UNPACK #-} !Int
+
+instance Var Variable where
+  var = Variable
+  varId (Variable _ i) = i
+
+instance Primal Variable where
+  primal (Variable a _) = a
+
+vary :: Var f => Variable a -> f a
+vary (Variable a i) = var a i
diff --git a/src/Numeric/AD/Mode/Chain.hs b/src/Numeric/AD/Mode/Chain.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/AD/Mode/Chain.hs
@@ -0,0 +1,194 @@
+{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, UndecidableInstances, ScopedTypeVariables #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD.Mode.Chain
+-- Copyright   :  (c) Edward Kmett 2010
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only
+--
+-- Reverse Automatic Differentiation using Data.Reflection
+--
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Mode.Chain
+    (
+    -- * Gradient
+      grad
+    , grad'
+    , gradWith
+    , gradWith'
+
+    -- * Jacobian
+    , jacobian
+    , jacobian'
+    , jacobianWith
+    , jacobianWith'
+
+    -- * Hessian
+    , hessian
+    , hessianF
+
+    -- * Derivatives
+    , diff
+    , diff'
+    , diffF
+    , diffF'
+    ) where
+
+import Control.Applicative ((<$>))
+import Data.Traversable (Traversable)
+
+import Numeric.AD.Types
+import Numeric.AD.Internal.Classes
+import Numeric.AD.Internal.Composition
+import Numeric.AD.Internal.Chain
+import Numeric.AD.Internal.Var
+
+-- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with reverse-mode AD in a single pass.
+--
+--
+-- >>> grad (\[x,y,z] -> x*y+z) [1,2,3]
+-- [2,1,1]
+grad :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a
+grad f as = reifyTape (snd bds) $ \p -> unbind vs $! partialArrayOf p bds $! f $ vary <$> vs
+  where (vs, bds) = bind as
+{-# INLINE grad #-}
+
+-- | The 'grad'' function calculates the result and gradient of a non-scalar-to-scalar function with reverse-mode AD in a single pass.
+--
+-- >>> grad' (\[x,y,z] -> x*y+z) [1,2,3]
+-- (5,[2,1,1])
+grad' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)
+grad' f as = reifyTape (snd bds) $ \p ->
+  let r = f (fmap vary vs) in (primal r, unbind vs $! partialArrayOf p bds $! r)
+  where (vs, bds) = bind as
+{-# INLINE grad' #-}
+
+-- | @'grad' g f@ function calculates the gradient of a non-scalar-to-scalar function @f@ with reverse-mode AD in a single pass.
+-- The gradient is combined element-wise with the argument using the function @g@.
+--
+-- @
+-- 'grad' == 'gradWith' (\_ dx -> dx)
+-- 'id' == 'gradWith' 'const'
+-- @
+gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f b
+gradWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs $! partialArrayOf p bds $! f $ vary <$> vs
+  where (vs,bds) = bind as
+{-# INLINE gradWith #-}
+
+-- | @'grad'' g f@ calculates the result and gradient of a non-scalar-to-scalar function @f@ with reverse-mode AD in a single pass
+-- the gradient is combined element-wise with the argument using the function @g@.
+--
+-- @
+-- 'grad'' == 'gradWith'' (\_ dx -> dx)
+-- @
+gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f b)
+gradWith' g f as = reifyTape (snd bds) $ \p ->
+   let r = f (fmap vary vs) in (primal r, unbindWith g vs $! partialArrayOf p bds $! r)
+    where (vs, bds) = bind as
+{-# INLINE gradWith' #-}
+
+-- | The 'jacobian' function calculates the jacobian of a non-scalar-to-non-scalar function with reverse AD lazily in @m@ passes for @m@ outputs.
+--
+-- >>> jacobian (\[x,y] -> [y,x,x*y]) [2,1]
+-- [[0,1],[1,0],[1,2]]
+jacobian :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)
+jacobian f as = reifyTape (snd bds) $ \p -> unbind vs . partialArrayOf p bds <$> f (fmap vary vs)
+  where (vs, bds) = bind as
+{-# INLINE jacobian #-}
+
+-- | The 'jacobian'' function calculates both the result and the Jacobian of a nonscalar-to-nonscalar function, using @m@ invocations of reverse AD,
+-- where @m@ is the output dimensionality. Applying @fmap snd@ to the result will recover the result of 'jacobian'
+-- | An alias for 'gradF''
+--
+-- >>> jacobian' (\[x,y] -> [y,x,x*y]) [2,1]
+-- [(1,[0,1]),(2,[1,0]),(2,[1,2])]
+jacobian' :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)
+jacobian' f as = reifyTape (snd bds) $ \p ->
+  let row a = (primal a, unbind vs $! partialArrayOf p bds $! a)
+  in row <$> f (vary <$> vs)
+  where (vs, bds) = bind as
+{-# INLINE jacobian' #-}
+
+-- | 'jacobianWith g f' calculates the Jacobian of a non-scalar-to-non-scalar function @f@ with reverse AD lazily in @m@ passes for @m@ outputs.
+--
+-- Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@.
+--
+-- @
+-- 'jacobian' == 'jacobianWith' (\_ dx -> dx)
+-- 'jacobianWith' 'const' == (\f x -> 'const' x '<$>' f x)
+-- @
+jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f b)
+jacobianWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs . partialArrayOf p bds <$> f (fmap vary vs) where
+    (vs, bds) = bind as
+{-# INLINE jacobianWith #-}
+
+-- | 'jacobianWith' g f' calculates both the result and the Jacobian of a nonscalar-to-nonscalar function @f@, using @m@ invocations of reverse AD,
+-- where @m@ is the output dimensionality. Applying @fmap snd@ to the result will recover the result of 'jacobianWith'
+--
+-- Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@.
+--
+-- @'jacobian'' == 'jacobianWith'' (\_ dx -> dx)@
+--
+jacobianWith' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f b)
+jacobianWith' g f as = reifyTape (snd bds) $ \p ->
+  let row a = (primal a, unbindWith g vs $! partialArrayOf p bds $! a)
+  in row <$> f (vary <$> vs)
+  where (vs, bds) = bind as
+{-# INLINE jacobianWith' #-}
+
+-- | Compute the derivative of a function.
+--
+-- >>> diff sin 0
+-- 1.0
+diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a
+diff f a = reifyTape 1 $ \p -> derivativeOf p $! f (var a 0)
+{-# INLINE diff #-}
+
+-- | The 'diff'' function calculates the result and derivative, as a pair, of a scalar-to-scalar function.
+--
+-- >>> diff' sin 0
+-- (0.0,1.0)
+--
+-- >>> diff' exp 0
+-- (1.0,1.0)
+diff' :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)
+diff' f a = reifyTape 1 $ \p -> derivativeOf' p $! f (var a 0)
+{-# INLINE diff' #-}
+
+-- | Compute the derivatives of each result of a scalar-to-vector function with regards to its input.
+--
+-- >>> diffF (\a -> [sin a, cos a]) 0
+-- [1.0,0.0]
+--
+diffF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a
+diffF f a = reifyTape 1 $ \p -> derivativeOf p <$> f (var a 0)
+{-# INLINE diffF #-}
+
+-- | Compute the derivatives of each result of a scalar-to-vector function with regards to its input along with the answer.
+--
+-- >>> diffF' (\a -> [sin a, cos a]) 0
+-- [(0.0,1.0),(1.0,0.0)]
+diffF' :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f (a, a)
+diffF' f a = reifyTape 1 $ \p -> derivativeOf' p <$> f (var a 0)
+{-# INLINE diffF' #-}
+
+-- | Compute the hessian via the jacobian of the gradient. gradient is computed in reverse mode and then the jacobian is computed in reverse mode.
+--
+-- However, since the @'grad' f :: f a -> f a@ is square this is not as fast as using the forward-mode Jacobian of a reverse mode gradient provided by 'Numeric.AD.hessian'.
+--
+-- >>> hessian (\[x,y] -> x*y) [1,2]
+-- [[0,1],[1,0]]
+hessian :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f (f a)
+hessian f = jacobian (grad (decomposeMode . f . fmap composeMode))
+
+-- | Compute the order 3 Hessian tensor on a non-scalar-to-non-scalar function via the reverse-mode Jacobian of the reverse-mode Jacobian of the function.
+--
+-- Less efficient than 'Numeric.AD.Mode.Mixed.hessianF'.
+--
+-- >>> hessianF (\[x,y] -> [x*y,x+y,exp x*cos y]) [1,2]
+-- [[[0.0,1.0],[1.0,0.0]],[[0.0,0.0],[0.0,0.0]],[[-1.1312043837568135,-2.4717266720048188],[-2.4717266720048188,1.1312043837568135]]]
+hessianF :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f (f a))
+hessianF f = decomposeFunctor . jacobian (ComposeFunctor . jacobian (fmap decomposeMode . f . fmap composeMode))
diff --git a/src/Numeric/AD/Mode/Directed.hs b/src/Numeric/AD/Mode/Directed.hs
--- a/src/Numeric/AD/Mode/Directed.hs
+++ b/src/Numeric/AD/Mode/Directed.hs
@@ -2,7 +2,7 @@
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.AD.Mode.Directed
--- Copyright   :  (c) Edward Kmett 2010
+-- Copyright   :  (c) Edward Kmett 2010-12
 -- License     :  BSD3
 -- Maintainer  :  ekmett@gmail.com
 -- Stability   :  experimental
@@ -33,6 +33,7 @@
 import qualified Numeric.AD.Mode.Reverse as R
 import qualified Numeric.AD.Mode.Forward as F
 import qualified Numeric.AD.Mode.Tower as T
+import qualified Numeric.AD.Mode.Chain as C
 import qualified Numeric.AD as M
 import Data.Ix
 
@@ -41,6 +42,7 @@
 data Direction
     = Forward
     | Reverse
+    | Chain
     | Tower
     | Mixed
     deriving (Show, Eq, Ord, Read, Bounded, Enum, Ix)
@@ -48,6 +50,7 @@
 diff :: Num a => Direction -> (forall s. Mode s => AD s a -> AD s a) -> a -> a
 diff Forward = F.diff
 diff Reverse = R.diff
+diff Chain = C.diff
 diff Tower = T.diff
 diff Mixed = F.diff
 {-# INLINE diff #-}
@@ -55,6 +58,7 @@
 diff' :: Num a => Direction -> (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)
 diff' Forward = F.diff'
 diff' Reverse = R.diff'
+diff' Chain = C.diff'
 diff' Tower = T.diff'
 diff' Mixed = F.diff'
 {-# INLINE diff' #-}
@@ -62,6 +66,7 @@
 jacobian :: (Traversable f, Traversable g, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)
 jacobian Forward = F.jacobian
 jacobian Reverse = R.jacobian
+jacobian Chain = C.jacobian
 jacobian Tower = F.jacobian -- error "jacobian Tower: unimplemented"
 jacobian Mixed = M.jacobian
 {-# INLINE jacobian #-}
@@ -69,6 +74,7 @@
 jacobian' :: (Traversable f, Traversable g, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)
 jacobian' Forward = F.jacobian'
 jacobian' Reverse = R.jacobian'
+jacobian' Chain = C.jacobian'
 jacobian' Tower = F.jacobian' -- error "jacobian' Tower: unimplemented"
 jacobian' Mixed = M.jacobian'
 {-# INLINE jacobian' #-}
@@ -76,6 +82,7 @@
 grad :: (Traversable f, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a
 grad Forward = F.grad
 grad Reverse = R.grad
+grad Chain   = C.grad
 grad Tower   = F.grad -- error "grad Tower: unimplemented"
 grad Mixed   = M.grad
 {-# INLINE grad #-}
@@ -83,7 +90,7 @@
 grad' :: (Traversable f, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)
 grad' Forward = F.grad'
 grad' Reverse = R.grad'
+grad' Chain   = C.grad'
 grad' Tower   = F.grad' -- error "grad' Tower: unimplemented"
 grad' Mixed   = M.grad'
 {-# INLINE grad' #-}
-
diff --git a/src/Numeric/AD/Mode/Forward.hs b/src/Numeric/AD/Mode/Forward.hs
--- a/src/Numeric/AD/Mode/Forward.hs
+++ b/src/Numeric/AD/Mode/Forward.hs
@@ -49,64 +49,90 @@
 import Numeric.AD.Internal.Composition
 import Numeric.AD.Internal.Forward
 
+-- | Compute the directional derivative of a function given a zipped up 'Functor' of the input values and their derivatives
 du :: (Functor f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> a
 du f = tangent . f . fmap (uncurry bundle)
 {-# INLINE du #-}
 
+-- | Compute the answer and directional derivative of a function given a zipped up 'Functor' of the input values and their derivatives
 du' :: (Functor f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> (a, a)
 du' f = unbundle . f . fmap (uncurry bundle)
 {-# INLINE du' #-}
 
+-- | Compute a vector of directional derivatives for a function given a zipped up 'Functor' of the input values and their derivatives.
 duF :: (Functor f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f (a, a) -> g a
 duF f = fmap tangent . f . fmap (uncurry bundle)
 {-# INLINE duF #-}
 
+-- | Compute a vector of answers and directional derivatives for a function given a zipped up 'Functor' of the input values and their derivatives.
 duF' :: (Functor f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f (a, a) -> g (a, a)
 duF' f = fmap unbundle . f . fmap (uncurry bundle)
 {-# INLINE duF' #-}
 
 -- | The 'diff' function calculates the first derivative of a scalar-to-scalar function by forward-mode 'AD'
 --
--- > diff sin == cos
+-- >>> diff sin 0
+-- 1.0
 diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a
 diff f a = tangent $ apply f a
 {-# INLINE diff #-}
 
--- | The 'd'' function calculates the result and first derivative of scalar-to-scalar function by F'orward' 'AD'
--- 
--- > d' sin == sin &&& cos
--- > d' f = f &&& d f
+-- | The 'diff'' function calculates the result and first derivative of scalar-to-scalar function by 'Forward' mode 'AD'
+--
+-- @
+-- 'diff'' 'sin' == 'sin' 'Control.Arrow.&&&' 'cos'
+-- 'diff'' f = f 'Control.Arrow.&&&' d f
+-- @
+--
+-- >>> diff' sin 0
+-- (0.0,1.0)
+--
+-- >>> diff' exp 0
+-- (1.0,1.0)
+
 diff' :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)
 diff' f a = unbundle $ apply f a
 {-# INLINE diff' #-}
 
--- | The 'diffF' function calculates the first derivative of scalar-to-nonscalar function by F'orward' 'AD'
+-- | The 'diffF' function calculates the first derivatives of scalar-to-nonscalar function by 'Forward' mode 'AD'
+--
+-- >>> diffF (\a -> [sin a, cos a]) 0
+-- [1.0,-0.0]
 diffF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a
 diffF f a = tangent <$> apply f a
 {-# INLINE diffF #-}
 
--- | The 'diffF'' function calculates the result and first derivative of a scalar-to-non-scalar function by F'orward' 'AD'
+-- | The 'diffF'' function calculates the result and first derivatives of a scalar-to-non-scalar function by 'Forward' mode 'AD'
+--
+-- >>> diffF' (\a -> [sin a, cos a]) 0
+-- [(0.0,1.0),(1.0,-0.0)]
 diffF' :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f (a, a)
 diffF' f a = unbundle <$> apply f a
 {-# INLINE diffF' #-}
 
--- | A fast, simple transposed Jacobian computed with forward-mode AD.
+-- | A fast, simple, transposed Jacobian computed with forward-mode AD.
 jacobianT :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> f (g a)
 jacobianT f = bind (fmap tangent . f)
 {-# INLINE jacobianT #-}
 
--- | A fast, simple transposed Jacobian computed with forward-mode AD.
+-- | A fast, simple, transposed Jacobian computed with 'Forward' mode 'AD' that combines the output with the input.
 jacobianWithT :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> f (g b)
 jacobianWithT g f = bindWith g' f
     where g' a ga = g a . tangent <$> ga
 {-# INLINE jacobianWithT #-}
 
+-- | Compute the Jacobian using 'Forward' mode 'AD'. This must transpose the result, so 'jacobianT' is faster and allows more result types.
+--
+--
+-- >>> jacobian (\[x,y] -> [y,x,x+y,x*y,exp x * sin y]) [pi,1]
+-- [[0.0,1.0],[1.0,0.0],[1.0,1.0],[1.0,3.141592653589793],[19.472221418841606,12.502969588876512]]
 jacobian :: (Traversable f, Traversable g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)
 jacobian f as = transposeWith (const id) t p
     where
         (p, t) = bind' (fmap tangent . f) as
 {-# INLINE jacobian #-}
 
+-- | Compute the Jacobian using 'Forward' mode 'AD' and combine the output with the input. This must transpose the result, so 'jacobianWithT' is faster, and allows more result types.
 jacobianWith :: (Traversable f, Traversable g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f b)
 jacobianWith g f as = transposeWith (const id) t p
     where
@@ -114,6 +140,7 @@
         g' a ga = g a . tangent <$> ga
 {-# INLINE jacobianWith #-}
 
+-- | Compute the Jacobian using 'Forward' mode 'AD' along with the actual answer.
 jacobian' :: (Traversable f, Traversable g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)
 jacobian' f as = transposeWith row t p
     where
@@ -121,6 +148,7 @@
         row x as' = (primal x, tangent <$> as')
 {-# INLINE jacobian' #-}
 
+-- | Compute the Jacobian using 'Forward' mode 'AD' combined with the input using a user specified function, along with the actual answer.
 jacobianWith' :: (Traversable f, Traversable g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f b)
 jacobianWith' g f as = transposeWith row t p
     where
@@ -129,29 +157,43 @@
         g' a ga = g a . tangent <$> ga
 {-# INLINE jacobianWith' #-}
 
+-- | Compute the gradient of a function using forward mode AD.
+--
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Chain.grad' for @n@ inputs, in exchange for better space utilization.
 grad :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a
 grad f = bind (tangent . f)
 {-# INLINE grad #-}
 
+-- | Compute the gradient and answer to a function using forward mode AD.
+--
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Chain.grad'' for @n@ inputs, in exchange for better space utilization.
 grad' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)
 grad' f as = (primal b, tangent <$> bs)
     where
         (b, bs) = bind' f as
 {-# INLINE grad' #-}
 
+-- | Compute the gradient of a function using forward mode AD and combine the result with the input using a user-specified function.
+--
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Chain.gradWith' for @n@ inputs, in exchange for better space utilization.
 gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f b
 gradWith g f = bindWith g (tangent . f)
 {-# INLINE gradWith #-}
 
+-- | Compute the gradient of a function using forward mode AD and the answer, and combine the result with the input using a
+-- user-specified function.
+--
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Chain.gradWith'' for @n@ inputs, in exchange for better space utilization.
 gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f b)
 gradWith' g f = bindWith' g (tangent . f)
 {-# INLINE gradWith' #-}
 
--- | Compute the product of a vector with the Hessian using forward-on-forward-mode AD. 
+-- | Compute the product of a vector with the Hessian using forward-on-forward-mode AD.
+--
 hessianProduct :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> f a
 hessianProduct f = duF $ grad $ decomposeMode . f . fmap composeMode
 
--- | Compute the gradient and hessian product using forward-on-forward-mode AD. 
+-- | Compute the gradient and hessian product using forward-on-forward-mode AD.
 hessianProduct' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> f (a, a)
 hessianProduct' f = duF' $ grad $ decomposeMode . f . fmap composeMode
 
diff --git a/src/Numeric/AD/Mode/Reverse.hs b/src/Numeric/AD/Mode/Reverse.hs
--- a/src/Numeric/AD/Mode/Reverse.hs
+++ b/src/Numeric/AD/Mode/Reverse.hs
@@ -50,14 +50,22 @@
 import Numeric.AD.Internal.Classes
 import Numeric.AD.Internal.Composition
 import Numeric.AD.Internal.Reverse
+import Numeric.AD.Internal.Var
 
 -- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with 'Reverse' AD in a single pass.
+--
+-- >>> grad (\[x,y,z] -> x*y+z) [1,2,3]
+-- [2,1,1]
+
 grad :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a
 grad f as = unbind vs (partialArray bds $ f vs)
     where (vs,bds) = bind as
 {-# INLINE grad #-}
 
 -- | The 'grad'' function calculates the result and gradient of a non-scalar-to-scalar function with 'Reverse' AD in a single pass.
+--
+-- >>> grad' (\[x,y,z] -> 4*x*exp y+cos z) [1,2,3]
+-- (28.566231899122155,[29.5562243957226,29.5562243957226,-0.1411200080598672])
 grad' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)
 grad' f as = (primal r, unbind vs $ partialArray bds r)
     where (vs, bds) = bind as
@@ -67,8 +75,12 @@
 -- | @'grad' g f@ function calculates the gradient of a non-scalar-to-scalar function @f@ with reverse-mode AD in a single pass.
 -- The gradient is combined element-wise with the argument using the function @g@.
 --
--- > grad == gradWith (\_ dx -> dx)
--- > id == gradWith const
+-- @
+-- 'grad' = 'gradWith' (\_ dx -> dx)
+-- 'id' = 'gradWith' const
+-- @
+--
+--
 gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f b
 gradWith g f as = unbindWith g vs (partialArray bds $ f vs)
     where (vs,bds) = bind as
@@ -77,7 +89,7 @@
 -- | @'grad'' g f@ calculates the result and gradient of a non-scalar-to-scalar function @f@ with 'Reverse' AD in a single pass
 -- the gradient is combined element-wise with the argument using the function @g@.
 --
--- > grad' == gradWith' (\_ dx -> dx)
+-- @'grad'' == 'gradWith'' (\_ dx -> dx)@
 gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f b)
 gradWith' g f as = (primal r, unbindWith g vs $ partialArray bds r)
     where (vs, bds) = bind as
@@ -85,6 +97,12 @@
 {-# INLINE gradWith' #-}
 
 -- | The 'jacobian' function calculates the jacobian of a non-scalar-to-non-scalar function with reverse AD lazily in @m@ passes for @m@ outputs.
+--
+-- >>> jacobian (\[x,y] -> [y,x,x*y]) [2,1]
+-- [[0,1],[1,0],[1,2]]
+--
+-- >>> jacobian (\[x,y] -> [exp y,cos x,x+y]) [1,2]
+-- [[0.0,7.38905609893065],[-0.8414709848078965,0.0],[1.0,1.0]]
 jacobian :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)
 jacobian f as = unbind vs . partialArray bds <$> f vs where
     (vs, bds) = bind as
@@ -93,6 +111,9 @@
 -- | The 'jacobian'' function calculates both the result and the Jacobian of a nonscalar-to-nonscalar function, using @m@ invocations of reverse AD,
 -- where @m@ is the output dimensionality. Applying @fmap snd@ to the result will recover the result of 'jacobian'
 -- | An alias for 'gradF''
+--
+-- ghci> jacobian' (\[x,y] -> [y,x,x*y]) [2,1]
+-- [(1,[0,1]),(2,[1,0]),(2,[1,2])]
 jacobian' :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)
 jacobian' f as = row <$> f vs where
     (vs, bds) = bind as
@@ -103,9 +124,10 @@
 --
 -- Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@.
 --
--- > jacobian == jacobianWith (\_ dx -> dx)
--- > jacobianWith const == (\f x -> const x <$> f x)
---
+-- @
+-- 'jacobian' = 'jacobianWith' (\_ dx -> dx)
+-- 'jacobianWith' 'const' = (\f x -> 'const' x '<$>' f x)
+-- @
 jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f b)
 jacobianWith g f as = unbindWith g vs . partialArray bds <$> f vs where
     (vs, bds) = bind as
@@ -116,41 +138,66 @@
 --
 -- Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@.
 --
--- > jacobian' == jacobianWith' (\_ dx -> dx)
---
+-- @'jacobian'' == 'jacobianWith'' (\_ dx -> dx)@
 jacobianWith' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f b)
 jacobianWith' g f as = row <$> f vs where
     (vs, bds) = bind as
     row a = (primal a, unbindWith g vs (partialArray bds a))
 {-# INLINE jacobianWith' #-}
 
+-- | Compute the derivative of a function.
+--
+-- >>> diff sin 0
+-- 1.0
+--
+-- >>> cos 0
+-- 1.0
 diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a
 diff f a = derivative $ f (var a 0)
 {-# INLINE diff #-}
 
--- | The 'd'' function calculates the value and derivative, as a
+-- | The 'diff'' function calculates the value and derivative, as a
 -- pair, of a scalar-to-scalar function.
+--
+--
+-- >>> diff' sin 0
+-- (0.0,1.0)
 diff' :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)
 diff' f a = derivative' $ f (var a 0)
 {-# INLINE diff' #-}
 
+-- | Compute the derivatives of a function that returns a vector with regards to its single input.
+--
+-- >>> diffF (\a -> [sin a, cos a]) 0
+-- [1.0,0.0]
 diffF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a
 diffF f a = derivative <$> f (var a 0)
 {-# INLINE diffF #-}
 
+-- | Compute the derivatives of a function that returns a vector with regards to its single input
+-- as well as the primal answer.
+--
+-- >>> diffF' (\a -> [sin a, cos a]) 0
+-- [(0.0,1.0),(1.0,0.0)]
 diffF' :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f (a, a)
 diffF' f a = derivative' <$> f (var a 0)
 {-# INLINE diffF' #-}
 
--- | Compute the hessian via the jacobian of the gradient. gradient is computed in reverse mode and then the jacobian is computed in reverse mode.
+-- | Compute the 'hessian' via the 'jacobian' of the gradient. gradient is computed in reverse mode and then the 'jacobian' is computed in reverse mode.
 --
--- However, since the @'grad f :: f a -> f a'@ is square this is not as fast as using the forward-mode Jacobian of a reverse mode gradient provided by 'Numeric.AD.hessian'.
+-- However, since the @'grad' f :: f a -> f a@ is square this is not as fast as using the forward-mode 'jacobian' of a reverse mode gradient provided by 'Numeric.AD.hessian'.
+--
+-- >>> hessian (\[x,y] -> x*y) [1,2]
+-- [[0,1],[1,0]]
 hessian :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f (f a)
 hessian f = jacobian (grad (decomposeMode . f . fmap composeMode))
 
 -- | Compute the order 3 Hessian tensor on a non-scalar-to-non-scalar function via the reverse-mode Jacobian of the reverse-mode Jacobian of the function.
 --
 -- Less efficient than 'Numeric.AD.Mode.Mixed.hessianF'.
+--
+-- >>> hessianF (\[x,y] -> [x*y,x+y,exp x*cos y]) [1,2]
+-- [[[0.0,1.0],[1.0,0.0]],[[0.0,0.0],[0.0,0.0]],[[-1.1312043837568135,-2.4717266720048188],[-2.4717266720048188,1.1312043837568135]]]
 hessianF :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f (f a))
 hessianF f = decomposeFunctor . jacobian (ComposeFunctor . jacobian (fmap decomposeMode . f . fmap composeMode))
 
diff --git a/src/Numeric/AD/Newton.hs b/src/Numeric/AD/Newton.hs
--- a/src/Numeric/AD/Newton.hs
+++ b/src/Numeric/AD/Newton.hs
@@ -36,11 +36,12 @@
 --
 -- Examples:
 --
---  > take 10 $ findZero (\\x->x^2-4) 1  -- converge to 2.0
---
---  > module Data.Complex
---  > take 10 $ findZero ((+1).(^2)) (1 :+ 1)  -- converge to (0 :+ 1)@
+-- >>> take 10 $ findZero (\x->x^2-4) 1
+-- [1.0,2.5,2.05,2.000609756097561,2.0000000929222947,2.000000000000002,2.0]
 --
+-- >>> import Data.Complex
+-- >>> last $ take 10 $ findZero ((+1).(^2)) (1 :+ 1)
+-- 0.0 :+ 1.0
 findZero :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
 findZero f = go
     where
@@ -49,14 +50,14 @@
                 (y,y') = diff' f x
 {-# INLINE findZero #-}
 
--- | The 'inverseNewton' function inverts a scalar function using
+-- | The 'inverse' function inverts a scalar function using
 -- Newton's method; its output is a stream of increasingly accurate
 -- results.  (Modulo the usual caveats.)
 --
 -- Example:
 --
--- > take 10 $ inverseNewton sqrt 1 (sqrt 10)  -- converges to 10
---
+-- >>> last $ take 10 $ inverse sqrt 1 (sqrt 10)
+-- 10.0
 inverse :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> a -> [a]
 inverse f x0 y = findZero (\x -> f x - lift y) x0
 {-# INLINE inverse  #-}
@@ -65,7 +66,8 @@
 -- function using Newton's method; its output is a stream of
 -- increasingly accurate results.  (Modulo the usual caveats.)
 --
--- > take 10 $ fixedPoint cos 1 -- converges to 0.7390851332151607
+-- >>> last $ take 10 $ fixedPoint cos 1
+-- 0.7390851332151607
 fixedPoint :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
 fixedPoint f = findZero (\x -> f x - x)
 {-# INLINE fixedPoint #-}
@@ -74,7 +76,8 @@
 -- function using Newton's method; produces a stream of increasingly
 -- accurate results.  (Modulo the usual caveats.)
 --
--- > take 10 $ extremum cos 1 -- convert to 0
+-- >>> last $ take 10 $ extremum cos 1
+-- 0.0
 extremum :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
 extremum f = findZero (diff (decomposeMode . f . composeMode))
 {-# INLINE extremum #-}
diff --git a/src/Numeric/AD/Types.hs b/src/Numeric/AD/Types.hs
--- a/src/Numeric/AD/Types.hs
+++ b/src/Numeric/AD/Types.hs
@@ -2,14 +2,13 @@
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.AD.Types
--- Copyright   :  (c) Edward Kmett 2010
+-- Copyright   :  (c) Edward Kmett 2010-12
 -- License     :  BSD3
 -- Maintainer  :  ekmett@gmail.com
 -- Stability   :  experimental
 -- Portability :  GHC only
 --
 -----------------------------------------------------------------------------
-
 module Numeric.AD.Types
     (
     -- * AD modes
@@ -30,20 +29,22 @@
 import Numeric.AD.Internal.Jet
 import Numeric.AD.Internal.Classes
 
--- these exploit the 'magic' that is probed to avoid the need for Functor, etc.
-
+-- | Evaluate a scalar-to-scalar function in the trivial identity AD mode.
 lowerUU :: (forall s. Mode s => AD s a -> AD s a) -> a -> a
 lowerUU f = unprobe . f . probe
 {-# INLINE lowerUU #-}
 
+-- | Evaluate a scalar-to-nonscalar function in the trivial identity AD mode.
 lowerUF :: (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a
 lowerUF f = unprobed . f . probe
 {-# INLINE lowerUF #-}
 
+-- | Evaluate a nonscalar-to-scalar function in the trivial identity AD mode.
 lowerFU :: (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> a
 lowerFU f = unprobe . f . probed
 {-# INLINE lowerFU #-}
 
+-- | Evaluate a nonscalar-to-nonscalar function in the trivial identity AD mode.
 lowerFF :: (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g a
 lowerFF f = unprobed . f . probed
 {-# INLINE lowerFF #-}