diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown
--- a/CHANGELOG.markdown
+++ b/CHANGELOG.markdown
@@ -1,3 +1,9 @@
+3.3
+---
+* Renamed `Reverse` to `Kahn` and `Wengert` to `Reverse`. We use Arthur Kahn's topological sorting algorithm to
+  sort the tape after the fact in Kahn mode, while the stock Reverse mode builds a Wengert list as it goes, which
+  is more efficient in practice.
+
 3.2.2
 -----
 * Export of the `conjugateGradientDescent` and `gradientDescent` from `Numeric.AD`
diff --git a/README.markdown b/README.markdown
--- a/README.markdown
+++ b/README.markdown
@@ -87,7 +87,7 @@
  * `Numeric.AD.Mode.Forward` provides basic forward-mode AD. It is good for computing simple derivatives.
  * `Numeric.AD.Mode.Sparse` computes a sparse forward-mode AD tower. It is good for higher derivatives or large numbers of outputs.
  * `Numeric.AD.Mode.Reverse` computes with reverse-mode AD. It is good for computing a few outputs given many inputs.
- * `Numeric.AD.Mode.Chain` computes with reverse-mode AD. It is good for computing a few outputs given many inputs, when not using sparks.
+ * `Numeric.AD.Mode.Wengert` computes with reverse-mode AD. It is good for computing a few outputs given many inputs, when not using sparks.
  * `Numeric.AD.Mode.Tower` computes a dense forward-mode AD tower useful for higher derivatives of single input functions.
 
  * `Numeric.AD.Newton` provides a number of combinators for root finding using Newton's method with quadratic convergence.
diff --git a/ad.cabal b/ad.cabal
--- a/ad.cabal
+++ b/ad.cabal
@@ -1,5 +1,5 @@
 name:         ad
-version:      3.2.2
+version:      3.3.0.1
 license:      BSD3
 license-File: LICENSE
 copyright:    (c) Edward Kmett 2010-2012,
@@ -24,9 +24,9 @@
     .
     * @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. It generates a tree-like tape that needs to be topologically sorted in the end.
+    * @Numeric.AD.Mode.Reverse@ uses benign side-effects to compute reverse-mode AD. It is good for computing gradients in one pass. It generates a Wengert list (linear tape) using @Data.Reflection@.
     .
-    * @Numeric.AD.Mode.Wengert@ uses benign side-effects to compute reverse-mode AD. It is good for computing gradients in one pass. It generates a Wengert list (linear tape) using @Data.Reflection@.
+    * @Numeric.AD.Mode.Kahn@ uses benign side-effects to compute reverse-mode AD. It is good for computing gradients in one pass. It generates a tree-like tape that needs to be topologically sorted in the end.
     .
     * @Numeric.AD.Mode.Sparse@ computes a sparse forward-mode AD tower. It is good for higher derivatives or large numbers of outputs.
     .
@@ -91,7 +91,7 @@
     comonad          == 3.0.*,
     containers       >= 0.2 && < 0.6,
     data-reify       >= 0.6 && < 0.7,
-    free             >= 3.0 && <= 3.3,
+    free             >= 3.0 && <= 3.4,
     mtl,
     reflection       >= 1.1.6 && < 1.2,
     tagged           >= 0.4.2.1 && < 0.5,
@@ -104,24 +104,24 @@
     Numeric.AD.Newton
     Numeric.AD.Halley
 
-    Numeric.AD.Mode.Wengert
     Numeric.AD.Mode.Directed
     Numeric.AD.Mode.Forward
+    Numeric.AD.Mode.Kahn
     Numeric.AD.Mode.Reverse
     Numeric.AD.Mode.Tower
     Numeric.AD.Mode.Sparse
 
     Numeric.AD.Variadic
-    Numeric.AD.Variadic.Reverse
+    Numeric.AD.Variadic.Kahn
     Numeric.AD.Variadic.Sparse
 
     Numeric.AD.Internal.Classes
     Numeric.AD.Internal.Combinators
     Numeric.AD.Internal.Forward
     Numeric.AD.Internal.Tower
+    Numeric.AD.Internal.Kahn
     Numeric.AD.Internal.Reverse
     Numeric.AD.Internal.Var
-    Numeric.AD.Internal.Wengert
     Numeric.AD.Internal.Sparse
     Numeric.AD.Internal.Dense
     Numeric.AD.Internal.Composition
