diff --git a/Numeric/AD/Internal/Identity.hs b/Numeric/AD/Internal/Identity.hs
--- a/Numeric/AD/Internal/Identity.hs
+++ b/Numeric/AD/Internal/Identity.hs
@@ -1,4 +1,5 @@
 {-# LANGUAGE GeneralizedNewtypeDeriving, FlexibleContexts, FlexibleInstances, MultiParamTypeClasses, DeriveDataTypeable #-}
+{-# OPTIONS_HADDOCK hide #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.AD.Internal.Identity
diff --git a/Numeric/AD/Internal/Reverse.hs b/Numeric/AD/Internal/Reverse.hs
--- a/Numeric/AD/Internal/Reverse.hs
+++ b/Numeric/AD/Internal/Reverse.hs
@@ -32,6 +32,8 @@
     , unbindMap
     , unbindWith
     , unbindMapWithDefault
+    , vgrad, vgrad'
+    , Grad(..)
     ) where
 
 import Prelude hiding (mapM)
@@ -215,3 +217,35 @@
 
 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
+    unpack :: ([a] -> [a]) -> o
+    unpack' :: ([a] -> (a, [a])) -> o'
+
+instance Num a => Grad (AD Reverse a) [a] (a, [a]) a where
+    pack i _ = i
+    unpack f = f []
+    unpack' f = f []
+
+instance Grad i o o' a => Grad (AD Reverse a -> i) (a -> o) (a -> o') a where
+    pack f (a:as) = pack (f a) as
+    pack _ [] = error "Grad.pack: logic error"
+    unpack f a = unpack (f . (a:))
+    unpack' f a = unpack' (f . (a:))
+
+vgrad :: Grad i o o' a => i -> o
+vgrad i = unpack (unsafeGrad (pack i))
+    where
+        unsafeGrad f as = unbind vs (partialArray bds $ f vs)
+            where
+                (vs,bds) = bind as
+
+vgrad' :: Grad i o o' a => i -> o'
+vgrad' i = unpack' (unsafeGrad' (pack i))
+    where
+        unsafeGrad' f as = (primal r, unbind vs (partialArray bds r))
+            where
+                r = f vs
+                (vs,bds) = bind as
+
diff --git a/Numeric/AD/Internal/Sparse.hs b/Numeric/AD/Internal/Sparse.hs
--- a/Numeric/AD/Internal/Sparse.hs
+++ b/Numeric/AD/Internal/Sparse.hs
@@ -1,17 +1,21 @@
-{-# LANGUAGE BangPatterns, TemplateHaskell, TypeFamilies, TypeOperators, FlexibleContexts, UndecidableInstances, DeriveDataTypeable #-}
+{-# LANGUAGE BangPatterns, TemplateHaskell, TypeFamilies, TypeOperators, FlexibleContexts, UndecidableInstances, DeriveDataTypeable, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances #-}
 module Numeric.AD.Internal.Sparse 
     ( Index(..)
     , emptyIndex
     , addToIndex
     , indices
     , Sparse(..)
+    , apply
     , vars
-    , d
-    , d'
-    , ds
+    , d, d', ds
     , skeleton
     , spartial
     , partial
+    , vgrad
+    , vgrad'
+    , vgrads
+    , Grad(..)
+    , Grads(..)
     ) where
 
 import Prelude hiding (lookup)
@@ -65,6 +69,10 @@
         var !n a = (n + 1, AD $ Sparse a $ singleton n $ lift 1)
 {-# INLINE vars #-}
 
+apply :: (Traversable f, Num a) => (f (AD Sparse a) -> b) -> f a -> b
+apply f = f . vars 
+{-# INLINE apply #-}
+
 skeleton :: Traversable f => f a -> f Int
 skeleton = snd . mapAccumL (\ !n _ -> (n + 1, n)) 0
 {-# INLINE skeleton #-}
@@ -86,6 +94,32 @@
             where ix' = addToIndex i ix
 {-# INLINE ds #-}
 
+{-
+vvars :: Num a => Vector a -> Vector (AD Sparse a)
+vvars = Vector.imap (\n a -> AD $ Sparse a $ singleton n $ lift 1)
+{-# INLINE vvars #-}
+
+vapply :: Num a => (Vector (AD Sparse a) -> b) -> Vector a -> b
+vapply f = f . vvars 
+{-# INLINE vapply #-}
+
+
+vd :: Num a => Int -> AD Sparse a -> Vector a
+vd n (AD (Sparse _ da)) = Vector.generate n $ \i -> maybe 0 primal $ lookup i da
+{-# INLINE vd #-}
+
+vd' :: Num a => Int -> AD Sparse a -> (a, Vector a)
+vd' n (AD (Sparse a da)) = (a , Vector.generate n $ \i -> maybe 0 primal $ lookup i da)
+{-# INLINE vd' #-}
+
+vds :: Num a => Int -> AD Sparse a -> Stream Vector a
+vds n (AD as@(Sparse a _)) = a :< Vector.generate n (go emptyIndex)
+    where
+        go ix i = partial (indices ix') as :< Vector.generate n (go ix')
+            where ix' = addToIndex i ix
+{-# INLINE vds #-}
+-}
+
 partial :: Num a => [Int] -> Sparse a -> a
 partial []     (Sparse a _)  = a
 partial (n:ns) (Sparse _ da) = partial ns $ findWithDefault (lift 0) n da
@@ -98,6 +132,8 @@
     spartial ns a'
 {-# INLINE spartial #-}
 
+
+
 instance Primal Sparse where
     primal (Sparse a _) = a
 
@@ -134,3 +170,52 @@
             (mapWithKey (times dadc) dc)
 
 deriveLifted id $ conT ''Sparse
+
+
+class Num a => Grad i o o' a | i -> a o o', o -> a i o', o' -> a i o where
+    pack :: i -> [AD Sparse a] -> AD Sparse a
+    unpack :: ([a] -> [a]) -> o
+    unpack' :: ([a] -> (a, [a])) -> o'
+
+instance Num a => Grad (AD Sparse a) [a] (a, [a]) a where
+    pack i _ = i
+    unpack f = f []
+    unpack' f = f []
+
+instance Grad i o o' a => Grad (AD Sparse a -> i) (a -> o) (a -> o') a where
+    pack f (a:as) = pack (f a) as
+    pack _ [] = error "Grad.pack: logic error"
+    unpack f a = unpack (f . (a:))
+    unpack' f a = unpack' (f . (a:))
+
+vgrad :: Grad i o o' a => i -> o
+vgrad i = unpack (unsafeGrad (pack i))
+    where
+        unsafeGrad f as = d as $ apply f as
+{-# INLINE vgrad #-}
+
+vgrad' :: Grad i o o' a => i -> o'
+vgrad' i = unpack' (unsafeGrad' (pack i))
+    where
+        unsafeGrad' f as = d' as $ apply f as
+{-# INLINE vgrad' #-}
+
+class Num a => Grads i o a | i -> a o, o -> a i where
+    packs :: i -> [AD Sparse a] -> AD Sparse a
+    unpacks :: ([a] -> Stream [] a) -> o
+
+instance Num a => Grads (AD Sparse a) (Stream [] a) a where
+    packs i _ = i
+    unpacks f = f []
+
+instance Grads i o a => Grads (AD Sparse a -> i) (a -> o) a where
+    packs f (a:as) = packs (f a) as
+    packs _ [] = error "Grad.pack: logic error"
+    unpacks f a = unpacks (f . (a:))
+
+vgrads :: Grads i o a => i -> o
+vgrads i = unpacks (unsafeGrads (packs i))
+    where
+        unsafeGrads f as = ds as $ apply f as
+{-# INLINE vgrads #-}
+
diff --git a/Numeric/AD/Internal/Stream.hs b/Numeric/AD/Internal/Stream.hs
--- a/Numeric/AD/Internal/Stream.hs
+++ b/Numeric/AD/Internal/Stream.hs
@@ -59,6 +59,7 @@
 tailS (_ :< as) = as
 {-# INLINE tailS #-}
 
+
 unfoldS :: Functor f => (a -> (b, f a)) -> a -> Stream f b
 unfoldS f a = h :< unfoldS f <$> t 
     where
diff --git a/Numeric/AD/Internal/Tensors.hs b/Numeric/AD/Internal/Tensors.hs
--- a/Numeric/AD/Internal/Tensors.hs
+++ b/Numeric/AD/Internal/Tensors.hs
@@ -62,7 +62,6 @@
         distribute :: Functor f => f (Tensors f a) -> Tensors f (f a)
         distribute x = (headT <$> x) :- distribute (tailT <$> x)
 
-
 instance Typeable1 f => Typeable1 (Tensors f) where
     typeOf1 tfa = mkTyConApp tensorsTyCon [typeOf1 (undefined `asArgsType` tfa)]
         where asArgsType :: f a -> t f a -> f a
diff --git a/Numeric/AD/Mode/Forward.hs b/Numeric/AD/Mode/Forward.hs
--- a/Numeric/AD/Mode/Forward.hs
+++ b/Numeric/AD/Mode/Forward.hs
@@ -40,9 +40,6 @@
     , du'
     , duF
     , duF'
-    -- * Monadic Combinators
-    , diffM
-    , diffM'
     -- * Exposed Types
     , UU, UF, FU, FF
     , AD(..)
@@ -51,7 +48,6 @@
 
 import Data.Traversable (Traversable)
 import Control.Applicative
-import Control.Monad (liftM)
 import Numeric.AD.Types
 import Numeric.AD.Internal.Classes
 import Numeric.AD.Internal.Composition
@@ -97,16 +93,6 @@
 diffF' :: (Functor f, Num a) => UF f a -> a -> f (a, a)
 diffF' f a = unbundle <$> apply f a
 {-# INLINE diffF' #-}
-
--- | The 'dUM' function calculates the first derivative of scalar-to-scalar monadic function by F'orward' 'AD'
-diffM :: (Monad m, Num a) => UF m a -> a -> m a
-diffM f a = tangent `liftM` apply f a
-{-# INLINE diffM #-}
-
--- | The 'd'UM' function calculates the result and first derivative of a scalar-to-scalar monadic function by F'orward' 'AD'
-diffM' :: (Monad m, Num a) => UF m a -> a -> m (a, a)
-diffM' f a = unbundle `liftM` apply f a
-{-# INLINE diffM' #-}
 
 -- | A fast, simple transposed Jacobian computed with forward-mode AD.
 jacobianT :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> f (g a)
diff --git a/Numeric/AD/Mode/Mixed.hs b/Numeric/AD/Mode/Mixed.hs
--- a/Numeric/AD/Mode/Mixed.hs
+++ b/Numeric/AD/Mode/Mixed.hs
@@ -1,4 +1,4 @@
-{-# LANGUAGE Rank2Types, TypeFamilies #-}
+{-# LANGUAGE Rank2Types, TypeFamilies, PatternGuards #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.AD.Mode.Mixed
@@ -11,6 +11,31 @@
 -- Mixed-Mode Automatic Differentiation.
 --
 -- Each combinator exported from this module chooses an appropriate AD mode.
+-- The following basic operations are supported, modified as appropriate by the suffixes below:
+-- 
+-- * 'grad' computes the gradient (partial derivatives) of a function at a point
+--
+-- * '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' compute the Hessian matrix (matrix of second partial derivatives) of a function at a point
+-- 
+-- The suffixes have the following meanings:
+-- 
+-- * @\'@ -- 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.
+--
+-- * @0@ means that the resulting derivative list is padded with 0s at the end.
 -----------------------------------------------------------------------------
 
 module Numeric.AD.Mode.Mixed
@@ -21,41 +46,35 @@
     , gradWith
     , gradWith'
 
-    -- * Jacobians (Mixed Mode)
+    -- * Higher Order Gradients (Sparse-on-Reverse)
+    , grads
+
+    -- * Jacobians (Sparse or Reverse)
     , jacobian
     , jacobian'
     , jacobianWith
     , jacobianWith'
 
-    -- * Monadic Gradient/Jacobian (Reverse Mode)
-    , gradM
-    , gradM'
-    , gradWithM
-    , gradWithM'
-
-    -- * Functorial Gradient/Jacobian (Reverse Mode)
-    , gradF
-    , gradF'
-    , gradWithF
-    , gradWithF'
+    -- * Higher Order Jacobian (Sparse-on-Reverse)
+    , jacobians
 
     -- * Transposed Jacobians (Forward Mode)
     , jacobianT
     , jacobianWithT
 
-    -- * Hessian (Forward-On-Reverse)
+    -- * Hessian (Sparse-On-Reverse)
     , hessian
+    , hessian'
 
-    -- * Hessian Tensors (Forward-On-Mixed)
-    , hessianTensor
+    -- * Hessian Tensors (Sparse or Sparse-On-Reverse)
+    , hessianF
+    -- * Hessian Tensors (Sparse)
+    , hessianF'
 
     -- * Hessian Vector Products (Forward-On-Reverse)
     , hessianProduct
     , hessianProduct'
 
-    -- * Higher Order Gradients/Hessians (Sparse Forward)
-    , gradients
-
     -- * Derivatives (Forward Mode)
     , diff
     , diffF
@@ -90,13 +109,16 @@
     , maclaurin
     , maclaurin0
 
-    -- * Monadic Combinators (Forward Mode)
-    , diffM
-    , diffM'
+    -- * Unsafe Variadic Grad
+    , vgrad
+    , vgrad'
+    , vgrads
 
     -- * Exposed Types
     , module Numeric.AD.Types
     , Mode(..)
+    , Grad
+    , Grads
     ) where
 
 import Data.Traversable (Traversable)
@@ -104,88 +126,84 @@
 import Control.Applicative
 
 import Numeric.AD.Types
-import Numeric.AD.Internal.Identity (probed, unprobe)
 import Numeric.AD.Internal.Composition
 import Numeric.AD.Classes (Mode(..))
 
-import qualified Numeric.AD.Mode.Forward as Forward
 import Numeric.AD.Mode.Forward 
     ( diff, diff', diffF, diffF'
     , du, du', duF, duF'
-    , diffM, diffM'
-    , jacobianT, jacobianWithT
-    ) 
+    , jacobianT, jacobianWithT ) 
 
 import Numeric.AD.Mode.Tower 
     ( diffsF, diffs0F, diffs, diffs0
     , taylor, taylor0, maclaurin, maclaurin0
-    , dus, dus0, dusF, dus0F
-    )
+    , dus, dus0, dusF, dus0F )
 
 import qualified Numeric.AD.Mode.Reverse as Reverse
 import Numeric.AD.Mode.Reverse 
-    ( grad, grad', gradWith, gradWith'
-    , gradM, gradM', gradWithM, gradWithM'
-    , gradF, gradF', gradWithF, gradWithF'
-    )
+    ( grad, grad', gradWith, gradWith', vgrad, vgrad', Grad)
 
 -- temporary until we make a full sparse mode
-import qualified Numeric.AD.Internal.Sparse as Sparse
-
+import qualified Numeric.AD.Mode.Sparse as Sparse
+import Numeric.AD.Mode.Sparse
+    ( grads, jacobians, hessian', hessianF', vgrads, Grads)
+    
 -- | Calculate the Jacobian of a non-scalar-to-non-scalar function, automatically choosing between forward and reverse mode AD based on the number of inputs and outputs.
 --
--- If you need to support functions where the output is only a 'Functor' or 'Monad', consider 'Numeric.AD.Reverse.jacobian' or 'Numeric.AD.Reverse.gradM' from "Numeric.AD.Reverse".
-jacobian :: (Traversable f, Traversable g, Num a) => FF f g a -> f a -> g (f a)
+-- If you know the relative number of inputs and outputs, consider 'Numeric.AD.Reverse.jacobian' or 'Nuneric.AD.Sparse.jacobian'.
+jacobian :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (f a)
 jacobian f bs = snd <$> jacobian' f bs
 {-# INLINE jacobian #-}
 
--- | Calculate both the answer and Jacobian of a non-scalar-to-non-scalar function, automatically choosing between forward- and reverse- mode AD based on the relative, number of inputs and outputs.
+data Nat = Z | S Nat deriving (Eq, Ord)
+
+size :: Foldable f => f a -> Nat
+size = foldr' (\_ b -> S b) Z 
+
+big :: Nat -> Bool
+big (S (S (S (S (S (S (S (S (S (S _)))))))))) = True
+big _ = False
+
+-- | Calculate both the answer and Jacobian of a non-scalar-to-non-scalar function, automatically choosing between forward- and reverse- mode AD based on the relative, based on the number of inputs
 --
--- If you need to support functions where the output is only a 'Functor' or 'Monad', consider 'Numeric.AD.Reverse.jacobian'' or 'Numeric.AD.Reverse.gradM'' from "Numeric.AD.Reverse".
-jacobian' :: (Traversable f, Traversable g, Num a) => FF f g a -> f a -> g (a, f a)
-jacobian' f bs | n == 0    = fmap (\x -> (unprobe x, bs)) as
-               | n > m     = Reverse.jacobian' f bs
-               | otherwise = Forward.jacobian' f bs
+-- If you know the relative number of inputs and outputs, consider 'Numeric.AD.Reverse.jacobian'' or 'Nuneric.AD.Sparse.jacobian''.
+jacobian' :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (a, f a)
+jacobian' f bs | Z <- n = fmap (\x -> (unprobe x, bs)) (f (probed bs))
+               | big n  = Reverse.jacobian' f bs
+               | otherwise = Sparse.jacobian' f bs
     where
-        as = f (probed bs)
         n = size bs
-        m = size as
-        size :: Foldable f => f a -> Int
-        size = foldr' (\_ b -> 1 + b) 0
 {-# INLINE jacobian' #-}
 
 -- | @'jacobianWith' g f@ calculates the Jacobian of a non-scalar-to-non-scalar function, automatically choosing between forward and reverse mode AD based on the number of inputs and outputs.
 --
 -- The resulting Jacobian matrix is then recombined element-wise with the input using @g@.
 --
--- If you need to support functions where the output is only a 'Functor' or 'Monad', consider 'Numeric.AD.Reverse.jacobianWith' or 'Numeric.AD.Reverse.gradWithM' from "Numeric.AD.Reverse".
-jacobianWith :: (Traversable f, Traversable g, Num a) => (a -> a -> b) -> FF f g a -> f a -> g (f b)
+-- If you know the relative number of inputs and outputs, consider 'Numeric.AD.Reverse.jacobianWith' or 'Nuneric.AD.Sparse.jacobianWith'.
+jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> FF f g a -> f a -> g (f b)
 jacobianWith g f bs = snd <$> jacobianWith' g f bs
 {-# INLINE jacobianWith #-}
 
--- | @'jacobianWith'' g f@ calculates the answer and Jacobian of a non-scalar-to-non-scalar function, automatically choosing between forward and reverse mode AD based on the number of inputs and outputs.
+-- | @'jacobianWith'' g f@ calculates the answer and Jacobian of a non-scalar-to-non-scalar function, automatically choosing between sparse and reverse mode AD based on the number of inputs and outputs.
 --
 -- The resulting Jacobian matrix is then recombined element-wise with the input using @g@.
 --
--- If you need to support functions where the output is only a 'Functor' or 'Monad', consider 'Numeric.AD.Reverse.jacobianWith'' or 'Numeric.AD.Reverse.gradWithM'' from "Numeric.AD.Reverse".
-jacobianWith' :: (Traversable f, Traversable g, Num a) => (a -> a -> b) -> FF f g a -> f a -> g (a, f b)
+-- If you know the relative number of inputs and outputs, consider 'Numeric.AD.Reverse.jacobianWith'' or 'Nuneric.AD.Sparse.jacobianWith''.
+jacobianWith' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> FF f g a -> f a -> g (a, f b)
 jacobianWith' g f bs
-    | n == 0    = fmap (\x -> (unprobe x, undefined <$> bs)) as
-    | n > m     = Reverse.jacobianWith' g f bs
-    | otherwise = Forward.jacobianWith' g f bs
+    | Z <- n = fmap (\x -> (unprobe x, undefined <$> bs)) (f (probed bs))
+    | big n  = Reverse.jacobianWith' g f bs
+    | otherwise = Sparse.jacobianWith' g f bs
     where
-        as = f (probed bs)
         n = size bs
-        m = size as
-        size :: Foldable f => f a -> Int
-        size = foldr' (\_ b -> 1 + b) 0
 {-# INLINE jacobianWith' #-}
 
 -- | @'hessianProduct' f wv@ computes the product of the hessian @H@ of a non-scalar-to-scalar function @f@ at @w = 'fst' <$> wv@ with a vector @v = snd <$> wv@ using \"Pearlmutter\'s method\" from <http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.29.6143>, which states:
 --
 -- > H v = (d/dr) grad_w (w + r v) | r = 0
 -- 
--- Or in other words, we take the directional derivative of the gradient.
+-- Or in other words, we take the directional derivative of the gradient. The gradient is calculated in reverse mode, then the directional derivative is calculated in forward mode.
+--
 hessianProduct :: (Traversable f, Num a) => FU f a -> f (a, a) -> f a
 hessianProduct f = duF (grad (decomposeMode . f . fmap composeMode))
 
@@ -193,32 +211,16 @@
 --
 -- > H v = (d/dr) grad_w (w + r v) | r = 0
 -- 
--- Or in other words, we take the directional derivative of the gradient.
--- 
+-- Or in other words, we return the gradient and the directional derivative of the gradient. The gradient is calculated in reverse mode, then the directional derivative is calculated in forward mode.
 hessianProduct' :: (Traversable f, Num a) => FU f a -> f (a, a) -> f (a, a)
 hessianProduct' f = duF' (grad (decomposeMode . f . fmap composeMode))
 
--- hessianProductWith' :: (Traversable f, Num a) => (a -> a -> a -> a -> b) -> (forall s. Mode s. f (AD s a) -> AD s a) -> f (a, a) -> f b
-
--- | Compute the hessian via the jacobian of the gradient. gradient is computed in reverse mode and then the jacobian is computed in forward mode.
+-- | Compute the Hessian via the Jacobian of the gradient. gradient is computed in reverse mode and then the Jacobian is computed in sparse (forward) mode.
 hessian :: (Traversable f, Num a) => FU f a -> f a -> f (f a)
-hessian f = Forward.jacobian (grad (decomposeMode . f . fmap composeMode))
-
--- | Compute the order 3 Hessian tensor on a non-scalar-to-non-scalar function via the forward-mode Jacobian of the mixed-mode Jacobian of the function.
-hessianTensor :: (Traversable f, Traversable g, Num a) => FF f g a -> f a -> g (f (f a))
-hessianTensor f = decomposeFunctor . Forward.jacobian (ComposeFunctor . jacobian (fmap decomposeMode . f . fmap composeMode))
-
--- data f :> a = a :< f (f :> a)
--- data f :- a = a :- (f :- f a) | Zero
-{-
-flatten :: (f :> a) -> (f :- a)
-grads :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (f :- a) 
-grads f b = a :- da :- d2a :- Zero
-    (a, da) = grad2 f a
-    dda = Forward.jacobian (grad (decomposeMode . f . fmap composeMode)
-    ddda = Forward
--}
+hessian f = Sparse.jacobian (grad (decomposeMode . f . fmap composeMode))
 
-gradients :: (Traversable f, Num a) => FU f a -> f a -> Stream f a
-gradients f as = Sparse.ds as $ f $ Sparse.vars as
-    
+-- | 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) => FF f g 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
+    | otherwise = Sparse.hessianF f as
diff --git a/Numeric/AD/Mode/Reverse.hs b/Numeric/AD/Mode/Reverse.hs
--- a/Numeric/AD/Mode/Reverse.hs
+++ b/Numeric/AD/Mode/Reverse.hs
@@ -1,4 +1,5 @@
-{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns #-}
+-- {-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns #-}
+{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, UndecidableInstances, ScopedTypeVariables #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.AD.Mode.Reverse
@@ -24,6 +25,7 @@
     , grad'
     , gradWith
     , gradWith'
+
     -- * Jacobian
     , jacobian
     , jacobian'
@@ -31,34 +33,22 @@
     , jacobianWith'
     -- * Hessian
     , hessian
-    , hessianM
-    , hessianTensor
-    
+    , hessianF
     -- * Derivatives
     , diff
     , diff'
     , diffF
     , diffF'
-    -- * Monadic Combinators
-    , diffM
-    , diffM'
-    , gradM
-    , gradM'
-    , gradWithM
-    , gradWithM'
-    -- * Synonyms
-    , gradF
-    , gradF'
-    , gradWithF
-    , gradWithF'
+    -- * Unsafe Variadic Gradient
+    , vgrad, vgrad'
     -- * Exposed Types
     , UU, UF, FU, FF
     , AD(..)
     , Mode(..)
+    , Grad
     ) where
 
-import Control.Monad (liftM)
-import Control.Applicative (WrappedMonad(..),(<$>))
+import Control.Applicative ((<$>))
 import Data.Traversable (Traversable)
 
 import Numeric.AD.Types
@@ -99,23 +89,14 @@
           r = f vs
 {-# INLINE gradWith' #-}
 
--- | The 'gradF' function calculates the jacobian of a non-scalar-to-non-scalar function with reverse AD lazily in @m@ passes for @m@ outputs.
-gradF :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (f a)
-gradF = jacobian
-{-# INLINE gradF #-}
-
--- | An alias for 'gradF'
+-- | 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 :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (f a)
 jacobian f as = unbind vs . partialArray bds <$> f vs where
     (vs, bds) = bind as
 {-# INLINE jacobian #-}
 
--- | The 'gradF'' 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 'gradF'
-gradF' :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (a, f a)
-gradF' = jacobian' 
-{-# INLINE gradF' #-}
-
+-- | 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' :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (a, f a)
 jacobian' f as = row <$> f vs where
@@ -123,39 +104,29 @@
     row a = (primal a, unbind vs (partialArray bds a))
 {-# INLINE jacobian' #-}
 
--- | 'gradWithF g f' calculates the Jacobian of a non-scalar-to-non-scalar function @f@ with reverse AD lazily in @m@ passes for @m@ outputs.
+-- | '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@.
 --
--- > gradF == gradWithF (\_ dx -> dx)
--- > gradWithF const == (\f x -> const x <$> f x)
+-- > jacobian == jacobianWith (\_ dx -> dx)
+-- > jacobianWith const == (\f x -> const x <$> f x)
 --
-gradWithF :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> FF f g a -> f a -> g (f b)
-gradWithF g f as = unbindWith g vs . partialArray bds <$> f vs where
-    (vs, bds) = bind as
-{-# INLINE gradWithF #-}
-
--- | An alias for 'gradWithF'.
 jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> FF f g a -> f a -> g (f b)
-jacobianWith = gradWithF 
+jacobianWith g f as = unbindWith g vs . partialArray bds <$> f vs where
+    (vs, bds) = bind as
 {-# INLINE jacobianWith #-}
 
--- | 'gradWithF' 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 'gradWithF'
+-- | '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' == gradWithF' (\_ dx -> dx)
+-- > jacobian' == jacobianWith' (\_ dx -> dx)
 --
-gradWithF' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> FF f g a -> f a -> g (a, f b)
-gradWithF' g f as = row <$> f vs where
+jacobianWith' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> FF f g 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 gradWithF' #-}
-
--- | An alias for 'gradWithF''
-jacobianWith' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> FF f g a -> f a -> g (a, f b)
-jacobianWith' = gradWithF'
 {-# INLINE jacobianWith' #-}
 
 diff :: Num a => UU a -> a -> a
@@ -176,44 +147,15 @@
 diffF' f a = derivative' <$> f (var a 0)
 {-# INLINE diffF' #-}
 
--- * Monadic Combinators
-
-diffM :: (Monad m, Num a) => UF m a -> a -> m a
-diffM f a = liftM derivative $ f (var a 0)
-{-# INLINE diffM #-}
-
-diffM' :: (Monad m, Num a) => UF m a -> a -> m (a, a)
-diffM' f a = liftM derivative' $ f (var a 0)
-{-# INLINE diffM' #-}
-
-gradM :: (Traversable f, Monad m, Num a) => FF f m a -> f a -> m (f a)
-gradM f = unwrapMonad . jacobian (WrapMonad . f)
-{-# INLINE gradM #-}
-
-gradM' :: (Traversable f, Monad m, Num a) => FF f m a -> f a -> m (a, f a)
-gradM' f = unwrapMonad . jacobian' (WrapMonad . f)
-{-# INLINE gradM' #-}
-
-gradWithM :: (Traversable f, Monad m, Num a) => (a -> a -> b) -> FF f m a -> f a -> m (f b)
-gradWithM g f = unwrapMonad . jacobianWith g (WrapMonad . f)
-
-gradWithM' :: (Traversable f, Monad m, Num a) => (a -> a -> b) -> FF f m a -> f a -> m (a, f b)
-gradWithM' g f = unwrapMonad . jacobianWith' g (WrapMonad . f)
-
 -- | 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' in "Numeric.AD".
+-- 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 :: (Traversable f, Num a) => FU f 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 forward-mode Jacobian of the mixed-mode Jacobian of the function.
+-- | 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.
 --
--- While this is less efficient than 'Numeric.AD.hessianTensor' from "Numeric.AD" or 'Numeric.AD.Forward.hessianTensor' from "Numeric.AD.Forward", the type signature is more permissive with regards to the output non-scalar, and it may be more efficient if only a few coefficients of the result are consumed.
-hessianTensor :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (f (f a))
-hessianTensor f = decomposeFunctor . jacobian (ComposeFunctor . jacobian (fmap decomposeMode . f . fmap composeMode))
+-- Less efficient than 'Numeric.AD.Mode.Mixed.hessianF'.
+hessianF :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (f (f a))
+hessianF f = decomposeFunctor . jacobian (ComposeFunctor . jacobian (fmap decomposeMode . f . fmap composeMode))
 
--- | Compute the hessian via the reverse-mode jacobian of the reverse-mode gradient of a non-scalar-to-scalar monadic action. 
---
--- While this is less efficient than 'Numeric.AD.hessianTensor' from "Numeric.AD" or 'Numeric.AD.Forward.hessianTensor' from "Numeric.AD.Forward", the type signature is more permissive with regards to the output non-scalar, and it may be more efficient if only a few coefficients of the result are consumed.
-hessianM :: (Traversable f, Monad m, Num a) => FF f m a -> f a -> m (f (f a))
-hessianM f = unwrapMonad . hessianTensor (WrapMonad . f)
diff --git a/Numeric/AD/Mode/Sparse.hs b/Numeric/AD/Mode/Sparse.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/AD/Mode/Sparse.hs
@@ -0,0 +1,123 @@
+{-# LANGUAGE Rank2Types, BangPatterns #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      : Numeric.AD.Mode.Sparse
+-- Copyright   : (c) Edward Kmett 2010
+-- License     : BSD3
+-- Maintainer  : ekmett@gmail.com
+-- Stability   : experimental
+-- Portability : GHC only
+--
+-- Higher order derivatives via a \"dual number tower\".
+--
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Mode.Sparse
+    (
+    -- * Sparse Gradients
+      grad
+    , grad'
+    , gradWith
+    , gradWith'
+    , grads
+    
+    -- * Sparse Jacobians (synonyms)
+    , jacobian
+    , jacobian'
+    , jacobianWith
+    , jacobianWith'
+    , jacobians
+
+    -- * Sparse Hessians
+    , hessian
+    , hessian'
+
+    , hessianF
+    , hessianF'
+
+    -- * Unsafe gradients
+    , vgrad
+    , vgrads
+
+    -- * Exposed Types
+    , module Numeric.AD.Types
+    , Mode(..)
+    , Grad
+    , Grads
+    ) where
+
+import Control.Applicative ((<$>))
+import Data.Traversable
+import Numeric.AD.Types
+import Numeric.AD.Classes
+import Numeric.AD.Internal.Sparse
+import Numeric.AD.Internal.Combinators
+
+second :: (a -> b) -> (c, a) -> (c, b)
+second g (a,b) = (a, g b)
+{-# INLINE second #-}
+
+grad :: (Traversable f, Num a) => FU f a -> f a -> f a
+grad f as = d as $ apply f as
+{-# INLINE grad #-}
+
+grad' :: (Traversable f, Num a) => FU f a -> f a -> (a, f a)
+grad' f as = d' as $ apply f as
+{-# INLINE grad' #-}
+
+gradWith :: (Traversable f, Num a) => (a -> a -> b) -> FU f a -> f a -> f b
+gradWith g f as = zipWithT g as $ grad f as
+{-# INLINE gradWith #-}
+
+gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> FU f a -> f a -> (a, f b)
+gradWith' g f as = second (zipWithT g as) $ grad' f as
+{-# INLINE gradWith' #-}
+
+jacobian :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (f a)
+jacobian f as = d as <$> apply f as
+{-# INLINE jacobian #-}
+
+jacobian' :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (a, f a)
+jacobian' f as = d' as <$> apply f as
+{-# INLINE jacobian' #-}
+
+jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> FF f g a -> f a -> g (f b)
+jacobianWith g f as = zipWithT g as <$> jacobian f as
+{-# INLINE jacobianWith #-}
+
+jacobianWith' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> FF f g a -> f a -> g (a, f b)
+jacobianWith' g f as = second (zipWithT g as) <$> jacobian' f as
+{-# INLINE jacobianWith' #-}
+
+grads :: (Traversable f, Num a) => FU f a -> f a -> Stream f a
+grads f as = ds as $ apply f as
+{-# INLINE grads #-}
+
+jacobians :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (Stream f a)
+jacobians f as = ds as <$> apply f as
+{-# INLINE jacobians #-}
+
+d2 :: Functor f => Stream f a -> f (f a)
+d2 = headT . tailT . tailT . tensors 
+{-# INLINE d2 #-}
+
+d2' :: Functor f => Stream f a -> (a, f (a, f a))
+d2' (a :< as) = (a, fmap (\(da :< das) -> (da, fmap headS das)) as)
+{-# INLINE d2' #-}
+
+hessian :: (Traversable f, Num a) => FU f a -> f a -> f (f a)
+hessian f as = d2 $ grads f as
+{-# INLINE hessian #-}
+
+hessian' :: (Traversable f, Num a) => FU f a -> f a -> (a, f (a, f a))
+hessian' f as = d2' $ grads f as
+{-# INLINE hessian' #-}
+
+hessianF :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (f (f a))
+hessianF f as = d2 <$> jacobians f as
+{-# INLINE hessianF #-}
+
+hessianF' :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (a, f (a, f a))
+hessianF' f as = d2' <$> jacobians f as
+{-# INLINE hessianF' #-}
+
diff --git a/Numeric/AD/Mode/Tower.hs b/Numeric/AD/Mode/Tower.hs
--- a/Numeric/AD/Mode/Tower.hs
+++ b/Numeric/AD/Mode/Tower.hs
@@ -36,16 +36,12 @@
     , duF'    -- answer and directional derivative of (a -> f a)
     , dusF    -- answer and all directional derivatives of (a -> f a)
     , dus0F   -- answer and all zero padded directional derivatives of (a -> a)
-    -- * Monadic Combinators
-    , diffsM  -- answer and all derivatives of the monadic action (a -> m a)
-    , diffs0M -- answer and all zero padded derivatives of (a -> m a)
     -- * Exposed Types
     , UU, UF, FU, FF
     , Mode(..)
     , AD(..)
     ) where
 
-import Control.Monad (liftM)
 import Control.Applicative ((<$>))
 import Numeric.AD.Types
 import Numeric.AD.Classes
@@ -66,14 +62,6 @@
 diffs0F :: (Functor f, Num a) => UF f a -> a -> f [a]
 diffs0F f a = (zeroPad . getADTower) <$> apply f a
 {-# INLINE diffs0F #-}
-
-diffsM :: (Monad m, Num a) => UF m a -> a -> m [a]
-diffsM f a = getADTower `liftM` apply f a
-{-# INLINE diffsM #-}
-
-diffs0M :: (Monad m, Num a) => UF m a -> a -> m [a]
-diffs0M f a = (zeroPad . getADTower) `liftM` apply f a
-{-# INLINE diffs0M #-}
 
 taylor :: Fractional a => UU a -> a -> a -> [a]
 taylor f x dx = go 1 1 (diffs f x)
diff --git a/Numeric/AD/Newton.hs b/Numeric/AD/Newton.hs
--- a/Numeric/AD/Newton.hs
+++ b/Numeric/AD/Newton.hs
@@ -14,18 +14,12 @@
     (
     -- * Newton's Method (Forward AD)
       findZero
-    , findZeroM
     , inverse
-    , inverseM
     , fixedPoint
-    , fixedPointM
     , extremum
-    , extremumM
     -- * Gradient Ascent/Descent (Reverse AD)
     , gradientDescent
-    , gradientDescentM
     , gradientAscent
-    , gradientAscentM
     -- * Exposed Types
     , UU, UF, FU, FF
     , AD(..)
@@ -33,14 +27,12 @@
     ) where
 
 import Prelude hiding (all)
-import Control.Monad (liftM)
-import Data.MList
 import Data.Foldable (all)
 import Data.Traversable (Traversable)
 import Numeric.AD.Types
 import Numeric.AD.Classes
-import Numeric.AD.Mode.Forward (diff, diff', diffM, diffM')
-import Numeric.AD.Mode.Reverse (gradWith', gradWithM')
+import Numeric.AD.Mode.Forward (diff, diff')
+import Numeric.AD.Mode.Reverse (gradWith')
 import Numeric.AD.Internal.Composition
 
 -- | The 'findZero' function finds a zero of a scalar function using
@@ -62,16 +54,6 @@
                 (y,y') = diff' f x
 {-# INLINE findZero #-}
 
-findZeroM :: (Monad m, Fractional a) => UF m a -> a -> MList m a
-findZeroM f x0 = MList (go x0)
-    where
-        go x = return $ 
-               MCons x $ 
-               MList $ do
-                (y,y') <- diffM' f x
-                go (x - y/y')
-{-# INLINE findZeroM #-}
-
 -- | The 'inverseNewton' function inverts a scalar function using
 -- Newton's method; its output is a stream of increasingly accurate
 -- results.  (Modulo the usual caveats.)
@@ -84,10 +66,6 @@
 inverse f x0 y = findZero (\x -> f x - lift y) x0
 {-# INLINE inverse  #-}
 
-inverseM :: (Monad m, Fractional a) => UF m a -> a -> a -> MList m a
-inverseM f x0 y = findZeroM (\x -> subtract (lift y) `liftM` f x) x0
-{-# INLINE inverseM  #-}
-
 -- | The 'fixedPoint' function find a fixedpoint of a scalar
 -- function using Newton's method; its output is a stream of
 -- increasingly accurate results.  (Modulo the usual caveats.)
@@ -97,10 +75,6 @@
 fixedPoint f = findZero (\x -> f x - x)
 {-# INLINE fixedPoint #-}
 
-fixedPointM :: (Monad m, Fractional a) => UF m a -> a -> MList m a
-fixedPointM f = findZeroM (\x -> subtract x `liftM` f x)
-{-# INLINE fixedPointM #-}
-
 -- | The 'extremum' function finds an extremum of a scalar
 -- function using Newton's method; produces a stream of increasingly
 -- accurate results.  (Modulo the usual caveats.)
@@ -110,10 +84,6 @@
 extremum f = findZero (diff (decomposeMode . f . composeMode))
 {-# INLINE extremum #-}
 
-extremumM :: (Monad m, Fractional a) => UF m a -> a -> MList m a
-extremumM f = findZeroM (diffM (liftM decomposeMode . f . composeMode))
-{-# INLINE extremumM #-}
-
 -- | The 'gradientDescent' function performs a multivariate
 -- optimization, based on the naive-gradient-descent in the file
 -- @stalingrad\/examples\/flow-tests\/pre-saddle-1a.vlad@ from the
@@ -141,31 +111,3 @@
 gradientAscent :: (Traversable f, Fractional a, Ord a) => FU f a -> f a -> [f a]
 gradientAscent f = gradientDescent (negate . f)
 {-# INLINE gradientAscent #-}
-
--- monadic gradient descent
-gradientDescentM :: (Traversable f, Monad m, Fractional a, Ord a) => FF f m a -> f a -> MList m (f a)
-gradientDescentM f x0 = MList $ do
-        (fx0, xgx0) <- gradWithM' (,) f x0
-        go x0 fx0 xgx0 0.1 (0 :: Int)
-    where
-        go x fx xgx !eta !i
-            | eta == 0  = return MNil -- step size is 0
-            | otherwise = do
-                (fx1, xgx1) <- gradWithM' (,) f x1
-                case () of
-                 _ | fx1 > fx     -> go x fx xgx (eta/2) 0 -- we stepped too far
-                   | zeroGrad xgx -> return MNil -- gradient is 0
-                   | otherwise    -> return $ 
-                                     MCons x1 $ 
-                                     MList $
-                                     if i == 10
-                                     then go x1 fx1 xgx1 (eta*2) 0
-                                     else go x1 fx1 xgx1 eta (i+1)
-            where
-                x1 = fmap (\(xi,gxi) -> xi - eta * gxi) xgx
-                zeroGrad = all (\(_,g) -> g == 0)
-{-# INLINE gradientDescentM #-}
-
-gradientAscentM :: (Traversable f, Monad m, Fractional a, Ord a) => FF f m a -> f a -> MList m (f a)
-gradientAscentM f = gradientDescentM (liftM negate . f)
-{-# INLINE gradientAscentM #-}
diff --git a/Numeric/AD/Types.hs b/Numeric/AD/Types.hs
--- a/Numeric/AD/Types.hs
+++ b/Numeric/AD/Types.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE Rank2Types #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.AD.Types
@@ -24,8 +25,33 @@
     , headS
     , tailS
     , unfoldS
+    -- * An Identity Mode. 
+    , Id(..)
+    , probe, unprobe
+    , probed, unprobed
+    -- * Apply functions that use 'lift'
+    , lowerUU, lowerUF, lowerFU, lowerFF
     ) where
 
+import Numeric.AD.Internal.Identity
 import Numeric.AD.Internal.Types
 import Numeric.AD.Internal.Stream
 import Numeric.AD.Internal.Tensors
+
+-- these exploit the 'magic' that is probed to avoid the need for Functor, etc.
+
+lowerUU :: UU a -> a -> a
+lowerUU f = unprobe . f . probe
+{-# INLINE lowerUU #-}
+
+lowerUF :: UF f a -> a -> f a
+lowerUF f = unprobed . f . probe
+{-# INLINE lowerUF #-}
+
+lowerFU :: FU f a -> f a -> a
+lowerFU f = unprobe . f . probed
+{-# INLINE lowerFU #-}
+
+lowerFF :: FF f g a -> f a -> g a
+lowerFF f = unprobed . f . probed
+{-# INLINE lowerFF #-}
diff --git a/ad.cabal b/ad.cabal
--- a/ad.cabal
+++ b/ad.cabal
@@ -1,5 +1,5 @@
 name:         ad
-version:      0.40.1
+version:      0.44.0
 license:      BSD3
 license-File: LICENSE
 copyright:    (c) Edward Kmett 2010,
@@ -16,8 +16,61 @@
     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.
     .
-    The combinators in "Numeric.AD" choose from a variety of automatic differentiation modes,
-    based on the arity of their inputs and outputs.
+    Each mode has a separate module full of combinators.
+    .
+    * @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.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.
+    .
+    * @Numeric.AD.Mode.Mixed@ computes using whichever mode or combination thereof is suitable to each individual combinator. This mode is the default, re-exported by @Numeric.AD@
+    .
+    . 
+    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 (partial derivatives) of a function at a point.
+    .
+    * '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.
+    .
+    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.
+    .
+    * @0@ means that the resulting derivative list is padded with 0s at the end.
+    .
+    Changes since 0.40.0
+    .
+    * Bug fix in the derivative calculated for @'(/)' :: (Mode s, Fractional a) => AD s a@
+    .
+    * Improved documentation
+    .
+    * Regularized naming conventions
+    .
+    * Exposed 'Id', probe, and lower methods via @Numeric.AD.Types@
+    .
+    * Removed monadic combinators
+    .
+    * Retuned the 'Mixed' mode jacobian calculations to only require a 'Functor'-based result.
+    .
+    * Added unsafe variadic 'vgrad', 'vgrad'', and 'vgrads' combinators
 
 build-type:   Simple
 build-depends:       
@@ -25,7 +78,6 @@
     data-reify >= 0.5 && < 0.6, 
     containers >= 0.2 && < 0.4,
     template-haskell >= 2.4 && < 2.5,
-    mlist >= 0.0.2 && <= 0.1,
     array >= 0.2 && < 0.4
 
 exposed-modules:
@@ -43,20 +95,20 @@
     Numeric.AD.Internal.Sparse
     Numeric.AD.Internal.Dense
     Numeric.AD.Internal.Composition
-    Numeric.AD.Internal.Identity
 
     Numeric.AD.Mode.Directed
     Numeric.AD.Mode.Forward
     Numeric.AD.Mode.Mixed
     Numeric.AD.Mode.Reverse
     Numeric.AD.Mode.Tower
-
+    Numeric.AD.Mode.Sparse
 
 other-modules:
     Numeric.AD.Internal.Types
     Numeric.AD.Internal.Comonad
     Numeric.AD.Internal.Stream
     Numeric.AD.Internal.Tensors
+    Numeric.AD.Internal.Identity
 
 Extra-Source-Files: TODO
 GHC-Options: -Wall -fspec-constr -fdicts-cheap -O2