diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown
--- a/CHANGELOG.markdown
+++ b/CHANGELOG.markdown
@@ -1,3 +1,9 @@
+3.2
+---
+* Renamed `Chain` to `Wengert` to reflect its use of Wengert lists for reverse mode.
+* Renamed `lift` to `auto` to avoid conflict with the more prevalent `transformers` library.
+* Fixed a bug in `Numeric.AD.Forward.gradWith'`, which caused it to return the wrong value for the primal.
+
 3.1.4
 -----
 * Added a better "convergence" test for `findZero`
diff --git a/ad.cabal b/ad.cabal
--- a/ad.cabal
+++ b/ad.cabal
@@ -1,5 +1,5 @@
 name:         ad
-version:      3.1.4
+version:      3.2
 license:      BSD3
 license-File: LICENSE
 copyright:    (c) Edward Kmett 2010-2012,
@@ -26,7 +26,7 @@
     .
     * @Numeric.AD.Mode.Reverse@ uses benign side-effects to compute reverse-mode AD. It is good for computing gradients in one pass. It generates a tree-like tape that needs to be topologically sorted in the end.
     .
-    * @Numeric.AD.Mode.Chain@ uses benign side-effects to compute reverse-mode AD. It is good for computing gradients in one pass. It generates a linear tape using @Data.Reflection@.
+    * @Numeric.AD.Mode.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.Sparse@ computes a sparse forward-mode AD tower. It is good for higher derivatives or large numbers of outputs.
     .
@@ -104,7 +104,7 @@
     Numeric.AD.Newton
     Numeric.AD.Halley
 
-    Numeric.AD.Mode.Chain
+    Numeric.AD.Mode.Wengert
     Numeric.AD.Mode.Directed
     Numeric.AD.Mode.Forward
     Numeric.AD.Mode.Reverse
@@ -121,7 +121,7 @@
     Numeric.AD.Internal.Tower
     Numeric.AD.Internal.Reverse
     Numeric.AD.Internal.Var
-    Numeric.AD.Internal.Chain
+    Numeric.AD.Internal.Wengert
     Numeric.AD.Internal.Sparse
     Numeric.AD.Internal.Dense
     Numeric.AD.Internal.Composition
diff --git a/src/Numeric/AD/Halley.hs b/src/Numeric/AD/Halley.hs
--- a/src/Numeric/AD/Halley.hs
+++ b/src/Numeric/AD/Halley.hs
@@ -58,7 +58,7 @@
 -- Note: the @take 10 $ inverse sqrt 1 (sqrt 10)@ example that works for Newton's method
 -- fails with Halley's method because the preconditions do not hold!
 inverse :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> a -> [a]
