diff --git a/Numeric/AD.hs b/Numeric/AD.hs
--- a/Numeric/AD.hs
+++ b/Numeric/AD.hs
@@ -1,4 +1,3 @@
-{-# LANGUAGE Rank2Types, TypeFamilies #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.AD
@@ -8,192 +7,14 @@
 -- Stability   :  experimental
 -- Portability :  GHC only
 --
--- Mixed-Mode Automatic Differentiation.
+-- Mixed-mode automatic differentiation combinators.
 --
--- Each combinator exported from this module chooses an appropriate AD mode.
 -----------------------------------------------------------------------------
 
 module Numeric.AD
-    (
-    -- * Gradients (Reverse Mode)
-      grad
-    , grad'
-    , gradWith
-    , gradWith'
-
-    -- * Jacobians (Mixed Mode)
-    , jacobian
-    , jacobian'
-    , jacobianWith
-    , jacobianWith'
-
-    -- * Monadic Gradient/Jacobian (Reverse Mode)
-    , gradM
-    , gradM'
-    , gradWithM
-    , gradWithM'
-
-    -- * Functorial Gradient/Jacobian (Reverse Mode)
-    , gradF
-    , gradF'
-    , gradWithF
-    , gradWithF'
-
-    -- * Transposed Jacobians (Forward Mode)
-    , jacobianT
-    , jacobianWithT
-
-    -- * Hessian (Forward-On-Reverse)
-    , hessian
-
-    -- * Hessian Tensors (Forward-On-Mixed)
-    , hessianTensor
-
-    -- * Hessian Vector Products (Forward-On-Reverse)
-    , hessianProduct
-    , hessianProduct'
-
-    -- * Derivatives (Forward Mode)
-    , diff
-    , diffF
-
-    , diff'
-    , diffF'
-
-    -- * Derivatives (Tower)
-    , diffs
-    , diffsF
-
-    , diffs0
-    , diffs0F
-
-    -- * Directional Derivatives (Forward Mode)
-    , du
-    , du'
-    , duF
-    , duF'
-
-    -- * Directional Derivatives (Tower)
-    , dus
-    , dus0
-    , dusF
-    , dus0F
-
-    -- * Taylor Series (Tower)
-    , taylor
-    , taylor0
-
-    -- * Maclaurin Series (Tower)
-    , maclaurin
-    , maclaurin0
-
-    -- * Monadic Combinators (Forward Mode)
-    , diffM
-    , diffM'
-
-    -- * Exposed Types
-    , UU, UF, FU, FF
-    , AD(..)
-    , Mode(..)
+    ( module Numeric.AD.Mode.Mixed
+    , module Numeric.AD.Newton
     ) where
 
-import Data.Traversable (Traversable)
-import Data.Foldable (Foldable, foldr')
-import Control.Applicative
-import Numeric.AD.Internal (AD(..), UU, UF, FU, FF)
-import Numeric.AD.Internal.Identity (probed, unprobe)
-import Numeric.AD.Internal.Classes  (Mode(..))
-import Numeric.AD.Forward  (diff, diff', diffF, diffF', du, du', duF, duF', diffM, diffM', jacobianT, jacobianWithT) 
-import Numeric.AD.Tower    (diffsF, diffs0F , diffs, diffs0, taylor, taylor0, maclaurin, maclaurin0, dus, dus0, dusF, dus0F)
-import Numeric.AD.Reverse  (grad, grad', gradWith, gradWith', gradM, gradM', gradWithM, gradWithM', gradF, gradF', gradWithF, gradWithF')
-import Numeric.AD.Internal.Composition
-
-import qualified Numeric.AD.Forward as Forward
-import qualified Numeric.AD.Reverse as Reverse
-
--- | 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)
-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.
---
--- 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
-    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)
-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.
---
--- 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)
-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
-    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.
-hessianProduct :: (Traversable f, Num a) => FU f a -> f (a, a) -> f a
-hessianProduct f = duF (grad (decomposeMode . f . fmap composeMode))
-
--- | @'hessianProduct'' f wv@ computes both the gradient of a non-scalar-to-scalar @f@ at @w = 'fst' <$> wv@ and the product of the hessian @H@ at @w@ with a vector @v = snd <$> wv@ using \"Pearlmutter's method\". The outputs are returned wrapped in the same functor.
---
--- > H v = (d/dr) grad_w (w + r v) | r = 0
--- 
--- Or in other words, we take the directional derivative of the gradient.
--- 
-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.
-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
--}
+import Numeric.AD.Mode.Mixed
+import Numeric.AD.Newton hiding (Mode(..), AD(..), UU, UF, FU, FF)
diff --git a/Numeric/AD/Classes.hs b/Numeric/AD/Classes.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/AD/Classes.hs
@@ -0,0 +1,22 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD.Classes
+-- Copyright   :  (c) Edward Kmett 2010
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only
+--
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Classes
+    ( 
+    -- * AD Modes
+      Mode(..)
+    -- * Comonads
+    , Copointed(..)
+    , Comonad(..)
+    ) where
+
+import Numeric.AD.Internal.Classes
+import Numeric.AD.Internal.Comonad
diff --git a/Numeric/AD/Directed.hs b/Numeric/AD/Directed.hs
deleted file mode 100644
--- a/Numeric/AD/Directed.hs
+++ /dev/null
@@ -1,92 +0,0 @@
-{-# LANGUAGE Rank2Types #-}
------------------------------------------------------------------------------
--- |
--- Module      :  Numeric.AD.Directed
--- Copyright   :  (c) Edward Kmett 2010
--- License     :  BSD3
--- Maintainer  :  ekmett@gmail.com
--- Stability   :  experimental
--- Portability :  GHC only
---
--- Allows the choice of AD 'Mode' to be specified at the term level for
--- benchmarking or more complicated usage patterns.
------------------------------------------------------------------------------
-
-module Numeric.AD.Directed
-    (
-    -- * Gradients
-      grad
-    , grad'
-    -- * Jacobians
-    , jacobian
-    , jacobian'
-    -- * Derivatives
-    , diff
-    , diff'
-    -- * Exposed Types
-    , UU, UF, FU, FF
-    , Direction(..)
-    , Mode(..)
-    , AD(..)
-    ) where
-
-import Prelude hiding (reverse)
-import Numeric.AD.Internal
-import Data.Traversable (Traversable)
-import qualified Numeric.AD.Reverse as R
-import qualified Numeric.AD.Forward as F
-import qualified Numeric.AD.Tower as T
-import qualified Numeric.AD as M
-import Data.Ix
-
--- TODO: use a data types a la carte approach, so we can expose more methods here
--- rather than just the intersection of all of the functionality
-data Direction
-    = Forward
-    | Reverse
-    | Tower
-    | Mixed
-    deriving (Show, Eq, Ord, Read, Bounded, Enum, Ix)
-
-diff :: Num a => Direction -> UU a -> a -> a
-diff Forward = F.diff
-diff Reverse = R.diff
-diff Tower = T.diff
-diff Mixed = F.diff
-{-# INLINE diff #-}
-
-diff' :: Num a => Direction -> UU a -> a -> (a, a)
-diff' Forward = F.diff'
-diff' Reverse = R.diff'
-diff' Tower = T.diff'
-diff' Mixed = F.diff'
-{-# INLINE diff' #-}
-
-jacobian :: (Traversable f, Traversable g, Num a) => Direction -> FF f g a -> f a -> g (f a)
-jacobian Forward = F.jacobian
-jacobian Reverse = R.jacobian
-jacobian Tower = F.jacobian -- error "jacobian Tower: unimplemented"
-jacobian Mixed = M.jacobian
-{-# INLINE jacobian #-}
-
-jacobian' :: (Traversable f, Traversable g, Num a) => Direction -> FF f g a -> f a -> g (a, f a)
-jacobian' Forward = F.jacobian'
-jacobian' Reverse = R.jacobian'
-jacobian' Tower = F.jacobian' -- error "jacobian' Tower: unimplemented"
-jacobian' Mixed = M.jacobian'
-{-# INLINE jacobian' #-}
-
-grad :: (Traversable f, Num a) => Direction -> FU f a -> f a -> f a
-grad Forward = F.grad
-grad Reverse = R.grad
-grad Tower   = F.grad -- error "grad Tower: unimplemented"
-grad Mixed   = M.grad
-{-# INLINE grad #-}
-
-grad' :: (Traversable f, Num a) => Direction -> FU f a -> f a -> (a, f a)
-grad' Forward = F.grad'
-grad' Reverse = R.grad'
-grad' Tower   = F.grad' -- error "grad' Tower: unimplemented"
-grad' Mixed   = M.grad'
-{-# INLINE grad' #-}
-
diff --git a/Numeric/AD/Forward.hs b/Numeric/AD/Forward.hs
deleted file mode 100644
--- a/Numeric/AD/Forward.hs
+++ /dev/null
@@ -1,178 +0,0 @@
-{-# LANGUAGE Rank2Types #-}
------------------------------------------------------------------------------
--- |
--- Module      :  Numeric.AD.Forward
--- Copyright   :  (c) Edward Kmett 2010
--- License     :  BSD3
--- Maintainer  :  ekmett@gmail.com
--- Stability   :  experimental
--- Portability :  GHC only
---
--- Forward mode automatic differentiation
---
------------------------------------------------------------------------------
-
-module Numeric.AD.Forward
-    (
-    -- * Gradient
-      grad
-    , grad'
-    , gradWith
-    , gradWith'
-    -- * Jacobian
-    , jacobian
-    , jacobian'
-    , jacobianWith
-    , jacobianWith'
-    -- * Transposed Jacobian
-    , jacobianT
-    , jacobianWithT
-    -- * Hessian Product
-    , hessianProduct
-    , hessianProduct'
-    -- * Derivatives
-    , diff
-    , diff'
-    , diffF
-    , diffF'
-    -- * Directional Derivatives
-    , du
-    , du'
-    , duF
-    , duF'
-    -- * Monadic Combinators
-    , diffM
-    , diffM'
-    -- * Exposed Types
-    , UU, UF, FU, FF
-    , AD(..)
-    , Mode(..)
-    ) where
-
-import Data.Traversable (Traversable)
-import Control.Applicative
-import Control.Monad (liftM)
-import Numeric.AD.Internal
-import Numeric.AD.Internal.Composition
-import Numeric.AD.Internal.Forward
-
-du :: (Functor f, Num a) => FU f a -> f (a, a) -> a
-du f = tangent . f . fmap (uncurry bundle)
-{-# INLINE du #-}
-
-du' :: (Functor f, Num a) => FU f a -> f (a, a) -> (a, a)
-du' f = unbundle . f . fmap (uncurry bundle)
-{-# INLINE du' #-}
-
-duF :: (Functor f, Functor g, Num a) => FF f g a -> f (a, a) -> g a
-duF f = fmap tangent . f . fmap (uncurry bundle)
-{-# INLINE duF #-}
-
-duF' :: (Functor f, Functor g, Num a) => FF f g 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 :: Num a => UU a -> a -> a
-diff f a = tangent $ apply f a
-{-# INLINE diff #-}
-
--- | The 'd'UU' 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
-diff' :: Num a => UU 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'
-diffF :: (Functor f, Num a) => UF f 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'
-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)
-jacobianT f = bind (fmap tangent . f)
-{-# INLINE jacobianT #-}
-
--- | A fast, simple transposed Jacobian computed with forward-mode AD.
-jacobianWithT :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> FF f g a -> f a -> f (g b)
-jacobianWithT g f = bindWith g' f
-    where g' a ga = g a . tangent <$> ga
-{-# INLINE jacobianWithT #-}
-
-jacobian :: (Traversable f, Traversable g, Num a) => FF f g a -> f a -> g (f a)
-jacobian f as = transposeWith (const id) t p
-    where
-        (p, t) = bind' (fmap tangent . f) as
-{-# INLINE jacobian #-}
-
-jacobianWith :: (Traversable f, Traversable g, Num a) => (a -> a -> b) -> FF f g a -> f a -> g (f b)
-jacobianWith g f as = transposeWith (const id) t p
-    where
-        (p, t) = bindWith' g' f as
-        g' a ga = g a . tangent <$> ga
-{-# INLINE jacobianWith #-}
-
-jacobian' :: (Traversable f, Traversable g, Num a) => FF f g a -> f a -> g (a, f a)
-jacobian' f as = transposeWith row t p
-    where
-        (p, t) = bind' f as
-        row x as' = (primal x, tangent <$> as')
-{-# INLINE jacobian' #-}
-
-jacobianWith' :: (Traversable f, Traversable g, Num a) => (a -> a -> b) -> FF f g a -> f a -> g (a, f b)
-jacobianWith' g f as = transposeWith row t p
-    where
-        (p, t) = bindWith' g' f as
-        row x as' = (primal x, as')
-        g' a ga = g a . tangent <$> ga
-{-# INLINE jacobianWith' #-}
-
-grad :: (Traversable f, Num a) => FU f a -> f a -> f a
-grad f = bind (tangent . f)
-{-# INLINE grad #-}
-
-grad' :: (Traversable f, Num a) => FU f a -> f a -> (a, f a)
-grad' f as = (primal b, tangent <$> bs)
-    where
-        (b, bs) = bind' f as
-{-# INLINE grad' #-}
-
-gradWith :: (Traversable f, Num a) => (a -> a -> b) -> FU f a -> f a -> f b
-gradWith g f = bindWith g (tangent . f)
-{-# INLINE gradWith #-}
-
-gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> FU f 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. 
-hessianProduct :: (Traversable f, Num a) => FU f 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. 
-hessianProduct' :: (Traversable f, Num a) => FU f a -> f (a, a) -> f (a, a)
-hessianProduct' f = duF' $ grad $ decomposeMode . f . fmap composeMode
-
--- * Experimental
-
--- data f :> a = a :< f (f :> a)
--- gradients :: (Traversable f, Num a) => FU f a -> f a -> (f :> a)
diff --git a/Numeric/AD/Internal.hs b/Numeric/AD/Internal.hs
deleted file mode 100644
--- a/Numeric/AD/Internal.hs
+++ /dev/null
@@ -1,19 +0,0 @@
--- {-# OPTIONS_HADDOCK hide, prune #-}
------------------------------------------------------------------------------
--- |
--- Module      :  Numeric.AD.Internal
--- Copyright   :  (c) Edward Kmett 2010
--- License     :  BSD3
--- Maintainer  :  ekmett@gmail.com
--- Stability   :  experimental
--- Portability :  GHC only
---
------------------------------------------------------------------------------
-module Numeric.AD.Internal
-    ( module Numeric.AD.Internal.Classes
-    , module Numeric.AD.Internal.Types
-    ) where
-
-import Numeric.AD.Internal.Classes
-import Numeric.AD.Internal.Types
-
diff --git a/Numeric/AD/Internal/Classes.hs b/Numeric/AD/Internal/Classes.hs
--- a/Numeric/AD/Internal/Classes.hs
+++ b/Numeric/AD/Internal/Classes.hs
@@ -28,6 +28,7 @@
 import Control.Applicative
 import Data.Char
 import Language.Haskell.TH
+import Numeric.AD.Internal.Combinators (on)
 
 infixl 8 **!
 infixl 7 *!, /!, ^*, *^, ^/
@@ -155,9 +156,6 @@
 square1 :: (Lifted t, Num a) => t a -> t a
 square1 x = x *! x
 {-# INLINE square1 #-}
-
-on :: (a -> a -> c) -> (b -> a) -> b -> b -> c
-on f g a b = f (g a) (g b)
 
 discrete1 :: (Primal t, Num a) => (a -> c) -> t a -> c
 discrete1 f x = f (primal x)
diff --git a/Numeric/AD/Internal/Comonad.hs b/Numeric/AD/Internal/Comonad.hs
--- a/Numeric/AD/Internal/Comonad.hs
+++ b/Numeric/AD/Internal/Comonad.hs
@@ -1,5 +1,5 @@
 {-# LANGUAGE TypeOperators, TemplateHaskell, ScopedTypeVariables #-}
--- {-# OPTIONS_HADDOCK hide #-}
+{-# OPTIONS_HADDOCK hide #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.AD.Internal.Comonad
diff --git a/Numeric/AD/Internal/Composition.hs b/Numeric/AD/Internal/Composition.hs
--- a/Numeric/AD/Internal/Composition.hs
+++ b/Numeric/AD/Internal/Composition.hs
@@ -23,7 +23,8 @@
 import Data.Typeable (Typeable1(..), Typeable(..), TyCon, mkTyCon, mkTyConApp, typeOfDefault, gcast1)
 import Data.Foldable (Foldable(foldMap))
 import Data.Traversable (Traversable(traverse))
-import Numeric.AD.Internal
+import Numeric.AD.Internal.Classes
+import Numeric.AD.Internal.Types
 
 -- | Functor composition, used to nest the use of jacobian and grad
 newtype ComposeFunctor f g a = ComposeFunctor { decomposeFunctor :: f (g a) }
diff --git a/Numeric/AD/Internal/Forward.hs b/Numeric/AD/Internal/Forward.hs
--- a/Numeric/AD/Internal/Forward.hs
+++ b/Numeric/AD/Internal/Forward.hs
@@ -33,7 +33,8 @@
 import Data.Foldable (Foldable, toList)
 import Data.Data
 import Control.Applicative
-import Numeric.AD.Internal
+import Numeric.AD.Internal.Types
+import Numeric.AD.Internal.Classes
 import Numeric.AD.Internal.Identity
 
 data Forward a = Forward a a deriving (Show, Data, Typeable)
diff --git a/Numeric/AD/Internal/Iterated.hs b/Numeric/AD/Internal/Iterated.hs
--- a/Numeric/AD/Internal/Iterated.hs
+++ b/Numeric/AD/Internal/Iterated.hs
@@ -25,7 +25,8 @@
 import Data.Traversable
 import Data.Data (Data(..), mkDataType, DataType, mkConstr, Constr, constrIndex, Fixity(Infix))
 import Data.Typeable (Typeable1(..), Typeable(..), TyCon, mkTyCon, mkTyConApp, typeOfDefault, gcast1)
-import Numeric.AD.Internal
+import Numeric.AD.Internal.Types
+import Numeric.AD.Internal.Classes
 import Numeric.AD.Internal.Comonad
 import Numeric.AD.Internal.Combinators (on)
 -- import qualified Numeric.AD.Internal.Forward
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
@@ -50,7 +50,8 @@
 import Language.Haskell.TH
 import Data.Data (Data)
 import Data.Typeable (Typeable)
-import Numeric.AD.Internal
+import Numeric.AD.Internal.Types
+import Numeric.AD.Internal.Classes
 import Numeric.AD.Internal.Identity
 
 -- | A @Tape@ records the information needed back propagate from the output to each input during 'Reverse' 'Mode' AD.
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
@@ -1,5 +1,5 @@
 {-# LANGUAGE StandaloneDeriving, FlexibleContexts, UndecidableInstances #-}
--- {-# OPTIONS_HADDOCK hide #-}
+{-# OPTIONS_HADDOCK hide #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.AD.Internal.Stream
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
@@ -1,5 +1,5 @@
 {-# LANGUAGE TypeOperators, TemplateHaskell, ScopedTypeVariables #-}
--- {-# OPTIONS_HADDOCK hide #-}
+{-# OPTIONS_HADDOCK hide #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.AD.Internal.Tensors
@@ -22,7 +22,7 @@
 import Data.Foldable
 import Data.Traversable
 import Data.Monoid
-import Data.Typeable (Typeable1(..), Typeable(..), TyCon, mkTyCon, mkTyConApp, typeOfDefault, gcast1)
+import Data.Typeable (Typeable1(..), Typeable(..), TyCon, mkTyCon, mkTyConApp, typeOfDefault)
 import Numeric.AD.Internal.Comonad
 import Numeric.AD.Internal.Stream
 
diff --git a/Numeric/AD/Internal/Tower.hs b/Numeric/AD/Internal/Tower.hs
--- a/Numeric/AD/Internal/Tower.hs
+++ b/Numeric/AD/Internal/Tower.hs
@@ -32,7 +32,8 @@
 import Data.Data (Data)
 import Data.Typeable (Typeable)
 import Language.Haskell.TH
-import Numeric.AD.Internal
+import Numeric.AD.Internal.Types
+import Numeric.AD.Internal.Classes
 
 -- | @Tower@ is an AD 'Mode' that calculates a tangent tower by forward AD, and provides fast 'diffsUU', 'diffsUF'
 newtype Tower a = Tower { getTower :: [a] } deriving (Data, Typeable)
diff --git a/Numeric/AD/Internal/Types.hs b/Numeric/AD/Internal/Types.hs
--- a/Numeric/AD/Internal/Types.hs
+++ b/Numeric/AD/Internal/Types.hs
@@ -1,5 +1,5 @@
 {-# LANGUAGE Rank2Types, GeneralizedNewtypeDeriving, TemplateHaskell, FlexibleContexts, FlexibleInstances, MultiParamTypeClasses, UndecidableInstances #-}
--- {-# OPTIONS_HADDOCK hide, prune #-}
+{-# OPTIONS_HADDOCK hide #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.AD.Internal.Types
diff --git a/Numeric/AD/Mode/Directed.hs b/Numeric/AD/Mode/Directed.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/AD/Mode/Directed.hs
@@ -0,0 +1,93 @@
+{-# LANGUAGE Rank2Types #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD.Mode.Directed
+-- Copyright   :  (c) Edward Kmett 2010
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only
+--
+-- Allows the choice of AD 'Mode' to be specified at the term level for
+-- benchmarking or more complicated usage patterns.
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Mode.Directed
+    (
+    -- * Gradients
+      grad
+    , grad'
+    -- * Jacobians
+    , jacobian
+    , jacobian'
+    -- * Derivatives
+    , diff
+    , diff'
+    -- * Exposed Types
+    , UU, UF, FU, FF
+    , Direction(..)
+    , Mode(..)
+    , AD(..)
+    ) where
+
+import Prelude hiding (reverse)
+import Numeric.AD.Types
+import Numeric.AD.Classes
+import Data.Traversable (Traversable)
+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.Mixed as M
+import Data.Ix
+
+-- TODO: use a data types a la carte approach, so we can expose more methods here
+-- rather than just the intersection of all of the functionality
+data Direction
+    = Forward
+    | Reverse
+    | Tower
+    | Mixed
+    deriving (Show, Eq, Ord, Read, Bounded, Enum, Ix)
+
+diff :: Num a => Direction -> UU a -> a -> a
+diff Forward = F.diff
+diff Reverse = R.diff
+diff Tower = T.diff
+diff Mixed = F.diff
+{-# INLINE diff #-}
+
+diff' :: Num a => Direction -> UU a -> a -> (a, a)
+diff' Forward = F.diff'
+diff' Reverse = R.diff'
+diff' Tower = T.diff'
+diff' Mixed = F.diff'
+{-# INLINE diff' #-}
+
+jacobian :: (Traversable f, Traversable g, Num a) => Direction -> FF f g a -> f a -> g (f a)
+jacobian Forward = F.jacobian
+jacobian Reverse = R.jacobian
+jacobian Tower = F.jacobian -- error "jacobian Tower: unimplemented"
+jacobian Mixed = M.jacobian
+{-# INLINE jacobian #-}
+
+jacobian' :: (Traversable f, Traversable g, Num a) => Direction -> FF f g a -> f a -> g (a, f a)
+jacobian' Forward = F.jacobian'
+jacobian' Reverse = R.jacobian'
+jacobian' Tower = F.jacobian' -- error "jacobian' Tower: unimplemented"
+jacobian' Mixed = M.jacobian'
+{-# INLINE jacobian' #-}
+
+grad :: (Traversable f, Num a) => Direction -> FU f a -> f a -> f a
+grad Forward = F.grad
+grad Reverse = R.grad
+grad Tower   = F.grad -- error "grad Tower: unimplemented"
+grad Mixed   = M.grad
+{-# INLINE grad #-}
+
+grad' :: (Traversable f, Num a) => Direction -> FU f a -> f a -> (a, f a)
+grad' Forward = F.grad'
+grad' Reverse = R.grad'
+grad' Tower   = F.grad' -- error "grad' Tower: unimplemented"
+grad' Mixed   = M.grad'
+{-# INLINE grad' #-}
+
diff --git a/Numeric/AD/Mode/Forward.hs b/Numeric/AD/Mode/Forward.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/AD/Mode/Forward.hs
@@ -0,0 +1,179 @@
+{-# LANGUAGE Rank2Types #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD.Mode.Forward
+-- Copyright   :  (c) Edward Kmett 2010
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only
+--
+-- Forward mode automatic differentiation
+--
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Mode.Forward
+    (
+    -- * Gradient
+      grad
+    , grad'
+    , gradWith
+    , gradWith'
+    -- * Jacobian
+    , jacobian
+    , jacobian'
+    , jacobianWith
+    , jacobianWith'
+    -- * Transposed Jacobian
+    , jacobianT
+    , jacobianWithT
+    -- * Hessian Product
+    , hessianProduct
+    , hessianProduct'
+    -- * Derivatives
+    , diff
+    , diff'
+    , diffF
+    , diffF'
+    -- * Directional Derivatives
+    , du
+    , du'
+    , duF
+    , duF'
+    -- * Monadic Combinators
+    , diffM
+    , diffM'
+    -- * Exposed Types
+    , UU, UF, FU, FF
+    , AD(..)
+    , Mode(..)
+    ) where
+
+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
+import Numeric.AD.Internal.Forward
+
+du :: (Functor f, Num a) => FU f a -> f (a, a) -> a
+du f = tangent . f . fmap (uncurry bundle)
+{-# INLINE du #-}
+
+du' :: (Functor f, Num a) => FU f a -> f (a, a) -> (a, a)
+du' f = unbundle . f . fmap (uncurry bundle)
+{-# INLINE du' #-}
+
+duF :: (Functor f, Functor g, Num a) => FF f g a -> f (a, a) -> g a
+duF f = fmap tangent . f . fmap (uncurry bundle)
+{-# INLINE duF #-}
+
+duF' :: (Functor f, Functor g, Num a) => FF f g 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 :: Num a => UU a -> a -> a
+diff f a = tangent $ apply f a
+{-# INLINE diff #-}
+
+-- | The 'd'UU' 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
+diff' :: Num a => UU 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'
+diffF :: (Functor f, Num a) => UF f 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'
+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)
+jacobianT f = bind (fmap tangent . f)
+{-# INLINE jacobianT #-}
+
+-- | A fast, simple transposed Jacobian computed with forward-mode AD.
+jacobianWithT :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> FF f g a -> f a -> f (g b)
+jacobianWithT g f = bindWith g' f
+    where g' a ga = g a . tangent <$> ga
+{-# INLINE jacobianWithT #-}
+
+jacobian :: (Traversable f, Traversable g, Num a) => FF f g a -> f a -> g (f a)
+jacobian f as = transposeWith (const id) t p
+    where
+        (p, t) = bind' (fmap tangent . f) as
+{-# INLINE jacobian #-}
+
+jacobianWith :: (Traversable f, Traversable g, Num a) => (a -> a -> b) -> FF f g a -> f a -> g (f b)
+jacobianWith g f as = transposeWith (const id) t p
+    where
+        (p, t) = bindWith' g' f as
+        g' a ga = g a . tangent <$> ga
+{-# INLINE jacobianWith #-}
+
+jacobian' :: (Traversable f, Traversable g, Num a) => FF f g a -> f a -> g (a, f a)
+jacobian' f as = transposeWith row t p
+    where
+        (p, t) = bind' f as
+        row x as' = (primal x, tangent <$> as')
+{-# INLINE jacobian' #-}
+
+jacobianWith' :: (Traversable f, Traversable g, Num a) => (a -> a -> b) -> FF f g a -> f a -> g (a, f b)
+jacobianWith' g f as = transposeWith row t p
+    where
+        (p, t) = bindWith' g' f as
+        row x as' = (primal x, as')
+        g' a ga = g a . tangent <$> ga
+{-# INLINE jacobianWith' #-}
+
+grad :: (Traversable f, Num a) => FU f a -> f a -> f a
+grad f = bind (tangent . f)
+{-# INLINE grad #-}
+
+grad' :: (Traversable f, Num a) => FU f a -> f a -> (a, f a)
+grad' f as = (primal b, tangent <$> bs)
+    where
+        (b, bs) = bind' f as
+{-# INLINE grad' #-}
+
+gradWith :: (Traversable f, Num a) => (a -> a -> b) -> FU f a -> f a -> f b
+gradWith g f = bindWith g (tangent . f)
+{-# INLINE gradWith #-}
+
+gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> FU f 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. 
+hessianProduct :: (Traversable f, Num a) => FU f 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. 
+hessianProduct' :: (Traversable f, Num a) => FU f a -> f (a, a) -> f (a, a)
+hessianProduct' f = duF' $ grad $ decomposeMode . f . fmap composeMode
+
+-- * Experimental
+
+-- data f :> a = a :< f (f :> a)
+-- gradients :: (Traversable f, Num a) => FU f a -> f a -> (f :> a)
diff --git a/Numeric/AD/Mode/Mixed.hs b/Numeric/AD/Mode/Mixed.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/AD/Mode/Mixed.hs
@@ -0,0 +1,215 @@
+{-# LANGUAGE Rank2Types, TypeFamilies #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD.Mode.Mixed
+-- Copyright   :  (c) Edward Kmett 2010
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only
+--
+-- Mixed-Mode Automatic Differentiation.
+--
+-- Each combinator exported from this module chooses an appropriate AD mode.
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Mode.Mixed
+    (
+    -- * Gradients (Reverse Mode)
+      grad
+    , grad'
+    , gradWith
+    , gradWith'
+
+    -- * Jacobians (Mixed Mode)
+    , jacobian
+    , jacobian'
+    , jacobianWith
+    , jacobianWith'
+
+    -- * Monadic Gradient/Jacobian (Reverse Mode)
+    , gradM
+    , gradM'
+    , gradWithM
+    , gradWithM'
+
+    -- * Functorial Gradient/Jacobian (Reverse Mode)
+    , gradF
+    , gradF'
+    , gradWithF
+    , gradWithF'
+
+    -- * Transposed Jacobians (Forward Mode)
+    , jacobianT
+    , jacobianWithT
+
+    -- * Hessian (Forward-On-Reverse)
+    , hessian
+
+    -- * Hessian Tensors (Forward-On-Mixed)
+    , hessianTensor
+
+    -- * Hessian Vector Products (Forward-On-Reverse)
+    , hessianProduct
+    , hessianProduct'
+
+    -- * Derivatives (Forward Mode)
+    , diff
+    , diffF
+
+    , diff'
+    , diffF'
+
+    -- * Derivatives (Tower)
+    , diffs
+    , diffsF
+
+    , diffs0
+    , diffs0F
+
+    -- * Directional Derivatives (Forward Mode)
+    , du
+    , du'
+    , duF
+    , duF'
+
+    -- * Directional Derivatives (Tower)
+    , dus
+    , dus0
+    , dusF
+    , dus0F
+
+    -- * Taylor Series (Tower)
+    , taylor
+    , taylor0
+
+    -- * Maclaurin Series (Tower)
+    , maclaurin
+    , maclaurin0
+
+    -- * Monadic Combinators (Forward Mode)
+    , diffM
+    , diffM'
+
+    -- * Exposed Types
+    , UU, UF, FU, FF
+    , AD(..)
+    , Mode(..)
+    ) where
+
+import Data.Traversable (Traversable)
+import Data.Foldable (Foldable, foldr')
+import Control.Applicative
+
+import Numeric.AD.Types (AD(..), UU, UF, FU, FF)
+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
+    ) 
+
+import Numeric.AD.Mode.Tower 
+    ( diffsF, diffs0F, diffs, diffs0
+    , taylor, taylor0, maclaurin, maclaurin0
+    , 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'
+    )
+
+-- | 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)
+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.
+--
+-- 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
+    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)
+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.
+--
+-- 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)
+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
+    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.
+hessianProduct :: (Traversable f, Num a) => FU f a -> f (a, a) -> f a
+hessianProduct f = duF (grad (decomposeMode . f . fmap composeMode))
+
+-- | @'hessianProduct'' f wv@ computes both the gradient of a non-scalar-to-scalar @f@ at @w = 'fst' <$> wv@ and the product of the hessian @H@ at @w@ with a vector @v = snd <$> wv@ using \"Pearlmutter's method\". The outputs are returned wrapped in the same functor.
+--
+-- > H v = (d/dr) grad_w (w + r v) | r = 0
+-- 
+-- Or in other words, we take the directional derivative of the gradient.
+-- 
+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.
+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
+-}
diff --git a/Numeric/AD/Mode/Reverse.hs b/Numeric/AD/Mode/Reverse.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/AD/Mode/Reverse.hs
@@ -0,0 +1,219 @@
+{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD.Mode.Reverse
+-- Copyright   :  (c) Edward Kmett 2010
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only
+--
+-- Mixed-Mode Automatic Differentiation.
+--
+-- For reverse mode AD we use 'System.Mem.StableName.StableName' to recover sharing information from
+-- the tape to avoid combinatorial explosion, and thus run asymptotically faster
+-- than it could without such sharing information, but the use of side-effects
+-- contained herein is benign.
+--
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Mode.Reverse
+    (
+    -- * Gradient
+      grad
+    , grad'
+    , gradWith
+    , gradWith'
+    -- * Jacobian
+    , jacobian
+    , jacobian'
+    , jacobianWith
+    , jacobianWith'
+    -- * Hessian
+    , hessian
+    , hessianM
+    , hessianTensor
+    
+    -- * Derivatives
+    , diff
+    , diff'
+    , diffF
+    , diffF'
+    -- * Monadic Combinators
+    , diffM
+    , diffM'
+    , gradM
+    , gradM'
+    , gradWithM
+    , gradWithM'
+    -- * Synonyms
+    , gradF
+    , gradF'
+    , gradWithF
+    , gradWithF'
+    -- * Exposed Types
+    , UU, UF, FU, FF
+    , AD(..)
+    , Mode(..)
+    ) where
+
+import Control.Monad (liftM)
+import Control.Applicative (WrappedMonad(..),(<$>))
+import Data.Traversable (Traversable)
+
+import Numeric.AD.Types
+import Numeric.AD.Internal.Classes
+import Numeric.AD.Internal.Composition
+import Numeric.AD.Internal.Reverse
+
+-- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with 'Reverse' AD in a single pass.
+grad :: (Traversable f, Num a) => FU f 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' :: (Traversable f, Num a) => FU f a -> f a -> (a, f a)
+grad' f as = (primal r, unbind vs $ partialArray bds r)
+    where (vs, bds) = bind as
+          r = f vs
+{-# 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) -> FU f a -> f a -> f b
+gradWith g f as = unbindWith g vs (partialArray bds $ f 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' 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) -> FU f a -> f a -> (a, f b)
+gradWith' g f as = (primal r, unbindWith g vs $ partialArray bds r)
+    where (vs, bds) = bind as
+          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'
+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' #-}
+
+-- | 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
+    (vs, bds) = bind as
+    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.
+--
+-- 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)
+--
+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 
+{-# 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'
+--
+-- Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@.
+--
+-- > jacobian' == gradWithF' (\_ 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
+    (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
+diff f a = derivative $ f (var a 0)
+{-# INLINE diff #-}
+
+-- | The 'd'' function calculates the value and derivative, as a
+-- pair, of a scalar-to-scalar function.
+diff' :: Num a => UU a -> a -> (a, a)
+diff' f a = derivative' $ f (var a 0)
+{-# INLINE diff' #-}
+
+diffF :: (Functor f, Num a) => UF f a -> a -> f a
+diffF f a = derivative <$> f (var a 0)
+{-# INLINE diffF #-}
+
+diffF' :: (Functor f, Num a) => UF f a -> a -> f (a, a)
+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".
+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.
+--
+-- 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))
+
+-- | 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/Tower.hs b/Numeric/AD/Mode/Tower.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/AD/Mode/Tower.hs
@@ -0,0 +1,140 @@
+{-# LANGUAGE Rank2Types, BangPatterns #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      : Numeric.AD.Mode.Tower
+-- 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.Tower
+    (
+    -- * Taylor Series
+      taylor
+    , taylor0
+    -- * Maclaurin Series
+    , maclaurin
+    , maclaurin0
+    -- * Derivatives
+    , diff    -- first derivative of (a -> a) 
+    , diff'   -- answer and first derivative of (a -> a) 
+    , diffs   -- answer and all derivatives of (a -> a) 
+    , diffs0  -- zero padded derivatives of (a -> a)
+    , diffsF  -- answer and all derivatives of (a -> f a)
+    , diffs0F -- zero padded derivatives of (a -> f a)
+    -- * Directional Derivatives
+    , du      -- directional derivative of (a -> a)
+    , du'     -- answer and directional derivative of (a -> a)
+    , dus     -- answer and all directional derivatives of (a -> a) 
+    , dus0    -- answer and all zero padded directional derivatives of (a -> a)
+    , duF     -- directional derivative of (a -> f a)
+    , 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
+import Numeric.AD.Internal.Tower
+
+diffs :: Num a => UU a -> a -> [a]
+diffs f a = getADTower $ apply f a
+{-# INLINE diffs #-}
+
+diffs0 :: Num a => UU a -> a -> [a]
+diffs0 f a = zeroPad (diffs f a)
+{-# INLINE diffs0 #-}
+
+diffsF :: (Functor f, Num a) => UF f a -> a -> f [a]
+diffsF f a = getADTower <$> apply f a
+{-# INLINE diffsF #-}
+
+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)
+    where
+        go !n !acc (a:as) = a * acc : go (n + 1) (acc * dx / n) as
+        go _ _ [] = []
+
+taylor0 :: Fractional a => UU a -> a -> a -> [a]
+taylor0 f x dx = zeroPad (taylor f x dx)
+{-# INLINE taylor0 #-}
+
+maclaurin :: Fractional a => UU a -> a -> [a]
+maclaurin f = taylor f 0
+{-# INLINE maclaurin #-}
+
+maclaurin0 :: Fractional a => UU a -> a -> [a]
+maclaurin0 f = taylor0 f 0
+{-# INLINE maclaurin0 #-}
+
+diff :: Num a => UU a -> a -> a
+diff f = d . diffs f
+{-# INLINE diff #-}
+
+diff' :: Num a => UU a -> a -> (a, a)
+diff' f = d' . diffs f
+{-# INLINE diff' #-}
+
+du :: (Functor f, Num a) => FU f a -> f (a, a) -> a
+du f = d . getADTower . f . fmap withD
+{-# INLINE du #-}
+
+du' :: (Functor f, Num a) => FU f a -> f (a, a) -> (a, a)
+du' f = d' . getADTower . f . fmap withD
+{-# INLINE du' #-}
+
+duF :: (Functor f, Functor g, Num a) => FF f g a -> f (a, a) -> g a
+duF f = fmap (d . getADTower) . f . fmap withD
+{-# INLINE duF #-}
+
+duF' :: (Functor f, Functor g, Num a) => FF f g a -> f (a, a) -> g (a, a)
+duF' f = fmap (d' . getADTower) . f . fmap withD
+{-# INLINE duF' #-}
+
+dus :: (Functor f, Num a) => FU f a -> f [a] -> [a]
+dus f = getADTower . f . fmap tower
+{-# INLINE dus #-}
+
+dus0 :: (Functor f, Num a) => FU f a -> f [a] -> [a]
+dus0 f = zeroPad . getADTower . f . fmap tower
+{-# INLINE dus0 #-}
+
+dusF :: (Functor f, Functor g, Num a) => FF f g a -> f [a] -> g [a]
+dusF f = fmap getADTower . f . fmap tower
+{-# INLINE dusF #-}
+
+dus0F :: (Functor f, Functor g, Num a) => FF f g a -> f [a] -> g [a]
+dus0F f = fmap getADTower . f . fmap tower
+{-# INLINE dus0F #-}
+
+-- TODO: higher order gradients
+-- data f :> a = a :< f (f :> a) 
+-- gradients  :: (Traversable f, Num a) => FU f a -> f a -> f :> a
+-- gradientsF, jacobians :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (f :> a)
+-- gradientsM :: (Traversable f, Monad m, Num a) => FF f m a -> f a -> m (f :> a)
diff --git a/Numeric/AD/Newton.hs b/Numeric/AD/Newton.hs
--- a/Numeric/AD/Newton.hs
+++ b/Numeric/AD/Newton.hs
@@ -35,11 +35,12 @@
 import Prelude hiding (all)
 import Control.Monad (liftM)
 import Data.MList
-import Numeric.AD.Internal
 import Data.Foldable (all)
 import Data.Traversable (Traversable)
-import Numeric.AD.Forward (diff, diff', diffM, diffM')
-import Numeric.AD.Reverse (gradWith', gradWithM')
+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.Internal.Composition
 
 -- | The 'findZero' function finds a zero of a scalar function using
diff --git a/Numeric/AD/Reverse.hs b/Numeric/AD/Reverse.hs
deleted file mode 100644
--- a/Numeric/AD/Reverse.hs
+++ /dev/null
@@ -1,218 +0,0 @@
-{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns #-}
------------------------------------------------------------------------------
--- |
--- Module      :  Numeric.AD.Reverse
--- Copyright   :  (c) Edward Kmett 2010
--- License     :  BSD3
--- Maintainer  :  ekmett@gmail.com
--- Stability   :  experimental
--- Portability :  GHC only
---
--- Mixed-Mode Automatic Differentiation.
---
--- For reverse mode AD we use 'System.Mem.StableName.StableName' to recover sharing information from
--- the tape to avoid combinatorial explosion, and thus run asymptotically faster
--- than it could without such sharing information, but the use of side-effects
--- contained herein is benign.
---
------------------------------------------------------------------------------
-
-module Numeric.AD.Reverse
-    (
-    -- * Gradient
-      grad
-    , grad'
-    , gradWith
-    , gradWith'
-    -- * Jacobian
-    , jacobian
-    , jacobian'
-    , jacobianWith
-    , jacobianWith'
-    -- * Hessian
-    , hessian
-    , hessianM
-    , hessianTensor
-    
-    -- * Derivatives
-    , diff
-    , diff'
-    , diffF
-    , diffF'
-    -- * Monadic Combinators
-    , diffM
-    , diffM'
-    , gradM
-    , gradM'
-    , gradWithM
-    , gradWithM'
-    -- * Synonyms
-    , gradF
-    , gradF'
-    , gradWithF
-    , gradWithF'
-    -- * Exposed Types
-    , UU, UF, FU, FF
-    , AD(..)
-    , Mode(..)
-    ) where
-
-import Control.Monad (liftM)
-import Control.Applicative (WrappedMonad(..),(<$>))
-import Data.Traversable (Traversable)
-
-import Numeric.AD.Internal
-import Numeric.AD.Internal.Composition
-import Numeric.AD.Internal.Reverse
-
--- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with 'Reverse' AD in a single pass.
-grad :: (Traversable f, Num a) => FU f 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' :: (Traversable f, Num a) => FU f a -> f a -> (a, f a)
-grad' f as = (primal r, unbind vs $ partialArray bds r)
-    where (vs, bds) = bind as
-          r = f vs
-{-# 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) -> FU f a -> f a -> f b
-gradWith g f as = unbindWith g vs (partialArray bds $ f 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' 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) -> FU f a -> f a -> (a, f b)
-gradWith' g f as = (primal r, unbindWith g vs $ partialArray bds r)
-    where (vs, bds) = bind as
-          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'
-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' #-}
-
--- | 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
-    (vs, bds) = bind as
-    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.
---
--- 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)
---
-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 
-{-# 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'
---
--- Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@.
---
--- > jacobian' == gradWithF' (\_ 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
-    (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
-diff f a = derivative $ f (var a 0)
-{-# INLINE diff #-}
-
--- | The 'd'' function calculates the value and derivative, as a
--- pair, of a scalar-to-scalar function.
-diff' :: Num a => UU a -> a -> (a, a)
-diff' f a = derivative' $ f (var a 0)
-{-# INLINE diff' #-}
-
-diffF :: (Functor f, Num a) => UF f a -> a -> f a
-diffF f a = derivative <$> f (var a 0)
-{-# INLINE diffF #-}
-
-diffF' :: (Functor f, Num a) => UF f a -> a -> f (a, a)
-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".
-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.
---
--- 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))
-
--- | 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/Tensors.hs b/Numeric/AD/Tensors.hs
deleted file mode 100644
--- a/Numeric/AD/Tensors.hs
+++ /dev/null
@@ -1,32 +0,0 @@
-{-# LANGUAGE TypeOperators, TemplateHaskell, ScopedTypeVariables #-}
------------------------------------------------------------------------------
--- |
--- Module      :  Numeric.AD.Tensors
--- Copyright   :  (c) Edward Kmett 2010
--- License     :  BSD3
--- Maintainer  :  ekmett@gmail.com
--- Stability   :  experimental
--- Portability :  GHC only
---
------------------------------------------------------------------------------
-
-module Numeric.AD.Tensors
-    ( 
-    -- * Tensors
-      Tensors(..)
-    , headT
-    , tailT
-    , tensors
-    -- * f-Branching Streams
-    , Stream(..)
-    , headS
-    , tailS
-    , unfoldS
-    -- * Comonads
-    , Copointed(..)
-    , Comonad(..)
-    ) where
-
-import Numeric.AD.Internal.Comonad
-import Numeric.AD.Internal.Stream
-import Numeric.AD.Internal.Tensors
diff --git a/Numeric/AD/Tower.hs b/Numeric/AD/Tower.hs
deleted file mode 100644
--- a/Numeric/AD/Tower.hs
+++ /dev/null
@@ -1,139 +0,0 @@
-{-# LANGUAGE Rank2Types, BangPatterns #-}
------------------------------------------------------------------------------
--- |
--- Module      : Numeric.AD.Tower
--- 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.Tower
-    (
-    -- * Taylor Series
-      taylor
-    , taylor0
-    -- * Maclaurin Series
-    , maclaurin
-    , maclaurin0
-    -- * Derivatives
-    , diff    -- first derivative of (a -> a) 
-    , diff'   -- answer and first derivative of (a -> a) 
-    , diffs   -- answer and all derivatives of (a -> a) 
-    , diffs0  -- zero padded derivatives of (a -> a)
-    , diffsF  -- answer and all derivatives of (a -> f a)
-    , diffs0F -- zero padded derivatives of (a -> f a)
-    -- * Directional Derivatives
-    , du      -- directional derivative of (a -> a)
-    , du'     -- answer and directional derivative of (a -> a)
-    , dus     -- answer and all directional derivatives of (a -> a) 
-    , dus0    -- answer and all zero padded directional derivatives of (a -> a)
-    , duF     -- directional derivative of (a -> f a)
-    , 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.Internal
-import Numeric.AD.Internal.Tower
-
-diffs :: Num a => UU a -> a -> [a]
-diffs f a = getADTower $ apply f a
-{-# INLINE diffs #-}
-
-diffs0 :: Num a => UU a -> a -> [a]
-diffs0 f a = zeroPad (diffs f a)
-{-# INLINE diffs0 #-}
-
-diffsF :: (Functor f, Num a) => UF f a -> a -> f [a]
-diffsF f a = getADTower <$> apply f a
-{-# INLINE diffsF #-}
-
-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)
-    where
-        go !n !acc (a:as) = a * acc : go (n + 1) (acc * dx / n) as
-        go _ _ [] = []
-
-taylor0 :: Fractional a => UU a -> a -> a -> [a]
-taylor0 f x dx = zeroPad (taylor f x dx)
-{-# INLINE taylor0 #-}
-
-maclaurin :: Fractional a => UU a -> a -> [a]
-maclaurin f = taylor f 0
-{-# INLINE maclaurin #-}
-
-maclaurin0 :: Fractional a => UU a -> a -> [a]
-maclaurin0 f = taylor0 f 0
-{-# INLINE maclaurin0 #-}
-
-diff :: Num a => UU a -> a -> a
-diff f = d . diffs f
-{-# INLINE diff #-}
-
-diff' :: Num a => UU a -> a -> (a, a)
-diff' f = d' . diffs f
-{-# INLINE diff' #-}
-
-du :: (Functor f, Num a) => FU f a -> f (a, a) -> a
-du f = d . getADTower . f . fmap withD
-{-# INLINE du #-}
-
-du' :: (Functor f, Num a) => FU f a -> f (a, a) -> (a, a)
-du' f = d' . getADTower . f . fmap withD
-{-# INLINE du' #-}
-
-duF :: (Functor f, Functor g, Num a) => FF f g a -> f (a, a) -> g a
-duF f = fmap (d . getADTower) . f . fmap withD
-{-# INLINE duF #-}
-
-duF' :: (Functor f, Functor g, Num a) => FF f g a -> f (a, a) -> g (a, a)
-duF' f = fmap (d' . getADTower) . f . fmap withD
-{-# INLINE duF' #-}
-
-dus :: (Functor f, Num a) => FU f a -> f [a] -> [a]
-dus f = getADTower . f . fmap tower
-{-# INLINE dus #-}
-
-dus0 :: (Functor f, Num a) => FU f a -> f [a] -> [a]
-dus0 f = zeroPad . getADTower . f . fmap tower
-{-# INLINE dus0 #-}
-
-dusF :: (Functor f, Functor g, Num a) => FF f g a -> f [a] -> g [a]
-dusF f = fmap getADTower . f . fmap tower
-{-# INLINE dusF #-}
-
-dus0F :: (Functor f, Functor g, Num a) => FF f g a -> f [a] -> g [a]
-dus0F f = fmap getADTower . f . fmap tower
-{-# INLINE dus0F #-}
-
--- TODO: higher order gradients
--- data f :> a = a :< f (f :> a) 
--- gradients  :: (Traversable f, Num a) => FU f a -> f a -> f :> a
--- gradientsF, jacobians :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (f :> a)
--- gradientsM :: (Traversable f, Monad m, Num a) => FF f m a -> f a -> m (f :> a)
diff --git a/Numeric/AD/Types.hs b/Numeric/AD/Types.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/AD/Types.hs
@@ -0,0 +1,31 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD.Types
+-- Copyright   :  (c) Edward Kmett 2010
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only
+--
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Types
+    ( 
+      AD(..)
+    -- * Differentiable Functions
+    , UU, UF, FU, FF
+    -- * Tensors
+    , Tensors(..)
+    , headT
+    , tailT
+    , tensors
+    -- * f-Branching Streams
+    , Stream(..)
+    , headS
+    , tailS
+    , unfoldS
+    ) where
+
+import Numeric.AD.Internal.Types
+import Numeric.AD.Internal.Stream
+import Numeric.AD.Internal.Tensors
diff --git a/TODO b/TODO
--- a/TODO
+++ b/TODO
@@ -1,13 +1,43 @@
-* Implement the diffMF, etc. functionality from Numeric.FAD
+* Use the f-branching stream and tensors to generate:
 
+    gradients :: Traversable f => FU f a -> f a -> Stream f a 
+
 * Allow the type to vary within our AD data type container, in the same fashion as Numeric.FAD.
 
-    Although, while Pearlmutter and Siskind provided the functionality to permit it in derivative combinator, they provided
-    no combinators to convert, for instance, a Dual tag Float to a Dual tag Double, so that extra functionality cannot
+    Although, while Pearlmutter and Siskind provided the functionality to permit 
+    it in the derivative combinators, they provided no combinators to convert, 
+    say, @Dual tag Float@ to a @Dual tag Double@, so that extra functionality cannot
     currently be leveraged.
 
-* Provide the ability to use reverse mode locally on FAD inputs, i.e.
+  One approach: GADT'd Tape. 
+    Lets us use local matrix-valued jacobians as blackboxes.
+    However, this requires a custom higher-order data-reify.
 
-    reverseCheckpoint :: (forall s. AD s a -> AD s a) -> AD t a -> AD t a 
+* Do we need some kind of Array implementation? These'd be easy:
 
-* Provide forward-on-reverse computation of Hessians
+    (new)?type ADArray s i e = Array i (AD s e) 
+    (new)?type ADIOArray s i e = IOArray i (AD s e)
+    (new)?type ADSTArray s i e = STArray i (AD s e)
+
+  But how to handle possibly unboxed arrays where possible?
+
+    Mode s => ADUArray s i Float
+    Mode s => ADIOUArray s i Float
+    Mode s => ADSTUArray s i Float
+
+  Forward admits fast unboxed array access, but how to avoid violating quantification over s?
+
+  This seems to require baking the array type for Float and Double directly into the Mode.
+
+  class ( IArray (ADArray s) e
+        , MArray (ADSTArray s s') e (ST s')
+        , MArray (ADIOArray s) e IO
+        , Lifted s) => LiftedArray s e where
+
+     type ADArr s e  :: * -> * -> * 
+     type ADSTArr s e :: * -> * -> * 
+     type ADIOArr s e :: * -> * -> *
+
+  newtype ADArray s i e = ADArray (ADArr s e i e)
+
+  class ( LiftedArray s Float, LiftedArray s Double, ...) => Mode s
diff --git a/ad.cabal b/ad.cabal
--- a/ad.cabal
+++ b/ad.cabal
@@ -1,5 +1,5 @@
 Name:         ad
-Version:      0.31.0
+Version:      0.32.0
 License:      BSD3
 License-File: LICENSE
 Copyright:    (c) Edward Kmett 2010,
@@ -30,14 +30,16 @@
 
 Exposed-Modules:
     Numeric.AD
-    Numeric.AD.Forward
-    Numeric.AD.Reverse
-    Numeric.AD.Tower
-    Numeric.AD.Directed
+    Numeric.AD.Classes
+    Numeric.AD.Types
     Numeric.AD.Newton
-    Numeric.AD.Tensors
 
-    Numeric.AD.Internal
+    Numeric.AD.Mode.Mixed
+    Numeric.AD.Mode.Forward
+    Numeric.AD.Mode.Reverse
+    Numeric.AD.Mode.Tower
+    Numeric.AD.Mode.Directed
+
     Numeric.AD.Internal.Classes
     Numeric.AD.Internal.Types
     Numeric.AD.Internal.Combinators
@@ -55,3 +57,4 @@
 Extra-Source-Files: TODO
 GHC-Options: -Wall -fspec-constr 
 --           -O2
+