diff --git a/src/Numeric/AD/Internal/Kahn.hs b/src/Numeric/AD/Internal/Kahn.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/AD/Internal/Kahn.hs
@@ -0,0 +1,264 @@
+{-# LANGUAGE Rank2Types, TypeFamilies, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, FlexibleContexts, TemplateHaskell, UndecidableInstances, DeriveDataTypeable #-}
+-- {-# OPTIONS_HADDOCK hide, prune #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD.Internal.Kahn
+-- Copyright   :  (c) Edward Kmett 2010
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only
+--
+-- This module provides reverse-mode Automatic Differentiation implementation using
+-- linear time topological sorting after the fact.
+--
+-- For this form of 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.Internal.Kahn
+    ( Kahn(..)
+    , Tape(..)
+    , partials
+    , partialArray
+    , partialMap
+    , derivative
+    , derivative'
+    , vgrad, vgrad'
+    , Grad(..)
+    ) where
+
+import Prelude hiding (mapM)
+import Control.Applicative (Applicative(..),(<$>))
+import Control.Monad.ST
+import Control.Monad (forM_)
+import Data.List (foldl')
+import Data.Array.ST
+import Data.Array
+import Data.IntMap (IntMap, fromListWith)
+import Data.Graph (Vertex, transposeG, Graph)
+import Data.Reify (reifyGraph, MuRef(..))
+import qualified Data.Reify.Graph as Reified
+import System.IO.Unsafe (unsafePerformIO)
+import Language.Haskell.TH
+import Data.Data (Data)
+import Data.Typeable (Typeable)
+import Numeric.AD.Internal.Types
+import Numeric.AD.Internal.Classes
+import Numeric.AD.Internal.Identity
+import Numeric.AD.Internal.Var
+
+-- | A @Tape@ records the information needed back propagate from the output to each input during reverse 'Mode' AD.
+data Tape a t
+    = Zero
+    | Lift !a
+    | Var !a {-# UNPACK #-} !Int
+    | Binary !a a a t t
+    | Unary !a a t
+    deriving (Show, Data, Typeable)
+
+-- | @Kahn@ is a 'Mode' using reverse-mode automatic differentiation that provides fast 'diffFU', 'diff2FU', 'grad', 'grad2' and a fast 'jacobian' when you have a significantly smaller number of outputs than inputs.
+newtype Kahn a = Kahn (Tape a (Kahn a)) deriving (Show, Typeable)
+
+-- deriving instance (Data (Tape a (Kahn a)) => Data (Kahn a)
+
+instance MuRef (Kahn a) where
+    type DeRef (Kahn a) = Tape a
+
+    mapDeRef _ (Kahn Zero) = pure Zero
+    mapDeRef _ (Kahn (Lift a)) = pure (Lift a)
+    mapDeRef _ (Kahn (Var a v)) = pure (Var a v)
+    mapDeRef f (Kahn (Binary a dadb dadc b c)) = Binary a dadb dadc <$> f b <*> f c
+    mapDeRef f (Kahn (Unary a dadb b)) = Unary a dadb <$> f b
+
+instance Lifted Kahn => Mode Kahn where
+    isKnownZero (Kahn Zero) = True
+    isKnownZero _    = False
+
+    isKnownConstant (Kahn Zero) = True
+    isKnownConstant (Kahn (Lift _)) = True
+    isKnownConstant _ = False
+
+    auto a = Kahn (Lift a)
+    zero   = Kahn Zero
+    (<+>)  = binary (+) one one
+    a *^ b = lift1 (a *) (\_ -> auto a) b
+    a ^* b = lift1 (* b) (\_ -> auto b) a
+    a ^/ b = lift1 (/ b) (\_ -> auto (recip b)) a
+
+    Kahn Zero <**> y                = auto (0 ** primal y)
+    _            <**> Kahn Zero     = auto 1
+    x            <**> Kahn (Lift y) = lift1 (**y) (\z -> (y *^ z ** Id (y-1))) x
+    x            <**> y                = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y
+
+instance Primal Kahn where
+    primal (Kahn Zero) = 0
+    primal (Kahn (Lift a)) = a
+    primal (Kahn (Var a _)) = a
+    primal (Kahn (Binary a _ _ _ _)) = a
+    primal (Kahn (Unary a _ _)) = a
+
+instance Lifted Kahn => Jacobian Kahn where
+    type D Kahn = Id
+
+    unary f _         (Kahn Zero)     = Kahn (Lift (f 0))
+    unary f _         (Kahn (Lift a)) = Kahn (Lift (f a))
+    unary f (Id dadb) b                  = Kahn (Unary (f (primal b)) dadb b)
+
+    lift1 f df b = unary f (df (Id pb)) b
+        where pb = primal b
+
+    lift1_ f df b = unary (const a) (df (Id a) (Id pb)) b
+        where pb = primal b
+              a = f pb
+
+    binary f _         _         (Kahn Zero)     (Kahn Zero)     = Kahn (Lift (f 0 0))
+    binary f _         _         (Kahn Zero)     (Kahn (Lift c)) = Kahn (Lift (f 0 c))
+    binary f _         _         (Kahn (Lift b)) (Kahn Zero)     = Kahn (Lift (f b 0))
+    binary f _         _         (Kahn (Lift b)) (Kahn (Lift c)) = Kahn (Lift (f b c))
+    binary f _         (Id dadc) (Kahn Zero)     c                  = Kahn (Unary (f 0 (primal c)) dadc c)
+    binary f _         (Id dadc) (Kahn (Lift b)) c                  = Kahn (Unary (f b (primal c)) dadc c)
+    binary f (Id dadb) _         b                  (Kahn Zero)     = Kahn (Unary (f (primal b) 0) dadb b)
+    binary f (Id dadb) _         b                  (Kahn (Lift c)) = Kahn (Unary (f (primal b) c) dadb b)
+    binary f (Id dadb) (Id dadc) b                  c                  = Kahn (Binary (f (primal b) (primal c)) dadb dadc b c)
+
+    lift2 f df b c = binary f dadb dadc b c
+        where (dadb, dadc) = df (Id (primal b)) (Id (primal c))
+
+    lift2_ f df b c = binary (\_ _ -> a) dadb dadc b c
+        where
+            pb = primal b
+            pc = primal c
+            a = f pb pc
+            (dadb, dadc) = df (Id a) (Id pb) (Id pc)
+
+deriveLifted id (conT ''Kahn)
+
+derivative :: Num a => AD Kahn a -> a
+derivative = sum . map snd . partials
+{-# INLINE derivative #-}
+
+derivative' :: Num a => AD Kahn a -> (a, a)
+derivative' r = (primal r, derivative r)
+{-# INLINE derivative' #-}
+
+-- | back propagate sensitivities along a tape.
+backPropagate :: Num a => (Vertex -> (Tape a Int, Int, [Int])) -> STArray s Int a -> Vertex -> ST s ()
+backPropagate vmap ss v = do
+        case node of
+            Unary _ g b -> do
+                da <- readArray ss i
+                db <- readArray ss b
+                writeArray ss b (db + g*da)
+            Binary _ gb gc b c -> do
+                da <- readArray ss i
+                db <- readArray ss b
+                writeArray ss b (db + gb*da)
+                dc <- readArray ss c
+                writeArray ss c (dc + gc*da)
+            _ -> return ()
+    where
+        (node, i, _) = vmap v
+        -- this isn't _quite_ right, as it should allow negative zeros to multiply through
+
+topSortAcyclic :: Graph -> [Vertex]
+topSortAcyclic g = reverse $ runST $ do
+    del <- newArray (bounds g) False :: ST s (STUArray s Int Bool)
+    let tg = transposeG g
+        starters = [ n | (n, []) <- assocs tg ]
+        loop [] rs = return rs
+        loop (n:ns) rs = do
+            writeArray del n True
+            let add [] = return ns
+                add (m:ms) = do
+                    b <- ok (tg!m)
+                    ms' <- add ms
+                    if b then return (m:ms') else return ms'
+                ok [] = return True
+                ok (x:xs) = do b <- readArray del x; if b then ok xs else return False
+            ns' <- add (g!n)
+            loop ns' (n : rs)
+    loop starters []
+
+-- | This returns a list of contributions to the partials.
+-- The variable ids returned in the list are likely /not/ unique!
+{-# SPECIALIZE partials :: AD Kahn Double -> [(Int, Double)] #-}
+partials :: forall a . Num a => AD Kahn a -> [(Int, a)]
+partials (AD tape) = [ let v = sensitivities ! ix in seq v (ident, v) | (ix, Var _ ident) <- xs ]
+    where
+        Reified.Graph xs start = unsafePerformIO $ reifyGraph tape
+        g = array xsBounds [ (i, successors t) | (i, t) <- xs ]
+        vertexMap = array xsBounds xs
+        vmap i = (vertexMap ! i, i, [])
+        xsBounds = sbounds xs
+
+        sensitivities = runSTArray $ do
+            ss <- newArray xsBounds 0
+            writeArray ss start 1
+            forM_ (topSortAcyclic g) $
+                backPropagate vmap ss
+            return ss
+
+        sbounds ((a,_):as) = foldl' (\(lo,hi) (b,_) -> let lo' = min lo b; hi' = max hi b in lo' `seq` hi' `seq` (lo', hi')) (a,a) as
+        sbounds _ = undefined -- the graph can't be empty, it contains the output node!
+
+        successors :: Tape a t -> [t]
+        successors (Unary _ _ b) = [b]
+        successors (Binary _ _ _ b c) = [b,c]
+        successors _ = []
+
+-- | Return an 'Array' of 'partials' given bounds for the variable IDs.
+partialArray :: Num a => (Int, Int) -> AD Kahn a -> Array Int a
+partialArray vbounds tape = accumArray (+) 0 vbounds (partials tape)
+{-# INLINE partialArray #-}
+
+-- | Return an 'IntMap' of sparse partials
+partialMap :: Num a => AD Kahn a -> IntMap a
+partialMap = fromListWith (+) . partials
+{-# INLINE partialMap #-}
+
+-- A simple fresh variable supply monad
+newtype S a = S { runS :: Int -> (a,Int) }
+instance Monad S where
+    return a = S (\s -> (a,s))
+    S g >>= f = S (\s -> let (a,s') = g s in runS (f a) s')
+
+instance Var Kahn where
+    var a v = Kahn (Var a v)
+    varId (Kahn (Var _ v)) = v
+    varId _ = error "varId: not a Var"
+
+class Num a => Grad i o o' a | i -> a o o', o -> a i o', o' -> a i o where
+    pack :: i -> [AD Kahn a] -> AD Kahn a
+    unpack :: ([a] -> [a]) -> o
+    unpack' :: ([a] -> (a, [a])) -> o'
+
+instance Num a => Grad (AD Kahn a) [a] (a, [a]) a where
+    pack i _ = i
+    unpack f = f []
+    unpack' f = f []
+
+instance Grad i o o' a => Grad (AD Kahn 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/src/Numeric/AD/Internal/Reverse.hs b/src/Numeric/AD/Internal/Reverse.hs
--- a/src/Numeric/AD/Internal/Reverse.hs
+++ b/src/Numeric/AD/Internal/Reverse.hs
@@ -1,112 +1,137 @@
-{-# LANGUAGE Rank2Types, TypeFamilies, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, FlexibleContexts, TemplateHaskell, UndecidableInstances, DeriveDataTypeable #-}
+{-# LANGUAGE Rank2Types, TypeFamilies, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, FlexibleContexts, TemplateHaskell, UndecidableInstances, DeriveDataTypeable, GADTs, ScopedTypeVariables #-}
 -- {-# OPTIONS_HADDOCK hide, prune #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.AD.Internal.Reverse
--- Copyright   :  (c) Edward Kmett 2010
+-- Copyright   :  (c) Edward Kmett 2012
 -- License     :  BSD3
 -- Maintainer  :  ekmett@gmail.com
 -- Stability   :  experimental
 -- Portability :  GHC only
 --
--- Reverse-Mode Automatic Differentiation implementation details
+-- Reverse-Mode Automatic Differentiation using a single Wengert list (or \"tape\").
 --
--- 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.
+-- This version uses @Data.Reflection@ to find and update the tape.
 --
+-- This is asymptotically faster than using @Reverse@, which
+-- is forced to reify and topologically sort the graph, but it requires
+-- a fairly expensive rendezvous during construction when updated using
+-- multiple threads.
+--
 -----------------------------------------------------------------------------
 
 module Numeric.AD.Internal.Reverse
     ( Reverse(..)
     , Tape(..)
+    , Head(..)
+    , Cells(..)
+    , reifyTape
     , partials
-    , partialArray
-    , partialMap
-    , derivative
-    , derivative'
-    , vgrad, vgrad'
-    , Grad(..)
+    , partialArrayOf
+    , partialMapOf
+    , derivativeOf
+    , derivativeOf'
     ) where
 
-import Prelude hiding (mapM)
-import Control.Applicative (Applicative(..),(<$>))
 import Control.Monad.ST
-import Control.Monad (forM_)
-import Data.List (foldl')
 import Data.Array.ST
 import Data.Array
-import Data.IntMap (IntMap, fromListWith)
-import Data.Graph (Vertex, transposeG, Graph)
-import Data.Reify (reifyGraph, MuRef(..))
-import qualified Data.Reify.Graph as Reified
-import System.IO.Unsafe (unsafePerformIO)
-import Language.Haskell.TH
-import Data.Data (Data)
-import Data.Typeable (Typeable)
+import Data.Array.Unsafe as Unsafe
+import Data.IORef
+import Data.IntMap (IntMap, fromDistinctAscList)
+import Data.Proxy
+import Data.Reflection
+import Data.Typeable
+import Language.Haskell.TH hiding (reify)
 import Numeric.AD.Internal.Types
 import Numeric.AD.Internal.Classes
 import Numeric.AD.Internal.Identity
 import Numeric.AD.Internal.Var
+import Prelude hiding (mapM)
+import System.IO.Unsafe (unsafePerformIO)
+import Unsafe.Coerce
 
--- | A @Tape@ records the information needed back propagate from the output to each input during 'Reverse' 'Mode' AD.
-data Tape a t
-    = Zero
-    | Lift !a
-    | Var !a {-# UNPACK #-} !Int
-    | Binary !a a a t t
-    | Unary !a a t
-    deriving (Show, Data, Typeable)
+-- evil untyped tape
+data Cells where
+  Nil    :: Cells
+  Unary  :: {-# UNPACK #-} !Int -> a -> Cells -> Cells
+  Binary :: {-# UNPACK #-} !Int -> {-# UNPACK #-} !Int -> a -> a -> Cells -> Cells
 
--- | @Reverse@ is a 'Mode' using reverse-mode automatic differentiation that provides fast 'diffFU', 'diff2FU', 'grad', 'grad2' and a fast 'jacobian' when you have a significantly smaller number of outputs than inputs.
-newtype Reverse a = Reverse (Tape a (Reverse a)) deriving (Show, Typeable)
+dropCells :: Int -> Cells -> Cells
+dropCells 0 xs = xs
+dropCells _ Nil = Nil
+dropCells n (Unary _ _ xs)      = (dropCells $! n - 1) xs
+dropCells n (Binary _ _ _ _ xs) = (dropCells $! n - 1) xs
 
--- deriving instance (Data (Tape a (Reverse a)) => Data (Reverse a)
+data Head = Head {-# UNPACK #-} !Int Cells
 
-instance MuRef (Reverse a) where
-    type DeRef (Reverse a) = Tape a
+newtype Tape = Tape { getTape :: IORef Head }
 
-    mapDeRef _ (Reverse Zero) = pure Zero
-    mapDeRef _ (Reverse (Lift a)) = pure (Lift a)
-    mapDeRef _ (Reverse (Var a v)) = pure (Var a v)
-    mapDeRef f (Reverse (Binary a dadb dadc b c)) = Binary a dadb dadc <$> f b <*> f c
-    mapDeRef f (Reverse (Unary a dadb b)) = Unary a dadb <$> f b
+un :: Int -> a -> Head -> (Head, Int)
+un i di (Head r t) = h `seq` r' `seq` (h, r') where
+  r' = r + 1
+  h = Head r' (Unary i di t)
+{-# INLINE un #-}
 
-instance Lifted Reverse => Mode Reverse where
-    isKnownZero (Reverse Zero) = True
-    isKnownZero _    = False
+bin :: Int -> Int -> a -> a -> Head -> (Head, Int)
+bin i j di dj (Head r t) = h `seq` r' `seq` (h, r') where
+  r' = r + 1
+  h = Head r' (Binary i j di dj t)
+{-# INLINE bin #-}
 
-    isKnownConstant (Reverse Zero) = True
-    isKnownConstant (Reverse (Lift _)) = True
-    isKnownConstant _ = False
+modifyTape :: Reifies s Tape => p s -> (Head -> (Head, r)) -> IO r
+modifyTape p = atomicModifyIORef (getTape (reflect p))
+{-# INLINE modifyTape #-}
 
-    auto a = Reverse (Lift a)
-    zero   = Reverse Zero
-    (<+>)  = binary (+) one one
-    a *^ b = lift1 (a *) (\_ -> auto a) b
-    a ^* b = lift1 (* b) (\_ -> auto b) a
-    a ^/ b = lift1 (/ b) (\_ -> auto (recip b)) a
+-- | This is used to create a new entry on the chain given a unary function, its derivative with respect to its input,
+-- the variable ID of its input, and the value of its input. Used by 'unary' and 'binary' internally.
+unarily :: forall s a. Reifies s Tape => (a -> a) -> a -> Int -> a -> Reverse s a
+unarily f di i b = Reverse (unsafePerformIO (modifyTape (Proxy :: Proxy s) (un i di))) $! f b
+{-# INLINE unarily #-}
 
-    Reverse Zero <**> y                = auto (0 ** primal y)
-    _            <**> Reverse Zero     = auto 1
-    x            <**> Reverse (Lift y) = lift1 (**y) (\z -> (y *^ z ** Id (y-1))) x
-    x            <**> y                = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y
+-- | This is used to create a new entry on the chain given a binary function, its derivatives with respect to its inputs,
+-- their variable IDs and values. Used by 'binary' internally.
+binarily :: forall s a. Reifies s Tape => (a -> a -> a) -> a -> a -> Int -> a -> Int -> a -> Reverse s a
+binarily f di dj i b j c = Reverse (unsafePerformIO (modifyTape (Proxy :: Proxy s) (bin i j di dj))) $! f b c
+{-# INLINE binarily #-}
 
-instance Primal Reverse where
-    primal (Reverse Zero) = 0
-    primal (Reverse (Lift a)) = a
-    primal (Reverse (Var a _)) = a
-    primal (Reverse (Binary a _ _ _ _)) = a
-    primal (Reverse (Unary a _ _)) = a
+data Reverse s a where
+  Zero :: Reverse s a
+  Lift :: a -> Reverse s a
+  Reverse :: {-# UNPACK #-} !Int -> a -> Reverse s a
+  deriving (Show, Typeable)
 
-instance Lifted Reverse => Jacobian Reverse where
-    type D Reverse = Id
+instance (Reifies s Tape, Lifted (Reverse s)) => Mode (Reverse s) where
+  isKnownZero Zero = True
+  isKnownZero _    = False
 
-    unary f _         (Reverse Zero)     = Reverse (Lift (f 0))
-    unary f _         (Reverse (Lift a)) = Reverse (Lift (f a))
-    unary f (Id dadb) b                  = Reverse (Unary (f (primal b)) dadb b)
+  isKnownConstant Reverse{} = False
+  isKnownConstant _ = True
 
+  auto = Lift
+  zero = Zero
+  (<+>)  = binary (+) one one
+  a *^ b = lift1 (a *) (\_ -> auto a) b
+  a ^* b = lift1 (* b) (\_ -> auto b) a
+  a ^/ b = lift1 (/ b) (\_ -> auto (recip b)) a
+
+  Zero <**> y      = auto (0 ** primal y)
+  _    <**> Zero   = auto 1
+  x    <**> Lift y = lift1 (**y) (\z -> (y *^ z ** Id (y-1))) x
+  x    <**> y      = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y
+
+instance Primal (Reverse s) where
+    primal Zero = 0
+    primal (Lift a) = a
+    primal (Reverse _ a) = a
+
+instance (Reifies s Tape, Lifted (Reverse s)) => Jacobian (Reverse s) where
+    type D (Reverse s) = Id
+
+    unary f _         (Zero)   = Lift (f 0)
+    unary f _         (Lift a) = Lift (f a)
+    unary f (Id dadi) (Reverse i b) = unarily f dadi i b
+
     lift1 f df b = unary f (df (Id pb)) b
         where pb = primal b
 
@@ -114,16 +139,17 @@
         where pb = primal b
               a = f pb
 
-    binary f _         _         (Reverse Zero)     (Reverse Zero)     = Reverse (Lift (f 0 0))
-    binary f _         _         (Reverse Zero)     (Reverse (Lift c)) = Reverse (Lift (f 0 c))
-    binary f _         _         (Reverse (Lift b)) (Reverse Zero)     = Reverse (Lift (f b 0))
-    binary f _         _         (Reverse (Lift b)) (Reverse (Lift c)) = Reverse (Lift (f b c))
-    binary f _         (Id dadc) (Reverse Zero)     c                  = Reverse (Unary (f 0 (primal c)) dadc c)
-    binary f _         (Id dadc) (Reverse (Lift b)) c                  = Reverse (Unary (f b (primal c)) dadc c)
-    binary f (Id dadb) _         b                  (Reverse Zero)     = Reverse (Unary (f (primal b) 0) dadb b)
-    binary f (Id dadb) _         b                  (Reverse (Lift c)) = Reverse (Unary (f (primal b) c) dadb b)
-    binary f (Id dadb) (Id dadc) b                  c                  = Reverse (Binary (f (primal b) (primal c)) dadb dadc b c)
+    binary f _         _         Zero     Zero     = Lift (f 0 0)
+    binary f _         _         Zero     (Lift c) = Lift (f 0 c)
+    binary f _         _         (Lift b) Zero     = Lift (f b 0)
+    binary f _         _         (Lift b) (Lift c) = Lift (f b c)
 
+    binary f _         (Id dadc) Zero        (Reverse i c) = unarily (f 0) dadc i c
+    binary f _         (Id dadc) (Lift b)    (Reverse i c) = unarily (f b) dadc i c
+    binary f (Id dadb) _         (Reverse i b) Zero        = unarily (`f` 0) dadb i b
+    binary f (Id dadb) _         (Reverse i b) (Lift c)    = unarily (`f` c) dadb i b
+    binary f (Id dadb) (Id dadc) (Reverse i b) (Reverse j c) = binarily f dadb dadc i b j c
+
     lift2 f df b c = binary f dadb dadc b c
         where (dadb, dadc) = df (Id (primal b)) (Id (primal c))
 
@@ -134,130 +160,68 @@
             a = f pb pc
             (dadb, dadc) = df (Id a) (Id pb) (Id pc)
 
-deriveLifted id (conT ''Reverse)
-
-derivative :: Num a => AD Reverse a -> a
-derivative = sum . map snd . partials
-{-# INLINE derivative #-}
-
-derivative' :: Num a => AD Reverse a -> (a, a)
-derivative' r = (primal r, derivative r)
-{-# INLINE derivative' #-}
-
--- | back propagate sensitivities along a tape.
-backPropagate :: Num a => (Vertex -> (Tape a Int, Int, [Int])) -> STArray s Int a -> Vertex -> ST s ()
-backPropagate vmap ss v = do
-        case node of
-            Unary _ g b -> do
-                da <- readArray ss i
-                db <- readArray ss b
-                writeArray ss b (db + g*da)
-            Binary _ gb gc b c -> do
-                da <- readArray ss i
-                db <- readArray ss b
-                writeArray ss b (db + gb*da)
-                dc <- readArray ss c
-                writeArray ss c (dc + gc*da)
-            _ -> return ()
-    where
-        (node, i, _) = vmap v
-        -- this isn't _quite_ right, as it should allow negative zeros to multiply through
-
-topSortAcyclic :: Graph -> [Vertex]
-topSortAcyclic g = reverse $ runST $ do
-    del <- newArray (bounds g) False :: ST s (STUArray s Int Bool)
-    let tg = transposeG g
-        starters = [ n | (n, []) <- assocs tg ]
-        loop [] rs = return rs
-        loop (n:ns) rs = do
-            writeArray del n True
-            let add [] = return ns
-                add (m:ms) = do
-                    b <- ok (tg!m)
-                    ms' <- add ms
-                    if b then return (m:ms') else return ms'
-                ok [] = return True
-                ok (x:xs) = do b <- readArray del x; if b then ok xs else return False
-            ns' <- add (g!n)
-            loop ns' (n : rs)
-    loop starters []
+let s = varT (mkName "s") in
+  deriveLifted (classP ''Reifies [s, conT ''Tape] :) (conT ''Reverse `appT` s)
 
--- | This returns a list of contributions to the partials.
--- The variable ids returned in the list are likely /not/ unique!
-{-# SPECIALIZE partials :: AD Reverse Double -> [(Int, Double)] #-}
-partials :: forall a . Num a => AD Reverse a -> [(Int, a)]
-partials (AD tape) = [ let v = sensitivities ! ix in seq v (ident, v) | (ix, Var _ ident) <- xs ]
-    where
-        Reified.Graph xs start = unsafePerformIO $ reifyGraph tape
-        g = array xsBounds [ (i, successors t) | (i, t) <- xs ]
-        vertexMap = array xsBounds xs
-        vmap i = (vertexMap ! i, i, [])
-        xsBounds = sbounds xs
+-- | Helper that extracts the derivative of a chain when the chain was constructed with one variable.
+derivativeOf :: (Reifies s Tape, Num a) => Proxy s -> AD (Reverse s) a -> a
+derivativeOf _ = sum . partials
+{-# INLINE derivativeOf #-}
 
-        sensitivities = runSTArray $ do
-            ss <- newArray xsBounds 0
-            writeArray ss start 1
-            forM_ (topSortAcyclic g) $
-                backPropagate vmap ss
-            return ss
+-- | Helper that extracts both the primal and derivative of a chain when the chain was constructed with one variable.
+derivativeOf' :: (Reifies s Tape, Num a) => Proxy s -> AD (Reverse s) a -> (a, a)
+derivativeOf' p r = (primal r, derivativeOf p r)
+{-# INLINE derivativeOf' #-}
 
-        sbounds ((a,_):as) = foldl' (\(lo,hi) (b,_) -> let lo' = min lo b; hi' = max hi b in lo' `seq` hi' `seq` (lo', hi')) (a,a) as
-        sbounds _ = undefined -- the graph can't be empty, it contains the output node!
+-- | Used internally to push sensitivities down the chain.
+backPropagate :: Num a => Int -> Cells -> STArray s Int a -> ST s Int
+backPropagate k Nil _ = return k
+backPropagate k (Unary i g xs) ss = do
+  da <- readArray ss k
+  db <- readArray ss i
+  writeArray ss i $! db + unsafeCoerce g*da
+  (backPropagate $! k - 1) xs ss
+backPropagate k (Binary i j g h xs) ss = do
+  da <- readArray ss k
+  db <- readArray ss i
+  writeArray ss i $! db + unsafeCoerce g*da
+  dc <- readArray ss j
+  writeArray ss j $! dc + unsafeCoerce h*da
+  (backPropagate $! k - 1) xs ss
 
-        successors :: Tape a t -> [t]
-        successors (Unary _ _ b) = [b]
-        successors (Binary _ _ _ b c) = [b,c]
-        successors _ = []
+-- | Extract the partials from the current chain for a given AD variable.
+{-# SPECIALIZE partials :: Reifies s Tape => AD (Reverse s) Double -> [Double] #-}
+partials :: forall s a. (Reifies s Tape, Num a) => AD (Reverse s) a -> [a]
+partials (AD Zero)        = []
+partials (AD (Lift _))    = []
+partials (AD (Reverse k _)) = map (sensitivities !) [0..vs] where
+   Head n t = unsafePerformIO $ readIORef (getTape (reflect (Proxy :: Proxy s)))
+   tk = dropCells (n - k) t
+   (vs,sensitivities) = runST $ do
+     ss <- newArray (0, k) 0
+     writeArray ss k 1
+     v <- backPropagate k tk ss
+     as <- Unsafe.unsafeFreeze ss
+     return (v, as)
 
 -- | Return an 'Array' of 'partials' given bounds for the variable IDs.
-partialArray :: Num a => (Int, Int) -> AD Reverse a -> Array Int a
-partialArray vbounds tape = accumArray (+) 0 vbounds (partials tape)
-{-# INLINE partialArray #-}
+partialArrayOf :: (Reifies s Tape, Num a) => Proxy s -> (Int, Int) -> AD (Reverse s) a -> Array Int a
+partialArrayOf _ vbounds = accumArray (+) 0 vbounds . zip [0..] . partials
+{-# INLINE partialArrayOf #-}
 
 -- | Return an 'IntMap' of sparse partials
-partialMap :: Num a => AD Reverse a -> IntMap a
-partialMap = fromListWith (+) . partials
-{-# INLINE partialMap #-}
+partialMapOf :: (Reifies s Tape, Num a) => Proxy s -> AD (Reverse s) a -> IntMap a
+partialMapOf _ = fromDistinctAscList . zip [0..] . partials
+{-# INLINE partialMapOf #-}
 
--- A simple fresh variable supply monad
-newtype S a = S { runS :: Int -> (a,Int) }
-instance Monad S where
-    return a = S (\s -> (a,s))
-    S g >>= f = S (\s -> let (a,s') = g s in runS (f a) s')
+-- | Construct a tape that starts with @n@ variables.
+reifyTape :: Int -> (forall s. Reifies s Tape => Proxy s -> r) -> r
+reifyTape vs k = unsafePerformIO $ do
+  h <- newIORef (Head vs Nil)
+  return (reify (Tape h) k)
+{-# NOINLINE reifyTape #-}
 
-instance Var Reverse where
-    var a v = Reverse (Var a v)
-    varId (Reverse (Var _ v)) = v
+instance Var (Reverse s) where
+    var a v = Reverse v a
+    varId (Reverse v _) = v
     varId _ = error "varId: not a Var"
-
-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/src/Numeric/AD/Internal/Var.hs b/src/Numeric/AD/Internal/Var.hs
--- a/src/Numeric/AD/Internal/Var.hs
+++ b/src/Numeric/AD/Internal/Var.hs
@@ -8,13 +8,7 @@
 -- Stability   :  experimental
 -- Portability :  GHC only
 --
--- Reverse-Mode Automatic Differentiation using a single tape.
---
--- This version uses @Data.Reflection@ to update a single tape.
---
--- This is asymptotically faster than using @Reverse@, which
--- is forced to reify and topologically sort the graph, but it is
--- less friendly to the use of sparks.
+-- Variables used for reverse-mode automatic differentiation.
 -----------------------------------------------------------------------------
 
 module Numeric.AD.Internal.Var
diff --git a/src/Numeric/AD/Internal/Wengert.hs b/src/Numeric/AD/Internal/Wengert.hs
deleted file mode 100644
--- a/src/Numeric/AD/Internal/Wengert.hs
+++ /dev/null
@@ -1,227 +0,0 @@
-{-# LANGUAGE Rank2Types, TypeFamilies, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, FlexibleContexts, TemplateHaskell, UndecidableInstances, DeriveDataTypeable, GADTs, ScopedTypeVariables #-}
--- {-# OPTIONS_HADDOCK hide, prune #-}
------------------------------------------------------------------------------
--- |
--- Module      :  Numeric.AD.Internal.Wengert
--- Copyright   :  (c) Edward Kmett 2012
--- License     :  BSD3
--- Maintainer  :  ekmett@gmail.com
--- Stability   :  experimental
--- Portability :  GHC only
---
--- Reverse-Mode Automatic Differentiation using a single Wengert list (or \"tape\").
---
--- This version uses @Data.Reflection@ to find and update the tape.
---
--- This is asymptotically faster than using @Reverse@, which
--- is forced to reify and topologically sort the graph, but it requires
--- a fairly expensive rendezvous during construction when updated using
--- multiple threads.
---
------------------------------------------------------------------------------
-
-module Numeric.AD.Internal.Wengert
-    ( Wengert(..)
-    , Tape(..)
-    , Head(..)
-    , Cells(..)
-    , reifyTape
-    , partials
-    , partialArrayOf
-    , partialMapOf
-    , derivativeOf
-    , derivativeOf'
-    ) where
-
-import Control.Monad.ST
-import Data.Array.ST
-import Data.Array
-import Data.Array.Unsafe as Unsafe
-import Data.IORef
-import Data.IntMap (IntMap, fromDistinctAscList)
-import Data.Proxy
-import Data.Reflection
-import Data.Typeable
-import Language.Haskell.TH hiding (reify)
-import Numeric.AD.Internal.Types
-import Numeric.AD.Internal.Classes
-import Numeric.AD.Internal.Identity
-import Numeric.AD.Internal.Var
-import Prelude hiding (mapM)
-import System.IO.Unsafe (unsafePerformIO)
-import Unsafe.Coerce
-
--- evil untyped tape
-data Cells where
-  Nil    :: Cells
-  Unary  :: {-# UNPACK #-} !Int -> a -> Cells -> Cells
-  Binary :: {-# UNPACK #-} !Int -> {-# UNPACK #-} !Int -> a -> a -> Cells -> Cells
-
-dropCells :: Int -> Cells -> Cells
-dropCells 0 xs = xs
-dropCells _ Nil = Nil
-dropCells n (Unary _ _ xs)      = (dropCells $! n - 1) xs
-dropCells n (Binary _ _ _ _ xs) = (dropCells $! n - 1) xs
-
-data Head = Head {-# UNPACK #-} !Int Cells
-
-newtype Tape = Tape { getTape :: IORef Head }
-
-un :: Int -> a -> Head -> (Head, Int)
-un i di (Head r t) = h `seq` r' `seq` (h, r') where
-  r' = r + 1
-  h = Head r' (Unary i di t)
-{-# INLINE un #-}
-
-bin :: Int -> Int -> a -> a -> Head -> (Head, Int)
-bin i j di dj (Head r t) = h `seq` r' `seq` (h, r') where
-  r' = r + 1
-  h = Head r' (Binary i j di dj t)
-{-# INLINE bin #-}
-
-modifyTape :: Reifies s Tape => p s -> (Head -> (Head, r)) -> IO r
-modifyTape p = atomicModifyIORef (getTape (reflect p))
-{-# INLINE modifyTape #-}
-
--- | This is used to create a new entry on the chain given a unary function, its derivative with respect to its input,
--- the variable ID of its input, and the value of its input. Used by 'unary' and 'binary' internally.
-unarily :: forall s a. Reifies s Tape => (a -> a) -> a -> Int -> a -> Wengert s a
-unarily f di i b = Wengert (unsafePerformIO (modifyTape (Proxy :: Proxy s) (un i di))) $! f b
-{-# INLINE unarily #-}
-
--- | This is used to create a new entry on the chain given a binary function, its derivatives with respect to its inputs,
--- their variable IDs and values. Used by 'binary' internally.
-binarily :: forall s a. Reifies s Tape => (a -> a -> a) -> a -> a -> Int -> a -> Int -> a -> Wengert s a
-binarily f di dj i b j c = Wengert (unsafePerformIO (modifyTape (Proxy :: Proxy s) (bin i j di dj))) $! f b c
-{-# INLINE binarily #-}
-
-data Wengert s a where
-  Zero :: Wengert s a
-  Lift :: a -> Wengert s a
-  Wengert :: {-# UNPACK #-} !Int -> a -> Wengert s a
-  deriving (Show, Typeable)
-
-instance (Reifies s Tape, Lifted (Wengert s)) => Mode (Wengert s) where
-  isKnownZero Zero = True
-  isKnownZero _    = False
-
-  isKnownConstant Wengert{} = False
-  isKnownConstant _ = True
-
-  auto = Lift
-  zero = Zero
-  (<+>)  = binary (+) one one
-  a *^ b = lift1 (a *) (\_ -> auto a) b
-  a ^* b = lift1 (* b) (\_ -> auto b) a
-  a ^/ b = lift1 (/ b) (\_ -> auto (recip b)) a
-
-  Zero <**> y      = auto (0 ** primal y)
-  _    <**> Zero   = auto 1
-  x    <**> Lift y = lift1 (**y) (\z -> (y *^ z ** Id (y-1))) x
-  x    <**> y      = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y
-
-instance Primal (Wengert s) where
-    primal Zero = 0
-    primal (Lift a) = a
-    primal (Wengert _ a) = a
-
-instance (Reifies s Tape, Lifted (Wengert s)) => Jacobian (Wengert s) where
-    type D (Wengert s) = Id
-
-    unary f _         (Zero)   = Lift (f 0)
-    unary f _         (Lift a) = Lift (f a)
-    unary f (Id dadi) (Wengert i b) = unarily f dadi i b
-
-    lift1 f df b = unary f (df (Id pb)) b
-        where pb = primal b
-
-    lift1_ f df b = unary (const a) (df (Id a) (Id pb)) b
-        where pb = primal b
-              a = f pb
-
-    binary f _         _         Zero     Zero     = Lift (f 0 0)
-    binary f _         _         Zero     (Lift c) = Lift (f 0 c)
-    binary f _         _         (Lift b) Zero     = Lift (f b 0)
-    binary f _         _         (Lift b) (Lift c) = Lift (f b c)
-
-    binary f _         (Id dadc) Zero        (Wengert i c) = unarily (f 0) dadc i c
-    binary f _         (Id dadc) (Lift b)    (Wengert i c) = unarily (f b) dadc i c
-    binary f (Id dadb) _         (Wengert i b) Zero        = unarily (`f` 0) dadb i b
-    binary f (Id dadb) _         (Wengert i b) (Lift c)    = unarily (`f` c) dadb i b
-    binary f (Id dadb) (Id dadc) (Wengert i b) (Wengert j c) = binarily f dadb dadc i b j c
-
-    lift2 f df b c = binary f dadb dadc b c
-        where (dadb, dadc) = df (Id (primal b)) (Id (primal c))
-
-    lift2_ f df b c = binary (\_ _ -> a) dadb dadc b c
-        where
-            pb = primal b
-            pc = primal c
-            a = f pb pc
-            (dadb, dadc) = df (Id a) (Id pb) (Id pc)
-
-let s = varT (mkName "s") in
-  deriveLifted (classP ''Reifies [s, conT ''Tape] :) (conT ''Wengert `appT` s)
-
--- | Helper that extracts the derivative of a chain when the chain was constructed with one variable.
-derivativeOf :: (Reifies s Tape, Num a) => Proxy s -> AD (Wengert s) a -> a
-derivativeOf _ = sum . partials
-{-# INLINE derivativeOf #-}
-
--- | Helper that extracts both the primal and derivative of a chain when the chain was constructed with one variable.
-derivativeOf' :: (Reifies s Tape, Num a) => Proxy s -> AD (Wengert s) a -> (a, a)
-derivativeOf' p r = (primal r, derivativeOf p r)
-{-# INLINE derivativeOf' #-}
-
--- | Used internally to push sensitivities down the chain.
-backPropagate :: Num a => Int -> Cells -> STArray s Int a -> ST s Int
-backPropagate k Nil _ = return k
-backPropagate k (Unary i g xs) ss = do
-  da <- readArray ss k
-  db <- readArray ss i
-  writeArray ss i $! db + unsafeCoerce g*da
-  (backPropagate $! k - 1) xs ss
-backPropagate k (Binary i j g h xs) ss = do
-  da <- readArray ss k
-  db <- readArray ss i
-  writeArray ss i $! db + unsafeCoerce g*da
-  dc <- readArray ss j
-  writeArray ss j $! dc + unsafeCoerce h*da
-  (backPropagate $! k - 1) xs ss
-
--- | Extract the partials from the current chain for a given AD variable.
-{-# SPECIALIZE partials :: Reifies s Tape => AD (Wengert s) Double -> [Double] #-}
-partials :: forall s a. (Reifies s Tape, Num a) => AD (Wengert s) a -> [a]
-partials (AD Zero)        = []
-partials (AD (Lift _))    = []
-partials (AD (Wengert k _)) = map (sensitivities !) [0..vs] where
-   Head n t = unsafePerformIO $ readIORef (getTape (reflect (Proxy :: Proxy s)))
-   tk = dropCells (n - k) t
-   (vs,sensitivities) = runST $ do
-     ss <- newArray (0, k) 0
-     writeArray ss k 1
-     v <- backPropagate k tk ss
-     as <- Unsafe.unsafeFreeze ss
-     return (v, as)
-
--- | Return an 'Array' of 'partials' given bounds for the variable IDs.
-partialArrayOf :: (Reifies s Tape, Num a) => Proxy s -> (Int, Int) -> AD (Wengert s) a -> Array Int a
-partialArrayOf _ vbounds = accumArray (+) 0 vbounds . zip [0..] . partials
-{-# INLINE partialArrayOf #-}
-
--- | Return an 'IntMap' of sparse partials
-partialMapOf :: (Reifies s Tape, Num a) => Proxy s -> AD (Wengert s) a -> IntMap a
-partialMapOf _ = fromDistinctAscList . zip [0..] . partials
-{-# INLINE partialMapOf #-}
-
--- | Construct a tape that starts with @n@ variables.
-reifyTape :: Int -> (forall s. Reifies s Tape => Proxy s -> r) -> r
-reifyTape vs k = unsafePerformIO $ do
-  h <- newIORef (Head vs Nil)
-  return (reify (Tape h) k)
-{-# NOINLINE reifyTape #-}
-
-instance Var (Wengert s) where
-    var a v = Wengert v a
-    varId (Wengert v _) = v
-    varId _ = error "varId: not a Var"
diff --git a/src/Numeric/AD/Mode/Directed.hs b/src/Numeric/AD/Mode/Directed.hs
--- a/src/Numeric/AD/Mode/Directed.hs
+++ b/src/Numeric/AD/Mode/Directed.hs
@@ -30,67 +30,65 @@
 import Prelude hiding (reverse)
 import Numeric.AD.Types
 import Data.Traversable (Traversable)
-import qualified Numeric.AD.Mode.Reverse as R
+import qualified Numeric.AD.Mode.Kahn as K
 import qualified Numeric.AD.Mode.Forward as F
 import qualified Numeric.AD.Mode.Tower as T
-import qualified Numeric.AD.Mode.Wengert as W
+import qualified Numeric.AD.Mode.Reverse as R
 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
+    | Kahn
     | Reverse
-    | Wengert
     | Tower
     | Mixed
     deriving (Show, Eq, Ord, Read, Bounded, Enum, Ix)
 
 diff :: Num a => Direction -> (forall s. Mode s => AD s a -> AD s a) -> a -> a
 diff Forward = F.diff
+diff Kahn    = K.diff
 diff Reverse = R.diff
-diff Wengert = W.diff
-diff Tower = T.diff
-diff Mixed = F.diff
+diff Tower   = T.diff
+diff Mixed   = F.diff
 {-# INLINE diff #-}
 
 diff' :: Num a => Direction -> (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)
 diff' Forward = F.diff'
+diff' Kahn = K.diff'
 diff' Reverse = R.diff'
-diff' Wengert = W.diff'
 diff' Tower = T.diff'
 diff' Mixed = F.diff'
 {-# INLINE diff' #-}
 
 jacobian :: (Traversable f, Traversable g, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)
 jacobian Forward = F.jacobian
+jacobian Kahn    = K.jacobian
 jacobian Reverse = R.jacobian
-jacobian Wengert = W.jacobian
-jacobian Tower = F.jacobian -- error "jacobian Tower: unimplemented"
-jacobian Mixed = M.jacobian
+jacobian Tower   = F.jacobian -- error "jacobian Tower: unimplemented"
+jacobian Mixed   = M.jacobian
 {-# INLINE jacobian #-}
 
 jacobian' :: (Traversable f, Traversable g, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)
 jacobian' Forward = F.jacobian'
+jacobian' Kahn    = K.jacobian'
 jacobian' Reverse = R.jacobian'
-jacobian' Wengert = W.jacobian'
-jacobian' Tower = F.jacobian' -- error "jacobian' Tower: unimplemented"
-jacobian' Mixed = M.jacobian'
+jacobian' Tower   = F.jacobian' -- error "jacobian' Tower: unimplemented"
+jacobian' Mixed   = M.jacobian'
 {-# INLINE jacobian' #-}
 
 grad :: (Traversable f, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a
 grad Forward = F.grad
+grad Kahn    = K.grad
 grad Reverse = R.grad
-grad Wengert   = W.grad
 grad Tower   = F.grad -- error "grad Tower: unimplemented"
 grad Mixed   = M.grad
 {-# INLINE grad #-}
 
 grad' :: (Traversable f, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)
 grad' Forward = F.grad'
+grad' Kahn    = K.grad'
 grad' Reverse = R.grad'
-grad' Wengert   = W.grad'
 grad' Tower   = F.grad' -- error "grad' Tower: unimplemented"
 grad' Mixed   = M.grad'
 {-# INLINE grad' #-}
diff --git a/src/Numeric/AD/Mode/Kahn.hs b/src/Numeric/AD/Mode/Kahn.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/AD/Mode/Kahn.hs
@@ -0,0 +1,204 @@
+{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, UndecidableInstances, ScopedTypeVariables #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD.Mode.Kahn
+-- Copyright   :  (c) Edward Kmett 2010
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only
+--
+-- This module provides reverse-mode Automatic Differentiation using post-hoc linear time
+-- topological sorting.
+--
+-- 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.Kahn
+    (
+    -- * Gradient
+      grad
+    , grad'
+    , gradWith
+    , gradWith'
+
+    -- * Jacobian
+    , jacobian
+    , jacobian'
+    , jacobianWith
+    , jacobianWith'
+    -- * Hessian
+    , hessian
+    , hessianF
+    -- * Derivatives
+    , diff
+    , diff'
+    , diffF
+    , diffF'
+    -- * Unsafe Variadic Gradient
+    , vgrad, vgrad'
+    , Grad
+    ) where
+
+import Control.Applicative ((<$>))
+import Data.Traversable (Traversable)
+
+import Numeric.AD.Types
+import Numeric.AD.Internal.Classes
+import Numeric.AD.Internal.Composition
+import Numeric.AD.Internal.Kahn
+import Numeric.AD.Internal.Var
+
+-- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with reverse-mode AD in a single pass.
+--
+-- >>> grad (\[x,y,z] -> x*y+z) [1,2,3]
+-- [2,1,1]
+
+grad :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a
+grad f as = 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-mode AD in a single pass.
+--
+-- >>> grad' (\[x,y,z] -> 4*x*exp y+cos z) [1,2,3]
+-- (28.566231899122155,[29.5562243957226,29.5562243957226,-0.1411200080598672])
+grad' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)
+grad' f as = (primal r, unbind vs $ partialArray bds r)
+    where (vs, bds) = bind as
+          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) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f b
+gradWith g f as = unbindWith g vs (partialArray bds $ f vs)
+    where (vs,bds) = bind as
+{-# INLINE gradWith #-}
+
+-- | @'grad'' g f@ calculates the result and gradient of a non-scalar-to-scalar function @f@ with reverse-mode AD in a single pass
+-- the gradient is combined element-wise with the argument using the function @g@.
+--
+-- @'grad'' == 'gradWith'' (\_ dx -> dx)@
+gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f b)
+gradWith' g f as = (primal r, unbindWith g vs $ partialArray bds r)
+    where (vs, bds) = bind as
+          r = f vs
+{-# INLINE gradWith' #-}
+
+-- | The 'jacobian' function calculates the jacobian of a non-scalar-to-non-scalar function with reverse AD lazily in @m@ passes for @m@ outputs.
+--
+-- >>> jacobian (\[x,y] -> [y,x,x*y]) [2,1]
+-- [[0,1],[1,0],[1,2]]
+--
+-- >>> jacobian (\[x,y] -> [exp y,cos x,x+y]) [1,2]
+-- [[0.0,7.38905609893065],[-0.8414709848078965,0.0],[1.0,1.0]]
+jacobian :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)
+jacobian f as = unbind vs . partialArray bds <$> f vs where
+    (vs, bds) = bind as
+{-# INLINE jacobian #-}
+
+-- | The 'jacobian'' function calculates both the result and the Jacobian of a nonscalar-to-nonscalar function, using @m@ invocations of reverse AD,
+-- where @m@ is the output dimensionality. Applying @fmap snd@ to the result will recover the result of 'jacobian'
+-- | An alias for 'gradF''
+--
+-- ghci> jacobian' (\[x,y] -> [y,x,x*y]) [2,1]
+-- [(1,[0,1]),(2,[1,0]),(2,[1,2])]
+jacobian' :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)
+jacobian' f as = row <$> f vs where
+    (vs, bds) = bind as
+    row a = (primal a, unbind vs (partialArray bds a))
+{-# INLINE jacobian' #-}
+
+-- | 'jacobianWith g f' calculates the Jacobian of a non-scalar-to-non-scalar function @f@ with reverse AD lazily in @m@ passes for @m@ outputs.
+--
+-- Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@.
+--
+-- @
+-- 'jacobian' = 'jacobianWith' (\_ dx -> dx)
+-- 'jacobianWith' 'const' = (\f x -> 'const' x '<$>' f x)
+-- @
+jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f b)
+jacobianWith g f as = unbindWith g vs . partialArray bds <$> f vs where
+    (vs, bds) = bind as
+{-# INLINE jacobianWith #-}
+
+-- | 'jacobianWith' g f' calculates both the result and the Jacobian of a nonscalar-to-nonscalar function @f@, using @m@ invocations of reverse AD,
+-- where @m@ is the output dimensionality. Applying @fmap snd@ to the result will recover the result of 'jacobianWith'
+--
+-- Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@.
+--
+-- @'jacobian'' == 'jacobianWith'' (\_ dx -> dx)@
+jacobianWith' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f b)
+jacobianWith' g f as = row <$> f vs where
+    (vs, bds) = bind as
+    row a = (primal a, unbindWith g vs (partialArray bds a))
+{-# INLINE jacobianWith' #-}
+
+-- | Compute the derivative of a function.
+--
+-- >>> diff sin 0
+-- 1.0
+--
+-- >>> cos 0
+-- 1.0
+diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a
+diff f a = derivative $ f (var a 0)
+{-# INLINE diff #-}
+
+-- | The 'diff'' function calculates the value and derivative, as a
+-- pair, of a scalar-to-scalar function.
+--
+--
+-- >>> diff' sin 0
+-- (0.0,1.0)
+diff' :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)
+diff' f a = derivative' $ f (var a 0)
+{-# INLINE diff' #-}
+
+-- | Compute the derivatives of a function that returns a vector with regards to its single input.
+--
+-- >>> diffF (\a -> [sin a, cos a]) 0
+-- [1.0,0.0]
+diffF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a
+diffF f a = derivative <$> f (var a 0)
+{-# INLINE diffF #-}
+
+-- | Compute the derivatives of a function that returns a vector with regards to its single input
+-- as well as the primal answer.
+--
+-- >>> diffF' (\a -> [sin a, cos a]) 0
+-- [(0.0,1.0),(1.0,0.0)]
+diffF' :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f (a, a)
+diffF' f a = derivative' <$> f (var a 0)
+{-# INLINE diffF' #-}
+
+-- | Compute the 'hessian' via the 'jacobian' of the gradient. gradient is computed in reverse mode and then the 'jacobian' is computed in reverse mode.
+--
+-- However, since the @'grad' f :: f a -> f a@ is square this is not as fast as using the forward-mode 'jacobian' of a reverse mode gradient provided by 'Numeric.AD.hessian'.
+--
+-- >>> hessian (\[x,y] -> x*y) [1,2]
+-- [[0,1],[1,0]]
+hessian :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f (f a)
+hessian f = jacobian (grad (decomposeMode . f . fmap composeMode))
+
+-- | Compute the order 3 Hessian tensor on a non-scalar-to-non-scalar function via the reverse-mode Jacobian of the reverse-mode Jacobian of the function.
+--
+-- Less efficient than 'Numeric.AD.Mode.Mixed.hessianF'.
+--
+-- >>> hessianF (\[x,y] -> [x*y,x+y,exp x*cos y]) [1,2]
+-- [[[0.0,1.0],[1.0,0.0]],[[0.0,0.0],[0.0,0.0]],[[-1.1312043837568135,-2.4717266720048188],[-2.4717266720048188,1.1312043837568135]]]
+hessianF :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f (f a))
+hessianF f = decomposeFunctor . jacobian (ComposeFunctor . jacobian (fmap decomposeMode . f . fmap composeMode))
+
diff --git a/src/Numeric/AD/Mode/Reverse.hs b/src/Numeric/AD/Mode/Reverse.hs
--- a/src/Numeric/AD/Mode/Reverse.hs
+++ b/src/Numeric/AD/Mode/Reverse.hs
@@ -8,12 +8,8 @@
 -- 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.
+-- Reverse-mode automatic differentiation using Wengert lists and
+-- Data.Reflection
 --
 -----------------------------------------------------------------------------
 
@@ -30,17 +26,16 @@
     , jacobian'
     , jacobianWith
     , jacobianWith'
+
     -- * Hessian
     , hessian
     , hessianF
+
     -- * Derivatives
     , diff
     , diff'
     , diffF
     , diffF'
-    -- * Unsafe Variadic Gradient
-    , vgrad, vgrad'
-    , Grad
     ) where
 
 import Control.Applicative ((<$>))
@@ -52,72 +47,70 @@
 import Numeric.AD.Internal.Reverse
 import Numeric.AD.Internal.Var
 
--- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with 'Reverse' AD in a single pass.
+-- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with reverse-mode AD in a single pass.
 --
+--
 -- >>> grad (\[x,y,z] -> x*y+z) [1,2,3]
 -- [2,1,1]
-
 grad :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a
-grad f as = unbind vs (partialArray bds $ f vs)
-    where (vs,bds) = bind as
+grad f as = reifyTape (snd bds) $ \p -> unbind vs $! partialArrayOf p bds $! f $ vary <$> vs
+  where (vs, bds) = bind as
 {-# INLINE grad #-}
 
--- | The 'grad'' function calculates the result and gradient of a non-scalar-to-scalar function with 'Reverse' AD in a single pass.
+-- | The 'grad'' function calculates the result and gradient of a non-scalar-to-scalar function with reverse-mode AD in a single pass.
 --
--- >>> grad' (\[x,y,z] -> 4*x*exp y+cos z) [1,2,3]
--- (28.566231899122155,[29.5562243957226,29.5562243957226,-0.1411200080598672])
+-- >>> grad' (\[x,y,z] -> x*y+z) [1,2,3]
+-- (5,[2,1,1])
 grad' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)
-grad' f as = (primal r, unbind vs $ partialArray bds r)
-    where (vs, bds) = bind as
-          r = f vs
+grad' f as = reifyTape (snd bds) $ \p ->
+  let r = f (fmap vary vs) in (primal r, unbind vs $! partialArrayOf p bds $! r)
+  where (vs, bds) = bind as
 {-# INLINE grad' #-}
 
 -- | @'grad' g f@ function calculates the gradient of a non-scalar-to-scalar function @f@ with reverse-mode AD in a single pass.
 -- The gradient is combined element-wise with the argument using the function @g@.
 --
 -- @
--- 'grad' = 'gradWith' (\_ dx -> dx)
--- 'id' = 'gradWith' const
+-- 'grad' == 'gradWith' (\_ dx -> dx)
+-- 'id' == 'gradWith' 'const'
 -- @
---
---
 gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f b
-gradWith g f as = unbindWith g vs (partialArray bds $ f vs)
-    where (vs,bds) = bind as
+gradWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs $! partialArrayOf p bds $! f $ vary <$> vs
+  where (vs,bds) = bind as
 {-# INLINE gradWith #-}
 
--- | @'grad'' g f@ calculates the result and gradient of a non-scalar-to-scalar function @f@ with 'Reverse' AD in a single pass
+-- | @'grad'' g f@ calculates the result and gradient of a non-scalar-to-scalar function @f@ with reverse-mode AD in a single pass
 -- the gradient is combined element-wise with the argument using the function @g@.
 --
--- @'grad'' == 'gradWith'' (\_ dx -> dx)@
+-- @
+-- 'grad'' == 'gradWith'' (\_ dx -> dx)
+-- @
 gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f b)
-gradWith' g f as = (primal r, unbindWith g vs $ partialArray bds r)
+gradWith' g f as = reifyTape (snd bds) $ \p ->
+   let r = f (fmap vary vs) in (primal r, unbindWith g vs $! partialArrayOf p bds $! r)
     where (vs, bds) = bind as
-          r = f vs
 {-# INLINE gradWith' #-}
 
 -- | The 'jacobian' function calculates the jacobian of a non-scalar-to-non-scalar function with reverse AD lazily in @m@ passes for @m@ outputs.
 --
 -- >>> jacobian (\[x,y] -> [y,x,x*y]) [2,1]
 -- [[0,1],[1,0],[1,2]]
---
--- >>> jacobian (\[x,y] -> [exp y,cos x,x+y]) [1,2]
--- [[0.0,7.38905609893065],[-0.8414709848078965,0.0],[1.0,1.0]]
 jacobian :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)
-jacobian f as = unbind vs . partialArray bds <$> f vs where
-    (vs, bds) = bind as
+jacobian f as = reifyTape (snd bds) $ \p -> unbind vs . partialArrayOf p bds <$> f (fmap vary vs)
+  where (vs, bds) = bind as
 {-# INLINE jacobian #-}
 
 -- | The 'jacobian'' function calculates both the result and the Jacobian of a nonscalar-to-nonscalar function, using @m@ invocations of reverse AD,
 -- where @m@ is the output dimensionality. Applying @fmap snd@ to the result will recover the result of 'jacobian'
 -- | An alias for 'gradF''
 --
--- ghci> jacobian' (\[x,y] -> [y,x,x*y]) [2,1]
+-- >>> jacobian' (\[x,y] -> [y,x,x*y]) [2,1]
 -- [(1,[0,1]),(2,[1,0]),(2,[1,2])]
 jacobian' :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)
-jacobian' f as = row <$> f vs where
-    (vs, bds) = bind as
-    row a = (primal a, unbind vs (partialArray bds a))
+jacobian' f as = reifyTape (snd bds) $ \p ->
+  let row a = (primal a, unbind vs $! partialArrayOf p bds $! a)
+  in row <$> f (vary <$> vs)
+  where (vs, bds) = bind as
 {-# INLINE jacobian' #-}
 
 -- | 'jacobianWith g f' calculates the Jacobian of a non-scalar-to-non-scalar function @f@ with reverse AD lazily in @m@ passes for @m@ outputs.
@@ -125,11 +118,11 @@
 -- Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@.
 --
 -- @
--- 'jacobian' = 'jacobianWith' (\_ dx -> dx)
--- 'jacobianWith' 'const' = (\f x -> 'const' x '<$>' f x)
+-- 'jacobian' == 'jacobianWith' (\_ dx -> dx)
+-- 'jacobianWith' 'const' == (\f x -> 'const' x '<$>' f x)
 -- @
 jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f b)
-jacobianWith g f as = unbindWith g vs . partialArray bds <$> f vs where
+jacobianWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs . partialArrayOf p bds <$> f (fmap vary vs) where
     (vs, bds) = bind as
 {-# INLINE jacobianWith #-}
 
@@ -139,53 +132,53 @@
 -- Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@.
 --
 -- @'jacobian'' == 'jacobianWith'' (\_ dx -> dx)@
+--
 jacobianWith' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f b)
-jacobianWith' g f as = row <$> f vs where
-    (vs, bds) = bind as
-    row a = (primal a, unbindWith g vs (partialArray bds a))
+jacobianWith' g f as = reifyTape (snd bds) $ \p ->
+  let row a = (primal a, unbindWith g vs $! partialArrayOf p bds $! a)
+  in row <$> f (vary <$> vs)
+  where (vs, bds) = bind as
 {-# INLINE jacobianWith' #-}
 
 -- | Compute the derivative of a function.
 --
 -- >>> diff sin 0
 -- 1.0
---
--- >>> cos 0
--- 1.0
 diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a
-diff f a = derivative $ f (var a 0)
+diff f a = reifyTape 1 $ \p -> derivativeOf p $! f (var a 0)
 {-# INLINE diff #-}
 
--- | The 'diff'' function calculates the value and derivative, as a
--- pair, of a scalar-to-scalar function.
---
+-- | The 'diff'' function calculates the result and derivative, as a pair, of a scalar-to-scalar function.
 --
 -- >>> diff' sin 0
 -- (0.0,1.0)
+--
+-- >>> diff' exp 0
+-- (1.0,1.0)
 diff' :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)
-diff' f a = derivative' $ f (var a 0)
+diff' f a = reifyTape 1 $ \p -> derivativeOf' p $! f (var a 0)
 {-# INLINE diff' #-}
 
--- | Compute the derivatives of a function that returns a vector with regards to its single input.
+-- | Compute the derivatives of each result of a scalar-to-vector function with regards to its input.
 --
 -- >>> diffF (\a -> [sin a, cos a]) 0
 -- [1.0,0.0]
+--
 diffF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a
-diffF f a = derivative <$> f (var a 0)
+diffF f a = reifyTape 1 $ \p -> derivativeOf p <$> f (var a 0)
 {-# INLINE diffF #-}
 
--- | Compute the derivatives of a function that returns a vector with regards to its single input
--- as well as the primal answer.
+-- | Compute the derivatives of each result of a scalar-to-vector function with regards to its input along with the answer.
 --
 -- >>> diffF' (\a -> [sin a, cos a]) 0
 -- [(0.0,1.0),(1.0,0.0)]
 diffF' :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f (a, a)
-diffF' f a = derivative' <$> f (var a 0)
+diffF' f a = reifyTape 1 $ \p -> derivativeOf' p <$> f (var a 0)
 {-# INLINE diffF' #-}
 
--- | Compute the 'hessian' via the 'jacobian' of the gradient. gradient is computed in reverse mode and then the 'jacobian' is computed in reverse mode.
+-- | Compute the hessian via the jacobian of the gradient. gradient is computed in reverse mode and then the jacobian is computed in reverse mode.
 --
--- However, since the @'grad' f :: f a -> f a@ is square this is not as fast as using the forward-mode 'jacobian' of a reverse mode gradient provided by 'Numeric.AD.hessian'.
+-- However, since the @'grad' f :: f a -> f a@ is square this is not as fast as using the forward-mode Jacobian of a reverse mode gradient provided by 'Numeric.AD.hessian'.
 --
 -- >>> hessian (\[x,y] -> x*y) [1,2]
 -- [[0,1],[1,0]]
@@ -200,4 +193,3 @@
 -- [[[0.0,1.0],[1.0,0.0]],[[0.0,0.0],[0.0,0.0]],[[-1.1312043837568135,-2.4717266720048188],[-2.4717266720048188,1.1312043837568135]]]
 hessianF :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f (f a))
 hessianF f = decomposeFunctor . jacobian (ComposeFunctor . jacobian (fmap decomposeMode . f . fmap composeMode))
-
diff --git a/src/Numeric/AD/Mode/Wengert.hs b/src/Numeric/AD/Mode/Wengert.hs
deleted file mode 100644
--- a/src/Numeric/AD/Mode/Wengert.hs
+++ /dev/null
@@ -1,194 +0,0 @@
-{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, UndecidableInstances, ScopedTypeVariables #-}
------------------------------------------------------------------------------
--- |
--- Module      :  Numeric.AD.Mode.Wengert
--- Copyright   :  (c) Edward Kmett 2010
--- License     :  BSD3
--- Maintainer  :  ekmett@gmail.com
--- Stability   :  experimental
--- Portability :  GHC only
---
--- Reverse-mode automatic differentiation using Wengert lists and Data.Reflection
---
------------------------------------------------------------------------------
-
-module Numeric.AD.Mode.Wengert
-    (
-    -- * Gradient
-      grad
-    , grad'
-    , gradWith
-    , gradWith'
-
-    -- * Jacobian
-    , jacobian
-    , jacobian'
-    , jacobianWith
-    , jacobianWith'
-
-    -- * Hessian
-    , hessian
-    , hessianF
-
-    -- * Derivatives
-    , diff
-    , diff'
-    , diffF
-    , diffF'
-    ) where
-
-import Control.Applicative ((<$>))
-import Data.Traversable (Traversable)
-
-import Numeric.AD.Types
-import Numeric.AD.Internal.Classes
-import Numeric.AD.Internal.Composition
-import Numeric.AD.Internal.Wengert
-import Numeric.AD.Internal.Var
-
--- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with reverse-mode AD in a single pass.
---
---
--- >>> grad (\[x,y,z] -> x*y+z) [1,2,3]
--- [2,1,1]
-grad :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a
-grad f as = reifyTape (snd bds) $ \p -> unbind vs $! partialArrayOf p bds $! f $ vary <$> vs
-  where (vs, bds) = bind as
-{-# INLINE grad #-}
-
--- | The 'grad'' function calculates the result and gradient of a non-scalar-to-scalar function with reverse-mode AD in a single pass.
---
--- >>> grad' (\[x,y,z] -> x*y+z) [1,2,3]
--- (5,[2,1,1])
-grad' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)
-grad' f as = reifyTape (snd bds) $ \p ->
-  let r = f (fmap vary vs) in (primal r, unbind vs $! partialArrayOf p bds $! r)
-  where (vs, bds) = bind as
-{-# INLINE grad' #-}
-
--- | @'grad' g f@ function calculates the gradient of a non-scalar-to-scalar function @f@ with reverse-mode AD in a single pass.
--- The gradient is combined element-wise with the argument using the function @g@.
---
--- @
--- 'grad' == 'gradWith' (\_ dx -> dx)
--- 'id' == 'gradWith' 'const'
--- @
-gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f b
-gradWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs $! partialArrayOf p bds $! f $ vary <$> vs
-  where (vs,bds) = bind as
-{-# INLINE gradWith #-}
-
--- | @'grad'' g f@ calculates the result and gradient of a non-scalar-to-scalar function @f@ with reverse-mode AD in a single pass
--- the gradient is combined element-wise with the argument using the function @g@.
---
--- @
--- 'grad'' == 'gradWith'' (\_ dx -> dx)
--- @
-gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f b)
-gradWith' g f as = reifyTape (snd bds) $ \p ->
-   let r = f (fmap vary vs) in (primal r, unbindWith g vs $! partialArrayOf p bds $! r)
-    where (vs, bds) = bind as
-{-# INLINE gradWith' #-}
-
--- | The 'jacobian' function calculates the jacobian of a non-scalar-to-non-scalar function with reverse AD lazily in @m@ passes for @m@ outputs.
---
--- >>> jacobian (\[x,y] -> [y,x,x*y]) [2,1]
--- [[0,1],[1,0],[1,2]]
-jacobian :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)
-jacobian f as = reifyTape (snd bds) $ \p -> unbind vs . partialArrayOf p bds <$> f (fmap vary vs)
-  where (vs, bds) = bind as
-{-# INLINE jacobian #-}
-
--- | The 'jacobian'' function calculates both the result and the Jacobian of a nonscalar-to-nonscalar function, using @m@ invocations of reverse AD,
--- where @m@ is the output dimensionality. Applying @fmap snd@ to the result will recover the result of 'jacobian'
--- | An alias for 'gradF''
---
--- >>> jacobian' (\[x,y] -> [y,x,x*y]) [2,1]
--- [(1,[0,1]),(2,[1,0]),(2,[1,2])]
-jacobian' :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)
-jacobian' f as = reifyTape (snd bds) $ \p ->
-  let row a = (primal a, unbind vs $! partialArrayOf p bds $! a)
-  in row <$> f (vary <$> vs)
-  where (vs, bds) = bind as
-{-# INLINE jacobian' #-}
-
--- | 'jacobianWith g f' calculates the Jacobian of a non-scalar-to-non-scalar function @f@ with reverse AD lazily in @m@ passes for @m@ outputs.
---
--- Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@.
---
--- @
--- 'jacobian' == 'jacobianWith' (\_ dx -> dx)
--- 'jacobianWith' 'const' == (\f x -> 'const' x '<$>' f x)
--- @
-jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f b)
-jacobianWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs . partialArrayOf p bds <$> f (fmap vary vs) where
-    (vs, bds) = bind as
-{-# INLINE jacobianWith #-}
-
--- | 'jacobianWith' g f' calculates both the result and the Jacobian of a nonscalar-to-nonscalar function @f@, using @m@ invocations of reverse AD,
--- where @m@ is the output dimensionality. Applying @fmap snd@ to the result will recover the result of 'jacobianWith'
---
--- Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@.
---
--- @'jacobian'' == 'jacobianWith'' (\_ dx -> dx)@
---
-jacobianWith' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f b)
-jacobianWith' g f as = reifyTape (snd bds) $ \p ->
-  let row a = (primal a, unbindWith g vs $! partialArrayOf p bds $! a)
-  in row <$> f (vary <$> vs)
-  where (vs, bds) = bind as
-{-# INLINE jacobianWith' #-}
-
--- | Compute the derivative of a function.
---
--- >>> diff sin 0
--- 1.0
-diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a
-diff f a = reifyTape 1 $ \p -> derivativeOf p $! f (var a 0)
-{-# INLINE diff #-}
-
--- | The 'diff'' function calculates the result and derivative, as a pair, of a scalar-to-scalar function.
---
--- >>> diff' sin 0
--- (0.0,1.0)
---
--- >>> diff' exp 0
--- (1.0,1.0)
-diff' :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)
-diff' f a = reifyTape 1 $ \p -> derivativeOf' p $! f (var a 0)
-{-# INLINE diff' #-}
-
--- | Compute the derivatives of each result of a scalar-to-vector function with regards to its input.
---
--- >>> diffF (\a -> [sin a, cos a]) 0
--- [1.0,0.0]
---
-diffF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a
-diffF f a = reifyTape 1 $ \p -> derivativeOf p <$> f (var a 0)
-{-# INLINE diffF #-}
-
--- | Compute the derivatives of each result of a scalar-to-vector function with regards to its input along with the answer.
---
--- >>> diffF' (\a -> [sin a, cos a]) 0
--- [(0.0,1.0),(1.0,0.0)]
-diffF' :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f (a, a)
-diffF' f a = reifyTape 1 $ \p -> derivativeOf' p <$> f (var a 0)
-{-# INLINE diffF' #-}
-
--- | Compute the hessian via the jacobian of the gradient. gradient is computed in reverse mode and then the jacobian is computed in reverse mode.
---
--- However, since the @'grad' f :: f a -> f a@ is square this is not as fast as using the forward-mode Jacobian of a reverse mode gradient provided by 'Numeric.AD.hessian'.
---
--- >>> hessian (\[x,y] -> x*y) [1,2]
--- [[0,1],[1,0]]
-hessian :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f (f a)
-hessian f = jacobian (grad (decomposeMode . f . fmap composeMode))
-
--- | Compute the order 3 Hessian tensor on a non-scalar-to-non-scalar function via the reverse-mode Jacobian of the reverse-mode Jacobian of the function.
---
--- Less efficient than 'Numeric.AD.Mode.Mixed.hessianF'.
---
--- >>> hessianF (\[x,y] -> [x*y,x+y,exp x*cos y]) [1,2]
--- [[[0.0,1.0],[1.0,0.0]],[[0.0,0.0],[0.0,0.0]],[[-1.1312043837568135,-2.4717266720048188],[-2.4717266720048188,1.1312043837568135]]]
-hessianF :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f (f a))
-hessianF f = decomposeFunctor . jacobian (ComposeFunctor . jacobian (fmap decomposeMode . f . fmap composeMode))
diff --git a/src/Numeric/AD/Variadic.hs b/src/Numeric/AD/Variadic.hs
--- a/src/Numeric/AD/Variadic.hs
+++ b/src/Numeric/AD/Variadic.hs
@@ -25,5 +25,5 @@
     , Grads, vgrads
     ) where
 
-import Numeric.AD.Variadic.Reverse
+import Numeric.AD.Variadic.Kahn
 import Numeric.AD.Variadic.Sparse (Grads, vgrads)
diff --git a/src/Numeric/AD/Variadic/Kahn.hs b/src/Numeric/AD/Variadic/Kahn.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/AD/Variadic/Kahn.hs
@@ -0,0 +1,27 @@
+{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, UndecidableInstances, ScopedTypeVariables #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD.Variadic.Kahn
+-- Copyright   :  (c) Edward Kmett 2010-2012
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  non-portable
+--
+-- Variadic combinators for reverse-mode automatic differentiation.
+--
+-- Unfortunately, variadicity comes at the expense of being able to use
+-- quantification to avoid sensitivity confusion, so be careful when
+-- counting the number of @lift@ you use when taking the gradient of a
+-- function that takes gradients!
+--
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Variadic.Kahn
+    (
+    -- * Unsafe Variadic Gradient
+      vgrad, vgrad'
+    , Grad
+    ) where
+
+import Numeric.AD.Internal.Kahn
diff --git a/src/Numeric/AD/Variadic/Reverse.hs b/src/Numeric/AD/Variadic/Reverse.hs
deleted file mode 100644
--- a/src/Numeric/AD/Variadic/Reverse.hs
+++ /dev/null
@@ -1,27 +0,0 @@
-{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, UndecidableInstances, ScopedTypeVariables #-}
------------------------------------------------------------------------------
--- |
--- Module      :  Numeric.AD.Variadic.Reverse
--- Copyright   :  (c) Edward Kmett 2010-2012
--- License     :  BSD3
--- Maintainer  :  ekmett@gmail.com
--- Stability   :  experimental
--- Portability :  non-portable
---
--- Variadic combinators for reverse-mode automatic differentiation.
---
--- Unfortunately, variadicity comes at the expense of being able to use
--- quantification to avoid sensitivity confusion, so be careful when
--- counting the number of @lift@ you use when taking the gradient of a
--- function that takes gradients!
---
------------------------------------------------------------------------------
-
-module Numeric.AD.Variadic.Reverse
-    (
-    -- * Unsafe Variadic Gradient
-      vgrad, vgrad'
-    , Grad
-    ) where
-
-import Numeric.AD.Internal.Reverse