-inverse f x0 y = findZero (\x -> f x - lift y) x0
+inverse f x0 y = findZero (\x -> f x - auto y) x0
 {-# INLINE inverse  #-}
 
 -- | The 'fixedPoint' function find a fixedpoint of a scalar
diff --git a/src/Numeric/AD/Internal/Chain.hs b/src/Numeric/AD/Internal/Chain.hs
deleted file mode 100644
--- a/src/Numeric/AD/Internal/Chain.hs
+++ /dev/null
@@ -1,226 +0,0 @@
-{-# LANGUAGE Rank2Types, TypeFamilies, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, FlexibleContexts, TemplateHaskell, UndecidableInstances, DeriveDataTypeable, GADTs, ScopedTypeVariables #-}
--- {-# OPTIONS_HADDOCK hide, prune #-}
------------------------------------------------------------------------------
--- |
--- Module      :  Numeric.AD.Internal.Chain
--- Copyright   :  (c) Edward Kmett 2012
--- License     :  BSD3
--- Maintainer  :  ekmett@gmail.com
--- Stability   :  experimental
--- Portability :  GHC only
---
--- Reverse-Mode Automatic Differentiation using a single tape.
---
--- This version uses @Data.Reflection@ to find and update the tape
---
--- This is asymptotically faster than using @Reverse@, which
--- is forced to reify and topologically sort the graph, but it requires
--- a fairly expensive rendezvous during construction.
---
------------------------------------------------------------------------------
-
-module Numeric.AD.Internal.Chain
-    ( Chain(..)
-    , Tape(..)
-    , Head(..)
-    , Cells(..)
-    , reifyTape
-    , partials
-    , partialArrayOf
-    , partialMapOf
-    , derivativeOf
-    , derivativeOf'
-    ) where
-
-import Control.Monad.ST
-import Data.Array.ST
-import Data.Array
-import Data.Array.Unsafe as Unsafe
-import Data.IORef
-import Data.IntMap (IntMap, fromDistinctAscList)
-import Data.Proxy
-import Data.Reflection
-import Data.Typeable
-import Language.Haskell.TH hiding (reify)
-import Numeric.AD.Internal.Types
-import Numeric.AD.Internal.Classes
-import Numeric.AD.Internal.Identity
-import Numeric.AD.Internal.Var
-import Prelude hiding (mapM)
-import System.IO.Unsafe (unsafePerformIO)
-import Unsafe.Coerce
-
--- evil untyped tape
-data Cells where
-  Nil    :: Cells
-  Unary  :: {-# UNPACK #-} !Int -> a -> Cells -> Cells
-  Binary :: {-# UNPACK #-} !Int -> {-# UNPACK #-} !Int -> a -> a -> Cells -> Cells
-
-dropCells :: Int -> Cells -> Cells
-dropCells 0 xs = xs
-dropCells _ Nil = Nil
-dropCells n (Unary _ _ xs)      = (dropCells $! n - 1) xs
-dropCells n (Binary _ _ _ _ xs) = (dropCells $! n - 1) xs
-
-data Head = Head {-# UNPACK #-} !Int Cells
-
-newtype Tape = Tape { getTape :: IORef Head }
-
-un :: Int -> a -> Head -> (Head, Int)
-un i di (Head r t) = h `seq` r' `seq` (h, r') where
-  r' = r + 1
-  h = Head r' (Unary i di t)
-{-# INLINE un #-}
-
-bin :: Int -> Int -> a -> a -> Head -> (Head, Int)
-bin i j di dj (Head r t) = h `seq` r' `seq` (h, r') where
-  r' = r + 1
-  h = Head r' (Binary i j di dj t)
-{-# INLINE bin #-}
-
-modifyTape :: Reifies s Tape => p s -> (Head -> (Head, r)) -> IO r
-modifyTape p = atomicModifyIORef (getTape (reflect p))
-{-# INLINE modifyTape #-}
-
--- | This is used to create a new entry on the chain given a unary function, its derivative with respect to its input,
--- the variable ID of its input, and the value of its input. Used by 'unary' and 'binary' internally.
-unarily :: forall s a. Reifies s Tape => (a -> a) -> a -> Int -> a -> Chain s a
-unarily f di i b = Chain (unsafePerformIO (modifyTape (Proxy :: Proxy s) (un i di))) $! f b
-{-# INLINE unarily #-}
-
--- | This is used to create a new entry on the chain given a binary function, its derivatives with respect to its inputs,
--- their variable IDs and values. Used by 'binary' internally.
-binarily :: forall s a. Reifies s Tape => (a -> a -> a) -> a -> a -> Int -> a -> Int -> a -> Chain s a
-binarily f di dj i b j c = Chain (unsafePerformIO (modifyTape (Proxy :: Proxy s) (bin i j di dj))) $! f b c
-{-# INLINE binarily #-}
-
-data Chain s a where
-  Zero :: Chain s a
-  Lift :: a -> Chain s a
-  Chain :: {-# UNPACK #-} !Int -> a -> Chain s a
-  deriving (Show, Typeable)
-
-instance (Reifies s Tape, Lifted (Chain s)) => Mode (Chain s) where
-  isKnownZero Zero = True
-  isKnownZero _    = False
-
-  isKnownConstant Chain{} = False
-  isKnownConstant _ = True
-
-  lift = Lift
-  zero = Zero
-  (<+>)  = binary (+) one one
-  a *^ b = lift1 (a *) (\_ -> lift a) b
-  a ^* b = lift1 (* b) (\_ -> lift b) a
-  a ^/ b = lift1 (/ b) (\_ -> lift (recip b)) a
-
-  Zero <**> y      = lift (0 ** primal y)
-  _    <**> Zero   = lift 1
-  x    <**> Lift y = lift1 (**y) (\z -> (y *^ z ** Id (y-1))) x
-  x    <**> y      = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y
-
-instance Primal (Chain s) where
-    primal Zero = 0
-    primal (Lift a) = a
-    primal (Chain _ a) = a
-
-instance (Reifies s Tape, Lifted (Chain s)) => Jacobian (Chain s) where
-    type D (Chain s) = Id
-
-    unary f _         (Zero)   = Lift (f 0)
-    unary f _         (Lift a) = Lift (f a)
-    unary f (Id dadi) (Chain i b) = unarily f dadi i b
-
-    lift1 f df b = unary f (df (Id pb)) b
-        where pb = primal b
-
-    lift1_ f df b = unary (const a) (df (Id a) (Id pb)) b
-        where pb = primal b
-              a = f pb
-
-    binary f _         _         Zero     Zero     = Lift (f 0 0)
-    binary f _         _         Zero     (Lift c) = Lift (f 0 c)
-    binary f _         _         (Lift b) Zero     = Lift (f b 0)
-    binary f _         _         (Lift b) (Lift c) = Lift (f b c)
-
-    binary f _         (Id dadc) Zero        (Chain i c) = unarily (f 0) dadc i c
-    binary f _         (Id dadc) (Lift b)    (Chain i c) = unarily (f b) dadc i c
-    binary f (Id dadb) _         (Chain i b) Zero        = unarily (`f` 0) dadb i b
-    binary f (Id dadb) _         (Chain i b) (Lift c)    = unarily (`f` c) dadb i b
-    binary f (Id dadb) (Id dadc) (Chain i b) (Chain j c) = binarily f dadb dadc i b j c
-
-    lift2 f df b c = binary f dadb dadc b c
-        where (dadb, dadc) = df (Id (primal b)) (Id (primal c))
-
-    lift2_ f df b c = binary (\_ _ -> a) dadb dadc b c
-        where
-            pb = primal b
-            pc = primal c
-            a = f pb pc
-            (dadb, dadc) = df (Id a) (Id pb) (Id pc)
-
-let s = varT (mkName "s") in
-  deriveLifted (classP ''Reifies [s, conT ''Tape] :) (conT ''Chain `appT` s)
-
--- | Helper that extracts the derivative of a chain when the chain was constructed with one variable.
-derivativeOf :: (Reifies s Tape, Num a) => Proxy s -> AD (Chain s) a -> a
-derivativeOf _ = sum . partials
-{-# INLINE derivativeOf #-}
-
--- | Helper that extracts both the primal and derivative of a chain when the chain was constructed with one variable.
-derivativeOf' :: (Reifies s Tape, Num a) => Proxy s -> AD (Chain s) a -> (a, a)
-derivativeOf' p r = (primal r, derivativeOf p r)
-{-# INLINE derivativeOf' #-}
-
--- | Used internally to push sensitivities down the chain.
-backPropagate :: Num a => Int -> Cells -> STArray s Int a -> ST s Int
-backPropagate k Nil _ = return k
-backPropagate k (Unary i g xs) ss = do
-  da <- readArray ss k
-  db <- readArray ss i
-  writeArray ss i $! db + unsafeCoerce g*da
-  (backPropagate $! k - 1) xs ss
-backPropagate k (Binary i j g h xs) ss = do
-  da <- readArray ss k
-  db <- readArray ss i
-  writeArray ss i $! db + unsafeCoerce g*da
-  dc <- readArray ss j
-  writeArray ss j $! dc + unsafeCoerce h*da
-  (backPropagate $! k - 1) xs ss
-
--- | Extract the partials from the current chain for a given AD variable.
-{-# SPECIALIZE partials :: Reifies s Tape => AD (Chain s) Double -> [Double] #-}
-partials :: forall s a. (Reifies s Tape, Num a) => AD (Chain s) a -> [a]
-partials (AD Zero)        = []
-partials (AD (Lift _))    = []
-partials (AD (Chain k _)) = map (sensitivities !) [0..vs] where
-   Head n t = unsafePerformIO $ readIORef (getTape (reflect (Proxy :: Proxy s)))
-   tk = dropCells (n - k) t
-   (vs,sensitivities) = runST $ do
-     ss <- newArray (0, k) 0
-     writeArray ss k 1
-     v <- backPropagate k tk ss
-     as <- Unsafe.unsafeFreeze ss
-     return (v, as)
-
--- | Return an 'Array' of 'partials' given bounds for the variable IDs.
-partialArrayOf :: (Reifies s Tape, Num a) => Proxy s -> (Int, Int) -> AD (Chain s) a -> Array Int a
-partialArrayOf _ vbounds = accumArray (+) 0 vbounds . zip [0..] . partials
-{-# INLINE partialArrayOf #-}
-
--- | Return an 'IntMap' of sparse partials
-partialMapOf :: (Reifies s Tape, Num a) => Proxy s -> AD (Chain s) a -> IntMap a
-partialMapOf _ = fromDistinctAscList . zip [0..] . partials
-{-# INLINE partialMapOf #-}
-
--- | Construct a tape that starts with @n@ variables.
-reifyTape :: Int -> (forall s. Reifies s Tape => Proxy s -> r) -> r
-reifyTape vs k = unsafePerformIO $ do
-  h <- newIORef (Head vs Nil)
-  return (reify (Tape h) k)
-{-# NOINLINE reifyTape #-}
-
-instance Var (Chain s) where
-    var a v = Chain v a
-    varId (Chain v _) = v
-    varId _ = error "varId: not a Var"
diff --git a/src/Numeric/AD/Internal/Classes.hs b/src/Numeric/AD/Internal/Classes.hs
--- a/src/Numeric/AD/Internal/Classes.hs
+++ b/src/Numeric/AD/Internal/Classes.hs
@@ -92,7 +92,7 @@
     isKnownZero _ = False
 
     -- | Embed a constant
-    lift  :: Num a => a -> t a
+    auto  :: Num a => a -> t a
 
     -- | Vector sum
     (<+>) :: Num a => t a -> t a -> t a
@@ -113,19 +113,19 @@
     -- | > 'zero' = 'lift' 0
     zero :: Num a => t a
 
-    a *^ b = lift a *! b
-    a ^* b = a *! lift b
+    a *^ b = auto a *! b
+    a ^* b = a *! auto b
 
     a ^/ b = a ^* recip b
 
-    zero = lift 0
+    zero = auto 0
 
 one :: (Mode t, Num a) => t a
-one = lift 1
+one = auto 1
 {-# INLINE one #-}
 
 negOne :: (Mode t, Num a) => t a
-negOne = lift (-1)
+negOne = auto (-1)
 {-# INLINE negOne #-}
 
 -- | 'Primal' is used by 'deriveMode' but is not exposed
@@ -133,7 +133,7 @@
 -- via the AD data type.
 --
 -- It provides direct access to the result, stripped of its derivative information,
--- but this is unsafe in general as (lift . primal) would discard derivative
+-- but this is unsafe in general as (auto . primal) would discard derivative
 -- information. The end user is protected from accidentally using this function
 -- by the universal quantification on the various combinators we expose.
 
@@ -200,11 +200,11 @@
        instance Lifted $_t where
         (==!)         = (==) `on` primal
         compare1      = compare `on` primal
-        maxBound1     = lift maxBound
-        minBound1     = lift minBound
+        maxBound1     = auto maxBound
+        minBound1     = auto minBound
         showsPrec1 d  = showsPrec d . primal
         fromInteger1 0 = zero
-        fromInteger1 n = lift (fromInteger n)
+        fromInteger1 n = auto (fromInteger n)
         (+!)          = (<+>) -- binary (+) one one
         (-!)          = binary (-) one negOne -- TODO: <-> ? as it is, this might be pretty bad for Tower
         (*!)          = lift2 (*) (\x y -> (y, x))
@@ -212,14 +212,14 @@
         abs1          = lift1 abs signum1
         signum1       = lift1 signum (const zero)
         fromRational1 0 = zero
-        fromRational1 r = lift (fromRational r)
+        fromRational1 r = auto (fromRational r)
         x /! y        = x *! recip1 y
         recip1        = lift1_ recip (const . negate1 . square1)
-        pi1       = lift pi
+        pi1       = auto pi
         exp1      = lift1_ exp const
         log1      = lift1 log recip1
         logBase1 x y = log1 y /! log1 x
-        sqrt1     = lift1_ sqrt (\z _ -> recip1 (lift 2 *! z))
+        sqrt1     = lift1_ sqrt (\z _ -> recip1 (auto 2 *! z))
         (**!)     = (<**>)
         --x **! y
         --   | isKnownZero y     = 1
@@ -240,7 +240,7 @@
 
         succ1                 = lift1 succ (const one)
         pred1                 = lift1 pred (const one)
-        toEnum1               = lift . toEnum
+        toEnum1               = auto . toEnum
         fromEnum1             = discrete1 fromEnum
         enumFrom1 a           = withPrimal a <$> discrete1 enumFrom a
         enumFromTo1 a b       = withPrimal a <$> discrete2 enumFromTo a b
@@ -252,7 +252,7 @@
         floatDigits1     = discrete1 floatDigits
         floatRange1      = discrete1 floatRange
         decodeFloat1     = discrete1 decodeFloat
-        encodeFloat1 m e = lift (encodeFloat m e)
+        encodeFloat1 m e = auto (encodeFloat m e)
         isNaN1           = discrete1 isNaN
         isInfinite1      = discrete1 isInfinite
         isDenormalized1  = discrete1 isDenormalized
diff --git a/src/Numeric/AD/Internal/Composition.hs b/src/Numeric/AD/Internal/Composition.hs
--- a/src/Numeric/AD/Internal/Composition.hs
+++ b/src/Numeric/AD/Internal/Composition.hs
@@ -93,18 +93,18 @@
     primal = primal . primal . runComposeMode
 
 instance (Mode f, Mode g) => Mode (ComposeMode f g) where
-    lift = ComposeMode . lift . lift
+    auto = ComposeMode . auto . auto
     ComposeMode a <+> ComposeMode b = ComposeMode (a <+> b)
-    a *^ ComposeMode b = ComposeMode (lift a *^ b)
-    ComposeMode a ^* b = ComposeMode (a ^* lift b)
-    ComposeMode a ^/ b = ComposeMode (a ^/ lift b)
+    a *^ ComposeMode b = ComposeMode (auto a *^ b)
+    ComposeMode a ^* b = ComposeMode (a ^* auto b)
+    ComposeMode a ^/ b = ComposeMode (a ^/ auto b)
     ComposeMode a <**> ComposeMode b = ComposeMode (a <**> b)
 
 instance (Mode f, Mode g) => Lifted (ComposeMode f g) where
     showsPrec1 n (ComposeMode a) = showsPrec1 n a
     ComposeMode a ==! ComposeMode b  = a ==! b
     compare1 (ComposeMode a) (ComposeMode b) = compare1 a b
-    fromInteger1 = ComposeMode . lift . fromInteger1
+    fromInteger1 = ComposeMode . auto . fromInteger1
     ComposeMode a +! ComposeMode b = ComposeMode (a +! b)
     ComposeMode a -! ComposeMode b = ComposeMode (a -! b)
     ComposeMode a *! ComposeMode b = ComposeMode (a *! b)
@@ -113,7 +113,7 @@
     signum1 (ComposeMode a) = ComposeMode (signum1 a)
     ComposeMode a /! ComposeMode b = ComposeMode (a /! b)
     recip1 (ComposeMode a) = ComposeMode (recip1 a)
-    fromRational1 = ComposeMode . lift . fromRational1
+    fromRational1 = ComposeMode . auto . fromRational1
     toRational1 (ComposeMode a) = toRational1 a
     pi1 = ComposeMode pi1
     exp1 (ComposeMode a) = ComposeMode (exp1 a)
diff --git a/src/Numeric/AD/Internal/Dense.hs b/src/Numeric/AD/Internal/Dense.hs
--- a/src/Numeric/AD/Internal/Dense.hs
+++ b/src/Numeric/AD/Internal/Dense.hs
@@ -80,7 +80,7 @@
     primal (Dense a _) = a
 
 instance (Traversable f, Lifted (Dense f)) => Mode (Dense f) where
-    lift = Lift
+    auto = Lift
     zero = Zero
 
     Zero <+> a = a
@@ -90,8 +90,8 @@
     Dense a da <+> Lift b     = Dense (a + b) da
     Dense a da <+> Dense b db = Dense (a + b) $ zipWithT (+) da db
 
-    Zero <**> y      = lift (0 ** primal y)
-    _    <**> Zero   = lift 1
+    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
 
diff --git a/src/Numeric/AD/Internal/Forward.hs b/src/Numeric/AD/Internal/Forward.hs
--- a/src/Numeric/AD/Internal/Forward.hs
+++ b/src/Numeric/AD/Internal/Forward.hs
@@ -37,7 +37,7 @@
 import Numeric.AD.Internal.Classes
 import Numeric.AD.Internal.Identity
 
--- | 'Forward' mode AD.
+-- | 'Forward' mode AD
 data Forward a
   = Forward !a a
   | Lift !a
@@ -70,7 +70,7 @@
     primal Zero = 0
 
 instance Lifted Forward => Mode Forward where
-    lift = Lift
+    auto = Lift
     zero = Zero
 
     isKnownZero Zero = True
@@ -86,8 +86,8 @@
     Lift a <+> Forward b db = Forward (a + b) db
     Lift a <+> Lift b = Lift (a + b)
 
-    Zero <**> y      = lift (0 ** primal y)
-    _    <**> Zero   = lift 1
+    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
 
@@ -168,28 +168,28 @@
 bind f as = snd $ mapAccumL outer (0 :: Int) as
     where
         outer !i _ = (i + 1, f $ snd $ mapAccumL (inner i) 0 as)
-        inner !i !j a = (j + 1, if i == j then bundle a 1 else lift a)
+        inner !i !j a = (j + 1, if i == j then bundle a 1 else auto a)
 
 bind' :: (Traversable f, Num a) => (f (AD Forward a) -> b) -> f a -> (b, f b)
 bind' f as = dropIx $ mapAccumL outer (0 :: Int, b0) as
     where
         outer (!i, _) _ = let b = f $ snd $ mapAccumL (inner i) (0 :: Int) as in ((i + 1, b), b)
-        inner !i !j a = (j + 1, if i == j then bundle a 1 else lift a)
-        b0 = f (lift <$> as)
+        inner !i !j a = (j + 1, if i == j then bundle a 1 else auto a)
+        b0 = f (auto <$> as)
         dropIx ((_,b),bs) = (b,bs)
 
 bindWith :: (Traversable f, Num a) => (a -> b -> c) -> (f (AD Forward a) -> b) -> f a -> f c
 bindWith g f as = snd $ mapAccumL outer (0 :: Int) as
     where
         outer !i a = (i + 1, g a $ f $ snd $ mapAccumL (inner i) 0 as)
-        inner !i !j a = (j + 1, if i == j then bundle a 1 else AD Zero)
+        inner !i !j a = (j + 1, if i == j then bundle a 1 else auto a)
 
 bindWith' :: (Traversable f, Num a) => (a -> b -> c) -> (f (AD Forward a) -> b) -> f a -> (b, f c)
 bindWith' g f as = dropIx $ mapAccumL outer (0 :: Int, b0) as
     where
         outer (!i, _) a = let b = f $ snd $ mapAccumL (inner i) (0 :: Int) as in ((i + 1, b), g a b)
-        inner !i !j a = (j + 1, if i == j then bundle a 1 else AD Zero)
-        b0 = f (lift <$> as)
+        inner !i !j a = (j + 1, if i == j then bundle a 1 else auto a)
+        b0 = f (auto <$> as)
         dropIx ((_,b),bs) = (b,bs)
 
 -- we can't transpose arbitrary traversables, since we can't construct one out of whole cloth, and the outer
diff --git a/src/Numeric/AD/Internal/Identity.hs b/src/Numeric/AD/Internal/Identity.hs
--- a/src/Numeric/AD/Internal/Identity.hs
+++ b/src/Numeric/AD/Internal/Identity.hs
@@ -129,7 +129,7 @@
     maxBound1 = maxBound
 
 instance Mode Id where
-    lift = Id
+    auto = Id
     Id a ^* b = Id (a * b)
     a *^ Id b = Id (a * b)
     Id a <+> Id b = Id (a + b)
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
@@ -81,15 +81,15 @@
     isKnownConstant (Reverse (Lift _)) = True
     isKnownConstant _ = False
 
-    lift a = Reverse (Lift a)
+    auto a = Reverse (Lift a)
     zero   = Reverse Zero
     (<+>)  = binary (+) one one
-    a *^ b = lift1 (a *) (\_ -> lift a) b
-    a ^* b = lift1 (* b) (\_ -> lift b) a
-    a ^/ b = lift1 (/ b) (\_ -> lift (recip b)) a
+    a *^ b = lift1 (a *) (\_ -> auto a) b
+    a ^* b = lift1 (* b) (\_ -> auto b) a
+    a ^/ b = lift1 (/ b) (\_ -> auto (recip b)) a
 
-    Reverse Zero <**> y                = lift (0 ** primal y)
-    _            <**> Reverse Zero     = lift 1
+    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
 
diff --git a/src/Numeric/AD/Internal/Sparse.hs b/src/Numeric/AD/Internal/Sparse.hs
--- a/src/Numeric/AD/Internal/Sparse.hs
+++ b/src/Numeric/AD/Internal/Sparse.hs
@@ -72,7 +72,7 @@
 vars :: (Traversable f, Num a) => f a -> f (AD Sparse a)
 vars = snd . mapAccumL var 0
     where
-        var !n a = (n + 1, AD $ Sparse a $ singleton n $ lift 1)
+        var !n a = (n + 1, AD $ Sparse a $ singleton n $ auto 1)
 {-# INLINE vars #-}
 
 apply :: (Traversable f, Num a) => (f (AD Sparse a) -> b) -> f a -> b
@@ -105,7 +105,7 @@
 
 {-
 vvars :: Num a => Vector a -> Vector (AD Sparse a)
-vvars = Vector.imap (\n a -> AD $ Sparse a $ singleton n $ lift 1)
+vvars = Vector.imap (\n a -> AD $ Sparse a $ singleton n $ auto 1)
 {-# INLINE vvars #-}
 
 vapply :: Num a => (Vector (AD Sparse a) -> b) -> Vector a -> b
@@ -131,7 +131,7 @@
 
 partial :: Num a => [Int] -> Sparse a -> a
 partial []     (Sparse a _)  = a
-partial (n:ns) (Sparse _ da) = partial ns $ findWithDefault (lift 0) n da
+partial (n:ns) (Sparse _ da) = partial ns $ findWithDefault (auto 0) n da
 partial _      Zero          = 0
 {-# INLINE partial #-}
 
@@ -148,10 +148,10 @@
     primal Zero = 0
 
 instance Lifted Sparse => Mode Sparse where
-    lift a = Sparse a IntMap.empty
+    auto a = Sparse a IntMap.empty
     zero = Zero
-    Zero <**> y    = lift (0 ** primal y)
-    _    <**> Zero = lift 1
+    Zero <**> y    = auto (0 ** primal y)
+    _    <**> Zero = auto 1
     x    <**> y@(Sparse b bs)
       | IntMap.null bs = lift1 (**b) (\z -> (b *^ z <**> Sparse (b-1) IntMap.empty)) x
       | otherwise      = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y
@@ -167,17 +167,17 @@
 
 instance Lifted Sparse => Jacobian Sparse where
     type D Sparse = Sparse
-    unary f _ Zero = lift (f 0)
+    unary f _ Zero = auto (f 0)
     unary f dadb (Sparse pb bs) = Sparse (f pb) $ mapWithKey (times dadb) bs
 
-    lift1 f _ Zero = lift (f 0)
+    lift1 f _ Zero = auto (f 0)
     lift1 f df b@(Sparse pb bs) = Sparse (f pb) $ mapWithKey (times (df b)) bs
 
-    lift1_ f _  Zero = lift (f 0)
+    lift1_ f _  Zero = auto (f 0)
     lift1_ f df b@(Sparse pb bs) = a where
         a = Sparse (f pb) $ mapWithKey (times (df a b)) bs
 
-    binary f _    _    Zero           Zero           = lift (f 0 0)
+    binary f _    _    Zero           Zero           = auto (f 0 0)
     binary f _    dadc Zero           (Sparse pc dc) = Sparse (f 0  pc) $ mapWithKey (times dadc) dc
     binary f dadb _    (Sparse pb db) Zero           = Sparse (f pb 0 ) $ mapWithKey (times dadb) db
     binary f dadb dadc (Sparse pb db) (Sparse pc dc) = Sparse (f pb pc) $
@@ -185,7 +185,7 @@
             (mapWithKey (times dadb) db)
             (mapWithKey (times dadc) dc)
 
-    lift2 f _  Zero             Zero = lift (f 0 0)
+    lift2 f _  Zero             Zero = auto (f 0 0)
     lift2 f df Zero c@(Sparse pc dc) = Sparse (f 0 pc) $ mapWithKey (times dadc) dc where dadc = snd (df zero c)
     lift2 f df b@(Sparse pb db) Zero = Sparse (f pb 0) $ mapWithKey (times dadb) db where dadb = fst (df b zero)
     lift2 f df b@(Sparse pb db) c@(Sparse pc dc) = Sparse (f pb pc) da where
@@ -194,7 +194,7 @@
             (mapWithKey (times dadb) db)
             (mapWithKey (times dadc) dc)
 
-    lift2_ f _  Zero             Zero = lift (f 0 0)
+    lift2_ f _  Zero             Zero = auto (f 0 0)
     lift2_ f df b@(Sparse pb db) Zero = a where a = Sparse (f pb 0) (mapWithKey (times (fst (df a b zero))) db)
     lift2_ f df Zero c@(Sparse pc dc) = a where a = Sparse (f 0 pc) (mapWithKey (times (snd (df a zero c))) dc)
     lift2_ f df b@(Sparse pb db) c@(Sparse pc dc) = a where
diff --git a/src/Numeric/AD/Internal/Tower.hs b/src/Numeric/AD/Internal/Tower.hs
--- a/src/Numeric/AD/Internal/Tower.hs
+++ b/src/Numeric/AD/Internal/Tower.hs
@@ -103,10 +103,10 @@
     primal _ = 0
 
 instance Lifted Tower => Mode Tower where
-    lift a = Tower [a]
+    auto a = Tower [a]
     zero = Tower []
-    Tower [] <**> y         = lift (0 ** primal y)
-    _        <**> Tower []  = lift 1
+    Tower [] <**> y         = auto (0 ** primal y)
+    _        <**> Tower []  = auto 1
     x        <**> Tower [y] = lift1 (**y) (\z -> (y *^ z <**> Tower [y-1])) x
     x        <**> y         = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y
 
diff --git a/src/Numeric/AD/Internal/Wengert.hs b/src/Numeric/AD/Internal/Wengert.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/AD/Internal/Wengert.hs
@@ -0,0 +1,227 @@
+{-# 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/Chain.hs b/src/Numeric/AD/Mode/Chain.hs
deleted file mode 100644
--- a/src/Numeric/AD/Mode/Chain.hs
+++ /dev/null
@@ -1,194 +0,0 @@
-{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, UndecidableInstances, ScopedTypeVariables #-}
------------------------------------------------------------------------------
--- |
--- Module      :  Numeric.AD.Mode.Chain
--- Copyright   :  (c) Edward Kmett 2010
--- License     :  BSD3
--- Maintainer  :  ekmett@gmail.com
--- Stability   :  experimental
--- Portability :  GHC only
---
--- Reverse Automatic Differentiation using Data.Reflection
---
------------------------------------------------------------------------------
-
-module Numeric.AD.Mode.Chain
-    (
-    -- * Gradient
-      grad
-    , grad'
-    , gradWith
-    , gradWith'
-
-    -- * Jacobian
-    , jacobian
-    , jacobian'
-    , jacobianWith
-    , jacobianWith'
-
-    -- * Hessian
-    , hessian
-    , hessianF
-
-    -- * Derivatives
-    , diff
-    , diff'
-    , diffF
-    , diffF'
-    ) where
-
-import Control.Applicative ((<$>))
-import Data.Traversable (Traversable)
-
-import Numeric.AD.Types
-import Numeric.AD.Internal.Classes
-import Numeric.AD.Internal.Composition
-import Numeric.AD.Internal.Chain
-import Numeric.AD.Internal.Var
-
--- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with reverse-mode AD in a single pass.
---
---
--- >>> grad (\[x,y,z] -> x*y+z) [1,2,3]
--- [2,1,1]
-grad :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a
-grad f as = reifyTape (snd bds) $ \p -> unbind vs $! partialArrayOf p bds $! f $ vary <$> vs
-  where (vs, bds) = bind as
-{-# INLINE grad #-}
-
--- | The 'grad'' function calculates the result and gradient of a non-scalar-to-scalar function with reverse-mode AD in a single pass.
---
--- >>> grad' (\[x,y,z] -> x*y+z) [1,2,3]
--- (5,[2,1,1])
-grad' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)
-grad' f as = reifyTape (snd bds) $ \p ->
-  let r = f (fmap vary vs) in (primal r, unbind vs $! partialArrayOf p bds $! r)
-  where (vs, bds) = bind as
-{-# INLINE grad' #-}
-
--- | @'grad' g f@ function calculates the gradient of a non-scalar-to-scalar function @f@ with reverse-mode AD in a single pass.
--- The gradient is combined element-wise with the argument using the function @g@.
---
--- @
--- 'grad' == 'gradWith' (\_ dx -> dx)
--- 'id' == 'gradWith' 'const'
--- @
-gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f b
-gradWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs $! partialArrayOf p bds $! f $ vary <$> vs
-  where (vs,bds) = bind as
-{-# INLINE gradWith #-}
-
--- | @'grad'' g f@ calculates the result and gradient of a non-scalar-to-scalar function @f@ with reverse-mode AD in a single pass
--- the gradient is combined element-wise with the argument using the function @g@.
---
--- @
--- 'grad'' == 'gradWith'' (\_ dx -> dx)
--- @
-gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f b)
-gradWith' g f as = reifyTape (snd bds) $ \p ->
-   let r = f (fmap vary vs) in (primal r, unbindWith g vs $! partialArrayOf p bds $! r)
-    where (vs, bds) = bind as
-{-# INLINE gradWith' #-}
-
--- | The 'jacobian' function calculates the jacobian of a non-scalar-to-non-scalar function with reverse AD lazily in @m@ passes for @m@ outputs.
---
--- >>> jacobian (\[x,y] -> [y,x,x*y]) [2,1]
--- [[0,1],[1,0],[1,2]]
-jacobian :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)
-jacobian f as = reifyTape (snd bds) $ \p -> unbind vs . partialArrayOf p bds <$> f (fmap vary vs)
-  where (vs, bds) = bind as
-{-# INLINE jacobian #-}
-
--- | The 'jacobian'' function calculates both the result and the Jacobian of a nonscalar-to-nonscalar function, using @m@ invocations of reverse AD,
--- where @m@ is the output dimensionality. Applying @fmap snd@ to the result will recover the result of 'jacobian'
--- | An alias for 'gradF''
---
--- >>> jacobian' (\[x,y] -> [y,x,x*y]) [2,1]
--- [(1,[0,1]),(2,[1,0]),(2,[1,2])]
-jacobian' :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)
-jacobian' f as = reifyTape (snd bds) $ \p ->
-  let row a = (primal a, unbind vs $! partialArrayOf p bds $! a)
-  in row <$> f (vary <$> vs)
-  where (vs, bds) = bind as
-{-# INLINE jacobian' #-}
-
--- | 'jacobianWith g f' calculates the Jacobian of a non-scalar-to-non-scalar function @f@ with reverse AD lazily in @m@ passes for @m@ outputs.
---
--- Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@.
---
--- @
--- 'jacobian' == 'jacobianWith' (\_ dx -> dx)
--- 'jacobianWith' 'const' == (\f x -> 'const' x '<$>' f x)
--- @
-jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f b)
-jacobianWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs . partialArrayOf p bds <$> f (fmap vary vs) where
-    (vs, bds) = bind as
-{-# INLINE jacobianWith #-}
-
--- | 'jacobianWith' g f' calculates both the result and the Jacobian of a nonscalar-to-nonscalar function @f@, using @m@ invocations of reverse AD,
--- where @m@ is the output dimensionality. Applying @fmap snd@ to the result will recover the result of 'jacobianWith'
---
--- Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@.
---
--- @'jacobian'' == 'jacobianWith'' (\_ dx -> dx)@
---
-jacobianWith' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f b)
-jacobianWith' g f as = reifyTape (snd bds) $ \p ->
-  let row a = (primal a, unbindWith g vs $! partialArrayOf p bds $! a)
-  in row <$> f (vary <$> vs)
-  where (vs, bds) = bind as
-{-# INLINE jacobianWith' #-}
-
--- | Compute the derivative of a function.
---
--- >>> diff sin 0
--- 1.0
-diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a
-diff f a = reifyTape 1 $ \p -> derivativeOf p $! f (var a 0)
-{-# INLINE diff #-}
-
--- | The 'diff'' function calculates the result and derivative, as a pair, of a scalar-to-scalar function.
---
--- >>> diff' sin 0
--- (0.0,1.0)
---
--- >>> diff' exp 0
--- (1.0,1.0)
-diff' :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)
-diff' f a = reifyTape 1 $ \p -> derivativeOf' p $! f (var a 0)
-{-# INLINE diff' #-}
-
--- | Compute the derivatives of each result of a scalar-to-vector function with regards to its input.
---
--- >>> diffF (\a -> [sin a, cos a]) 0
--- [1.0,0.0]
---
-diffF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a
-diffF f a = reifyTape 1 $ \p -> derivativeOf p <$> f (var a 0)
-{-# INLINE diffF #-}
-
--- | Compute the derivatives of each result of a scalar-to-vector function with regards to its input along with the answer.
---
--- >>> diffF' (\a -> [sin a, cos a]) 0
--- [(0.0,1.0),(1.0,0.0)]
-diffF' :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f (a, a)
-diffF' f a = reifyTape 1 $ \p -> derivativeOf' p <$> f (var a 0)
-{-# INLINE diffF' #-}
-
--- | Compute the hessian via the jacobian of the gradient. gradient is computed in reverse mode and then the jacobian is computed in reverse mode.
---
--- However, since the @'grad' f :: f a -> f a@ is square this is not as fast as using the forward-mode Jacobian of a reverse mode gradient provided by 'Numeric.AD.hessian'.
---
--- >>> hessian (\[x,y] -> x*y) [1,2]
--- [[0,1],[1,0]]
-hessian :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f (f a)
-hessian f = jacobian (grad (decomposeMode . f . fmap composeMode))
-
--- | Compute the order 3 Hessian tensor on a non-scalar-to-non-scalar function via the reverse-mode Jacobian of the reverse-mode Jacobian of the function.
---
--- Less efficient than 'Numeric.AD.Mode.Mixed.hessianF'.
---
--- >>> hessianF (\[x,y] -> [x*y,x+y,exp x*cos y]) [1,2]
--- [[[0.0,1.0],[1.0,0.0]],[[0.0,0.0],[0.0,0.0]],[[-1.1312043837568135,-2.4717266720048188],[-2.4717266720048188,1.1312043837568135]]]
-hessianF :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f (f a))
-hessianF f = decomposeFunctor . jacobian (ComposeFunctor . jacobian (fmap decomposeMode . f . fmap composeMode))
diff --git a/src/Numeric/AD/Mode/Directed.hs b/src/Numeric/AD/Mode/Directed.hs
--- a/src/Numeric/AD/Mode/Directed.hs
+++ b/src/Numeric/AD/Mode/Directed.hs
@@ -33,7 +33,7 @@
 import qualified Numeric.AD.Mode.Reverse as R
 import qualified Numeric.AD.Mode.Forward as F
 import qualified Numeric.AD.Mode.Tower as T
-import qualified Numeric.AD.Mode.Chain as C
+import qualified Numeric.AD.Mode.Wengert as W
 import qualified Numeric.AD as M
 import Data.Ix
 
@@ -42,7 +42,7 @@
 data Direction
     = Forward
     | Reverse
-    | Chain
+    | Wengert
     | Tower
     | Mixed
     deriving (Show, Eq, Ord, Read, Bounded, Enum, Ix)
@@ -50,7 +50,7 @@
 diff :: Num a => Direction -> (forall s. Mode s => AD s a -> AD s a) -> a -> a
 diff Forward = F.diff
 diff Reverse = R.diff
-diff Chain = C.diff
+diff Wengert = W.diff
 diff Tower = T.diff
 diff Mixed = F.diff
 {-# INLINE diff #-}
@@ -58,7 +58,7 @@
 diff' :: Num a => Direction -> (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)
 diff' Forward = F.diff'
 diff' Reverse = R.diff'
-diff' Chain = C.diff'
+diff' Wengert = W.diff'
 diff' Tower = T.diff'
 diff' Mixed = F.diff'
 {-# INLINE diff' #-}
@@ -66,7 +66,7 @@
 jacobian :: (Traversable f, Traversable g, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)
 jacobian Forward = F.jacobian
 jacobian Reverse = R.jacobian
-jacobian Chain = C.jacobian
+jacobian Wengert = W.jacobian
 jacobian Tower = F.jacobian -- error "jacobian Tower: unimplemented"
 jacobian Mixed = M.jacobian
 {-# INLINE jacobian #-}
@@ -74,7 +74,7 @@
 jacobian' :: (Traversable f, Traversable g, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)
 jacobian' Forward = F.jacobian'
 jacobian' Reverse = R.jacobian'
-jacobian' Chain = C.jacobian'
+jacobian' Wengert = W.jacobian'
 jacobian' Tower = F.jacobian' -- error "jacobian' Tower: unimplemented"
 jacobian' Mixed = M.jacobian'
 {-# INLINE jacobian' #-}
@@ -82,7 +82,7 @@
 grad :: (Traversable f, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a
 grad Forward = F.grad
 grad Reverse = R.grad
-grad Chain   = C.grad
+grad Wengert   = W.grad
 grad Tower   = F.grad -- error "grad Tower: unimplemented"
 grad Mixed   = M.grad
 {-# INLINE grad #-}
@@ -90,7 +90,7 @@
 grad' :: (Traversable f, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)
 grad' Forward = F.grad'
 grad' Reverse = R.grad'
-grad' Chain   = C.grad'
+grad' Wengert   = W.grad'
 grad' Tower   = F.grad' -- error "grad' Tower: unimplemented"
 grad' Mixed   = M.grad'
 {-# INLINE grad' #-}
diff --git a/src/Numeric/AD/Mode/Forward.hs b/src/Numeric/AD/Mode/Forward.hs
--- a/src/Numeric/AD/Mode/Forward.hs
+++ b/src/Numeric/AD/Mode/Forward.hs
@@ -159,14 +159,14 @@
 
 -- | Compute the gradient of a function using forward mode AD.
 --
--- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Chain.grad' for @n@ inputs, in exchange for better space utilization.
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.grad' for @n@ inputs, in exchange for better space utilization.
 grad :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a
 grad f = bind (tangent . f)
 {-# INLINE grad #-}
 
 -- | Compute the gradient and answer to a function using forward mode AD.
 --
--- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Chain.grad'' for @n@ inputs, in exchange for better space utilization.
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.grad'' for @n@ inputs, in exchange for better space utilization.
 grad' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)
 grad' f as = (primal b, tangent <$> bs)
     where
@@ -175,7 +175,7 @@
 
 -- | Compute the gradient of a function using forward mode AD and combine the result with the input using a user-specified function.
 --
--- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Chain.gradWith' for @n@ inputs, in exchange for better space utilization.
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.gradWith' for @n@ inputs, in exchange for better space utilization.
 gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f b
 gradWith g f = bindWith g (tangent . f)
 {-# INLINE gradWith #-}
@@ -183,9 +183,11 @@
 -- | Compute the gradient of a function using forward mode AD and the answer, and combine the result with the input using a
 -- user-specified function.
 --
--- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Chain.gradWith'' for @n@ inputs, in exchange for better space utilization.
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.gradWith'' for @n@ inputs, in exchange for better space utilization.
+-- >>> gradWith' (,) sum [0..4]
+-- (10,[(0,1),(1,1),(2,1),(3,1),(4,1)])
 gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f b)
-gradWith' g f = bindWith' g (tangent . f)
+gradWith' g f as = (primal $ f (AD . Lift <$> as), bindWith g (tangent . f) as)
 {-# INLINE gradWith' #-}
 
 -- | Compute the product of a vector with the Hessian using forward-on-forward-mode AD.
diff --git a/src/Numeric/AD/Mode/Wengert.hs b/src/Numeric/AD/Mode/Wengert.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/AD/Mode/Wengert.hs
@@ -0,0 +1,194 @@
+{-# 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/Newton.hs b/src/Numeric/AD/Newton.hs
--- a/src/Numeric/AD/Newton.hs
+++ b/src/Numeric/AD/Newton.hs
@@ -64,7 +64,7 @@
 -- >>> last $ take 10 $ inverse sqrt 1 (sqrt 10)
 -- 10.0
 inverse :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> a -> [a]
-inverse f x0 y = findZero (\x -> f x - lift y) x0
+inverse f x0 y = findZero (\x -> f x - auto y) x0
 {-# INLINE inverse  #-}
 
 -- | The 'fixedPoint' function find a fixedpoint of a scalar
@@ -129,7 +129,7 @@
     d0 = negate <$> grad f x0
     go xi ri di = xi : go xi1 ri1 di1
       where
-        ai  = last $ take 20 $ extremum (\a -> f $ zipWithT (\x d -> lift x + a * lift d) xi di) 0
+        ai  = last $ take 20 $ extremum (\a -> f $ zipWithT (\x d -> auto x + a * auto d) xi di) 0
         xi1 = zipWithT (\x d -> x + ai*d) xi di
         ri1 = negate <$> grad f xi1
         bi1 = max 0 $ dot ri1 (zipWithT (-) ri1 ri) / dot ri1 ri1