diff --git a/.travis.yml b/.travis.yml
--- a/.travis.yml
+++ b/.travis.yml
@@ -2,11 +2,17 @@
 #
 #   haskell-ci '--output=.travis.yml' '--config=cabal.haskell-ci' 'cabal.project'
 #
+# To regenerate the script (for example after adjusting tested-with) run
+#
+#   haskell-ci regenerate
+#
 # For more information, see https://github.com/haskell-CI/haskell-ci
 #
-# version: 0.5.20180830
+# version: 0.10
 #
+version: ~> 1.0
 language: c
+os: linux
 dist: xenial
 git:
   # whether to recursively clone submodules
@@ -17,11 +23,12 @@
       - irc.freenode.org#haskell-lens
     skip_join: true
     template:
-      - "\"\\x0313ad\\x03/\\x0306%{branch}\\x03 \\x0314%{commit}\\x03 %{build_url} %{message}\""
+      - "\x0313ad\x03/\x0306%{branch}\x03 \x0314%{commit}\x03 %{build_url} %{message}"
 cache:
   directories:
     - $HOME/.cabal/packages
     - $HOME/.cabal/store
+    - $HOME/.hlint
 before_cache:
   - rm -fv $CABALHOME/packages/hackage.haskell.org/build-reports.log
   # remove files that are regenerated by 'cabal update'
@@ -31,30 +38,38 @@
   - rm -fv $CABALHOME/packages/hackage.haskell.org/01-index.tar
   - rm -fv $CABALHOME/packages/hackage.haskell.org/01-index.tar.idx
   - rm -rfv $CABALHOME/packages/head.hackage
-matrix:
+jobs:
   include:
-    - compiler: ghc-8.8.1
-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-8.8.1","cabal-install-3.0"]}}
+    - compiler: ghc-8.10.1
+      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-8.10.1","cabal-install-3.2"]}}
+      os: linux
+    - compiler: ghc-8.8.3
+      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-8.8.3","cabal-install-3.2"]}}
+      os: linux
     - compiler: ghc-8.6.5
-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-8.6.5","cabal-install-3.0"]}}
+      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-8.6.5","cabal-install-3.2"]}}
+      os: linux
     - compiler: ghc-8.4.4
-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-8.4.4","cabal-install-3.0"]}}
+      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-8.4.4","cabal-install-3.2"]}}
+      os: linux
     - compiler: ghc-8.2.2
-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-8.2.2","cabal-install-3.0"]}}
+      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-8.2.2","cabal-install-3.2"]}}
+      os: linux
     - compiler: ghc-8.0.2
-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-8.0.2","cabal-install-3.0"]}}
+      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-8.0.2","cabal-install-3.2"]}}
+      os: linux
     - compiler: ghc-7.10.3
-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-7.10.3","cabal-install-3.0"]}}
+      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-7.10.3","cabal-install-3.2"]}}
+      os: linux
     - compiler: ghc-7.8.4
-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-7.8.4","cabal-install-3.0"]}}
+      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-7.8.4","cabal-install-3.2"]}}
+      os: linux
     - compiler: ghc-7.6.3
-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-7.6.3","cabal-install-3.0"]}}
+      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-7.6.3","cabal-install-3.2"]}}
+      os: linux
     - compiler: ghc-7.4.2
-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-7.4.2","cabal-install-3.0"]}}
-    - compiler: ghc-head
-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-head","cabal-install-head"]}}
-  allow_failures:
-    - compiler: ghc-head
+      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-7.4.2","cabal-install-3.2"]}}
+      os: linux
 before_install:
   - HC=$(echo "/opt/$CC/bin/ghc" | sed 's/-/\//')
   - WITHCOMPILER="-w $HC"
@@ -67,33 +82,11 @@
   - TOP=$(pwd)
   - "HCNUMVER=$(${HC} --numeric-version|perl -ne '/^(\\d+)\\.(\\d+)\\.(\\d+)(\\.(\\d+))?$/; print(10000 * $1 + 100 * $2 + ($3 == 0 ? $5 != 1 : $3))')"
   - echo $HCNUMVER
-  - CABAL="$CABAL -vnormal+nowrap+markoutput"
+  - CABAL="$CABAL -vnormal+nowrap"
   - set -o pipefail
-  - |
-    echo 'function blue(s) { printf "\033[0;34m" s "\033[0m " }'           >> .colorful.awk
-    echo 'BEGIN { state = "output"; }'                                     >> .colorful.awk
-    echo '/^-----BEGIN CABAL OUTPUT-----$/ { state = "cabal" }'            >> .colorful.awk
-    echo '/^-----END CABAL OUTPUT-----$/ { state = "output" }'             >> .colorful.awk
-    echo '!/^(-----BEGIN CABAL OUTPUT-----|-----END CABAL OUTPUT-----)/ {' >> .colorful.awk
-    echo '  if (state == "cabal") {'                                       >> .colorful.awk
-    echo '    print blue($0)'                                              >> .colorful.awk
-    echo '  } else {'                                                      >> .colorful.awk
-    echo '    print $0'                                                    >> .colorful.awk
-    echo '  }'                                                             >> .colorful.awk
-    echo '}'                                                               >> .colorful.awk
-  - cat .colorful.awk
-  - |
-    color_cabal_output () {
-      awk -f $TOP/.colorful.awk
-    }
-  - echo text | color_cabal_output
-install:
-  - ${CABAL} --version
-  - echo "$(${HC} --version) [$(${HC} --print-project-git-commit-id 2> /dev/null || echo '?')]"
   - TEST=--enable-tests
   - BENCH=--enable-benchmarks
   - HEADHACKAGE=false
-  - if [ $HCNUMVER -gt 80801 ] ; then HEADHACKAGE=true ; fi
   - rm -f $CABALHOME/config
   - |
     echo "verbose: normal +nowrap +markoutput"          >> $CABALHOME/config
@@ -111,17 +104,12 @@
     echo "  prefix: $CABALHOME"                         >> $CABALHOME/config
     echo "repository hackage.haskell.org"               >> $CABALHOME/config
     echo "  url: http://hackage.haskell.org/"           >> $CABALHOME/config
+install:
+  - ${CABAL} --version
+  - echo "$(${HC} --version) [$(${HC} --print-project-git-commit-id 2> /dev/null || echo '?')]"
   - |
-    if $HEADHACKAGE; then
-    echo "allow-newer: $($HCPKG list --simple-output | sed -E 's/([a-zA-Z-]+)-[0-9.]+/*:\1/g')" >> $CABALHOME/config
-    echo "repository head.hackage.ghc.haskell.org"                                        >> $CABALHOME/config
-    echo "   url: https://ghc.gitlab.haskell.org/head.hackage/"                           >> $CABALHOME/config
-    echo "   secure: True"                                                                >> $CABALHOME/config
-    echo "   root-keys: 7541f32a4ccca4f97aea3b22f5e593ba2c0267546016b992dfadcd2fe944e55d" >> $CABALHOME/config
-    echo "              26021a13b401500c8eb2761ca95c61f2d625bfef951b939a8124ed12ecf07329" >> $CABALHOME/config
-    echo "              f76d08be13e9a61a377a85e2fb63f4c5435d40f8feb3e12eb05905edb8cdea89" >> $CABALHOME/config
-    echo "   key-threshold: 3"                                                            >> $CABALHOME/config
-    fi
+    echo "program-default-options"                >> $CABALHOME/config
+    echo "  ghc-options: $GHCJOBS +RTS -M6G -RTS" >> $CABALHOME/config
   - cat $CABALHOME/config
   - rm -fv cabal.project cabal.project.local cabal.project.freeze
   - travis_retry ${CABAL} v2-update -v
@@ -130,19 +118,21 @@
   - touch cabal.project
   - |
     echo "packages: ." >> cabal.project
+  - if [ $HCNUMVER -ge 80200 ] ; then echo 'package ad' >> cabal.project ; fi
+  - "if [ $HCNUMVER -ge 80200 ] ; then echo '  ghc-options: -Werror=missing-methods' >> cabal.project ; fi"
   - |
   - "for pkg in $($HCPKG list --simple-output); do echo $pkg | sed 's/-[^-]*$//' | (grep -vE -- '^(ad)$' || true) | sed 's/^/constraints: /' | sed 's/$/ installed/' >> cabal.project.local; done"
   - cat cabal.project || true
   - cat cabal.project.local || true
   - if [ -f "./configure.ac" ]; then (cd "." && autoreconf -i); fi
-  - ${CABAL} v2-freeze $WITHCOMPILER ${TEST} ${BENCH} | color_cabal_output
+  - ${CABAL} v2-freeze $WITHCOMPILER ${TEST} ${BENCH}
   - "cat cabal.project.freeze | sed -E 's/^(constraints: *| *)//' | sed 's/any.//'"
   - rm  cabal.project.freeze
-  - ${CABAL} v2-build $WITHCOMPILER ${TEST} ${BENCH} --dep -j2 all | color_cabal_output
+  - travis_wait 40 ${CABAL} v2-build $WITHCOMPILER ${TEST} ${BENCH} --dep -j2 all
 script:
   - DISTDIR=$(mktemp -d /tmp/dist-test.XXXX)
   # Packaging...
-  - ${CABAL} v2-sdist all | color_cabal_output
+  - ${CABAL} v2-sdist all
   # Unpacking...
   - mv dist-newstyle/sdist/*.tar.gz ${DISTDIR}/
   - cd ${DISTDIR} || false
@@ -154,19 +144,21 @@
   - touch cabal.project
   - |
     echo "packages: ${PKGDIR_ad}" >> cabal.project
+  - if [ $HCNUMVER -ge 80200 ] ; then echo 'package ad' >> cabal.project ; fi
+  - "if [ $HCNUMVER -ge 80200 ] ; then echo '  ghc-options: -Werror=missing-methods' >> cabal.project ; fi"
   - |
   - "for pkg in $($HCPKG list --simple-output); do echo $pkg | sed 's/-[^-]*$//' | (grep -vE -- '^(ad)$' || true) | sed 's/^/constraints: /' | sed 's/$/ installed/' >> cabal.project.local; done"
   - cat cabal.project || true
   - cat cabal.project.local || true
   # Building with tests and benchmarks...
   # build & run tests, build benchmarks
-  - ${CABAL} v2-build $WITHCOMPILER ${TEST} ${BENCH} all | color_cabal_output
+  - ${CABAL} v2-build $WITHCOMPILER ${TEST} ${BENCH} all
   # Testing...
-  - ${CABAL} v2-test $WITHCOMPILER ${TEST} ${BENCH} all | color_cabal_output
+  - ${CABAL} v2-test $WITHCOMPILER ${TEST} ${BENCH} all
   # cabal check...
   - (cd ${PKGDIR_ad} && ${CABAL} -vnormal check)
   # haddock...
-  - ${CABAL} v2-haddock $WITHCOMPILER --with-haddock $HADDOCK ${TEST} ${BENCH} all | color_cabal_output
+  - ${CABAL} v2-haddock $WITHCOMPILER --with-haddock $HADDOCK ${TEST} ${BENCH} all
 
-# REGENDATA ["--output=.travis.yml","--config=cabal.haskell-ci","cabal.project"]
+# REGENDATA ("0.10",["--output=.travis.yml","--config=cabal.haskell-ci","cabal.project"])
 # EOF
diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown
--- a/CHANGELOG.markdown
+++ b/CHANGELOG.markdown
@@ -1,3 +1,14 @@
+4.4.1 [2020.10.13]
+------------------
+* Change the fixity of `:-` in `Numeric.AD.Jet` to be right-associative.
+  Previously, it was `infixl`, which made things like `x :- y :- z` nearly
+  unusable.
+* Fix backpropagation error in Kahn mode.
+* Fix bugs in the `Erf` instance for `ForwardDouble`.
+* Add `Numeric.AD.Mode.Reverse.Double`, a variant of `Numeric.AD.Mode.Reverse`
+  that is specialized to `Double`.
+* Re-export `Jet(..)`, `headJet`, `tailJet` and `jet` from `Numeric.AD`.
+
 4.4 [2020.02.03]
 ----------------
 * Generalize the type of `stochasticGradientDescent`:
diff --git a/ad.cabal b/ad.cabal
--- a/ad.cabal
+++ b/ad.cabal
@@ -1,5 +1,5 @@
 name:          ad
-version:       4.4
+version:       4.4.1
 license:       BSD3
 license-File:  LICENSE
 copyright:     (c) Edward Kmett 2010-2015,
@@ -20,7 +20,8 @@
              , GHC == 8.2.2
              , GHC == 8.4.4
              , GHC == 8.6.5
-             , GHC == 8.8.1
+             , GHC == 8.8.3
+             , GHC == 8.10.1
 synopsis:      Automatic Differentiation
 extra-source-files:
   .gitignore
@@ -148,6 +149,7 @@
     Numeric.AD.Internal.On
     Numeric.AD.Internal.Or
     Numeric.AD.Internal.Reverse
+    Numeric.AD.Internal.Reverse.Double
     Numeric.AD.Internal.Sparse
     Numeric.AD.Internal.Tower
     Numeric.AD.Internal.Type
@@ -158,6 +160,7 @@
     Numeric.AD.Mode.Forward.Double
     Numeric.AD.Mode.Kahn
     Numeric.AD.Mode.Reverse
+    Numeric.AD.Mode.Reverse.Double
     Numeric.AD.Mode.Sparse
     Numeric.AD.Mode.Tower
     Numeric.AD.Newton
@@ -176,6 +179,8 @@
     Numeric.AD.Internal.Combinators
 
   ghc-options: -Wall
+  if impl(ghc >= 8.6)
+    ghc-options: -Wno-star-is-type
 
   ghc-options: -fspec-constr -fdicts-cheap -O2
 
@@ -188,7 +193,7 @@
     ad,
     base,
     directory,
-    doctest >= 0.9.0.1 && < 0.17,
+    doctest >= 0.9.0.1 && < 0.18,
     filepath
   ghc-options: -Wall -threaded
   hs-source-dirs: tests
diff --git a/src/Numeric/AD.hs b/src/Numeric/AD.hs
--- a/src/Numeric/AD.hs
+++ b/src/Numeric/AD.hs
@@ -129,6 +129,12 @@
   , conjugateGradientDescent
   , conjugateGradientAscent
   , stochasticGradientDescent
+
+  -- * Working with towers
+  , Jet(..)
+  , headJet
+  , tailJet
+  , jet
   ) where
 
 import Data.Functor.Compose
@@ -143,6 +149,7 @@
 import Numeric.AD.Internal.Sparse (Sparse, Grads, vgrads)
 
 import Numeric.AD.Internal.Type
+import Numeric.AD.Jet
 import Numeric.AD.Mode
 
 import qualified Numeric.AD.Rank1.Forward as Forward1
@@ -176,7 +183,7 @@
 --
 -- >>> import Numeric.AD.Internal.Doctest
 
--- | @'hessianProduct' f wv@ computes the product of the hessian @H@ of a non-scalar-to-scalar function @f@ at @w = 'fst' <$> wv@ with a vector @v = snd <$> wv@ using \"Pearlmutter\'s method\" from <http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.29.6143>, which states:
+-- | @'hessianProduct' f wv@ computes the product of the hessian @H@ of a non-scalar-to-scalar function @f@ at @w = 'fst' '<$>' wv@ with a vector @v = snd '<$>' wv@ using \"Pearlmutter\'s method\" from <http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.29.6143>, which states:
 --
 -- > H v = (d/dr) grad_w (w + r v) | r = 0
 --
@@ -185,7 +192,7 @@
 hessianProduct :: (Traversable f, Num a) => (forall s. Reifies s Tape => f (On (Reverse s (Forward a))) -> On (Reverse s (Forward a))) -> f (a, a) -> f a
 hessianProduct f = Forward1.duF (grad (off . f . fmap On))
 
--- | @'hessianProduct'' f wv@ computes both the gradient of a non-scalar-to-scalar @f@ at @w = 'fst' <$> wv@ and the product of the hessian @H@ at @w@ with a vector @v = snd <$> wv@ using \"Pearlmutter's method\". The outputs are returned wrapped in the same functor.
+-- | @'hessianProduct'' f wv@ computes both the gradient of a non-scalar-to-scalar @f@ at @w = 'fst' '<$>' wv@ and the product of the hessian @H@ at @w@ with a vector @v = snd '<$>' wv@ using \"Pearlmutter's method\". The outputs are returned wrapped in the same functor.
 --
 -- > H v = (d/dr) grad_w (w + r v) | r = 0
 --
diff --git a/src/Numeric/AD/Internal/Forward/Double.hs b/src/Numeric/AD/Internal/Forward/Double.hs
--- a/src/Numeric/AD/Internal/Forward/Double.hs
+++ b/src/Numeric/AD/Internal/Forward/Double.hs
@@ -194,12 +194,12 @@
 instance Erf ForwardDouble where
   erf = lift1 erf $ \x -> (2 / sqrt pi) * exp (negate x * x)
   erfc = lift1 erfc $ \x -> ((-2) / sqrt pi) * exp (negate x * x)
-  normcdf = lift1 normcdf $ \x -> ((-1) / sqrt pi) * exp (x * x * fromRational (- recip 2) / sqrt 2)
+  normcdf = lift1 normcdf $ \x -> (recip $ sqrt (2 * pi)) * exp (- x * x / 2)
 
 instance InvErf ForwardDouble where
-  inverf = lift1 inverfc $ \x -> recip $ (2 / sqrt pi) * exp (negate x * x)
-  inverfc = lift1 inverfc $ \x -> recip $ negate (2 / sqrt pi) * exp (negate x * x)
-  invnormcdf = lift1 invnormcdf $ \x -> recip $ ((-1) / sqrt pi) * exp (x * x * fromRational (- recip 2) / sqrt 2)
+  inverf = lift1_ inverf $ \x _ -> sqrt pi / 2 * exp (x * x)
+  inverfc = lift1_ inverfc $ \x _ -> negate (sqrt pi / 2) * exp (x * x)
+  invnormcdf = lift1_ invnormcdf $ \x _ -> sqrt (2 * pi) * exp (x * x / 2)
 
 bind :: Traversable f => (f ForwardDouble -> b) -> f Double -> f b
 bind f as = snd $ mapAccumL outer (0 :: Int) as where
diff --git a/src/Numeric/AD/Internal/Kahn.hs b/src/Numeric/AD/Internal/Kahn.hs
--- a/src/Numeric/AD/Internal/Kahn.hs
+++ b/src/Numeric/AD/Internal/Kahn.hs
@@ -231,9 +231,9 @@
   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 :: Tape a Int -> [Int]
   successors (Unary _ _ b) = [b]
-  successors (Binary _ _ _ b c) = [b,c]
+  successors (Binary _ _ _ b c) = if b == c then [b] else [b,c]
   successors _ = []
 
 -- | Return an 'Array' of 'partials' given bounds for the variable IDs.
diff --git a/src/Numeric/AD/Internal/Reverse/Double.hs b/src/Numeric/AD/Internal/Reverse/Double.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/AD/Internal/Reverse/Double.hs
@@ -0,0 +1,356 @@
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE DeriveDataTypeable #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE Rank2Types #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE UndecidableInstances #-}
+{-# OPTIONS_GHC -fno-full-laziness #-}
+{-# OPTIONS_HADDOCK not-home #-}
+
+module Numeric.AD.Internal.Reverse.Double
+  ( ReverseDouble(..)
+  , Tape(..)
+  , Head(..)
+  , Cells(..)
+  , reifyTape
+  , partials
+  , partialArrayOf
+  , partialMapOf
+  , derivativeOf
+  , derivativeOf'
+  , bind
+  , unbind
+  , unbindMap
+  , unbindWith
+  , unbindMapWithDefault
+  , var
+  , varId
+  , primal
+  ) where
+
+import Data.Functor
+import Control.Monad hiding (mapM)
+import Control.Monad.ST
+import Control.Monad.Trans.State
+import Data.Array.ST
+import Data.Array
+import Data.Array.Unsafe as Unsafe
+import Data.IORef
+import Data.IntMap (IntMap, fromDistinctAscList, findWithDefault)
+import Data.Number.Erf
+import Data.Proxy
+import Data.Reflection
+#if __GLASGOW_HASKELL__ < 710
+import Data.Traversable (Traversable, mapM)
+#else
+import Data.Traversable (mapM)
+#endif
+import Data.Typeable
+import Numeric.AD.Internal.Combinators
+import Numeric.AD.Internal.Identity
+import Numeric.AD.Jacobian
+import Numeric.AD.Mode
+import Prelude hiding (mapM)
+import System.IO.Unsafe (unsafePerformIO)
+import Unsafe.Coerce
+
+#ifdef HLINT
+{-# ANN module "HLint: ignore Reduce duplication" #-}
+#endif
+
+-- evil untyped tape
+#ifndef HLINT
+data Cells where
+  Nil    :: Cells
+  Unary  :: {-# UNPACK #-} !Int -> Double -> Cells -> Cells
+  Binary :: {-# UNPACK #-} !Int -> {-# UNPACK #-} !Int -> Double -> Double -> Cells -> Cells
+#endif
+
+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 -> Double -> 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 -> Double -> Double -> 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. Reifies s Tape => (Double -> Double) -> Double -> Int -> Double -> ReverseDouble s
+unarily f di i b = ReverseDouble (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. Reifies s Tape => (Double -> Double -> Double) -> Double -> Double -> Int -> Double -> Int -> Double -> ReverseDouble s
+binarily f di dj i b j c = ReverseDouble (unsafePerformIO (modifyTape (Proxy :: Proxy s) (bin i j di dj))) $! f b c
+{-# INLINE binarily #-}
+
+#ifndef HLINT
+data ReverseDouble s where
+  Zero :: ReverseDouble s
+  Lift :: Double -> ReverseDouble s
+  ReverseDouble :: {-# UNPACK #-} !Int -> Double -> ReverseDouble s
+  deriving (Show, Typeable)
+#endif
+
+instance (Reifies s Tape) => Mode (ReverseDouble s) where
+  type Scalar (ReverseDouble s) = Double
+
+  isKnownZero Zero = True
+  isKnownZero _    = False
+
+  isKnownConstant ReverseDouble{} = False
+  isKnownConstant _ = True
+
+  auto = Lift
+  zero = Zero
+  a *^ b = lift1 (a *) (\_ -> auto a) b
+  a ^* b = lift1 (* b) (\_ -> auto b) a
+  a ^/ b = lift1 (/ b) (\_ -> auto (recip b)) a
+
+(<+>) :: (Reifies s Tape) => ReverseDouble s -> ReverseDouble s -> ReverseDouble s
+(<+>)  = binary (+) 1 1
+
+(<**>) :: (Reifies s Tape) => ReverseDouble s -> ReverseDouble s -> ReverseDouble s
+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 * xi ** (yi - 1), z * log xi)) x y
+
+primal :: ReverseDouble s -> Double
+primal Zero = 0
+primal (Lift a) = a
+primal (ReverseDouble _ a) = a
+
+instance (Reifies s Tape) => Jacobian (ReverseDouble s) where
+  type D (ReverseDouble s) = Id Double
+
+  unary f _         (Zero)   = Lift (f 0)
+  unary f _         (Lift a) = Lift (f a)
+  unary f (Id dadi) (ReverseDouble 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        (ReverseDouble i c) = unarily (f 0) dadc i c
+  binary f _         (Id dadc) (Lift b)    (ReverseDouble i c) = unarily (f b) dadc i c
+  binary f (Id dadb) _         (ReverseDouble i b) Zero        = unarily (`f` 0) dadb i b
+  binary f (Id dadb) _         (ReverseDouble i b) (Lift c)    = unarily (`f` c) dadb i b
+  binary f (Id dadb) (Id dadc) (ReverseDouble i b) (ReverseDouble 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)
+
+
+instance (Reifies s Tape) => Eq (ReverseDouble s) where
+  a == b = primal a == primal b
+
+instance (Reifies s Tape) => Ord (ReverseDouble s) where
+  compare a b = compare (primal a) (primal b)
+
+instance (Reifies s Tape) => Num (ReverseDouble s) where
+  fromInteger 0  = zero
+  fromInteger n = auto (fromInteger n)
+  (+)          = (<+>) -- binary (+) 1 1
+  (-)          = binary (-) (auto 1) (auto (-1)) -- TODO: <-> ? as it is, this might be pretty bad for Tower
+  (*)          = lift2 (*) (\x y -> (y, x))
+  negate       = lift1 negate (const (auto (-1)))
+  abs          = lift1 abs signum
+  signum a     = lift1 signum (const zero) a
+
+instance (Reifies s Tape) => Fractional (ReverseDouble s) where
+  fromRational 0 = zero
+  fromRational r = auto (fromRational r)
+  x / y        = x * recip y
+  recip        = lift1_ recip (const . negate . join (*))
+
+instance (Reifies s Tape) => Floating (ReverseDouble s) where
+  pi       = auto pi
+  exp      = lift1_ exp const
+  log      = lift1 log recip
+  logBase x y = log y / log x
+  sqrt     = lift1_ sqrt (\z _ -> recip (auto 2 * z))
+  (**)     = (<**>)
+  --x ** y
+  --   | isKnownZero y     = 1
+  --   | isKnownConstant y, y' <- primal y = lift1 (** y') ((y'*) . (**(y'-1))) x
+  --   | otherwise         = lift2_ (**) (\z xi yi -> (yi * z / xi, z * log1 xi)) x y
+  sin      = lift1 sin cos
+  cos      = lift1 cos $ negate . sin
+  tan      = lift1 tan $ recip . join (*) . cos
+  asin     = lift1 asin $ \x -> recip (sqrt (auto 1 - join (*) x))
+  acos     = lift1 acos $ \x -> negate (recip (sqrt (1 - join (*) x)))
+  atan     = lift1 atan $ \x -> recip (1 + join (*) x)
+  sinh     = lift1 sinh cosh
+  cosh     = lift1 cosh sinh
+  tanh     = lift1 tanh $ recip . join (*) . cosh
+  asinh    = lift1 asinh $ \x -> recip (sqrt (1 + join (*) x))
+  acosh    = lift1 acosh $ \x -> recip (sqrt (join (*) x - 1))
+  atanh    = lift1 atanh $ \x -> recip (1 - join (*) x)
+
+instance (Reifies s Tape) => Enum (ReverseDouble s) where
+  succ             = lift1 succ (const 1)
+  pred             = lift1 pred (const 1)
+  toEnum           = auto . toEnum
+  fromEnum a       = fromEnum (primal a)
+  enumFrom a       = withPrimal a <$> enumFrom (primal a)
+  enumFromTo a b   = withPrimal a <$> enumFromTo (primal a) (primal b)
+  enumFromThen a b = zipWith (fromBy a delta) [0..] $ enumFromThen (primal a) (primal b) where delta = b - a
+  enumFromThenTo a b c = zipWith (fromBy a delta) [0..] $ enumFromThenTo (primal a) (primal b) (primal c) where delta = b - a
+
+instance (Reifies s Tape) => Real (ReverseDouble s) where
+  toRational      = toRational . primal
+
+instance (Reifies s Tape) => RealFloat (ReverseDouble s) where
+  floatRadix      = floatRadix . primal
+  floatDigits     = floatDigits . primal
+  floatRange      = floatRange . primal
+  decodeFloat     = decodeFloat . primal
+  encodeFloat m e = auto (encodeFloat m e)
+  isNaN           = isNaN . primal
+  isInfinite      = isInfinite . primal
+  isDenormalized  = isDenormalized . primal
+  isNegativeZero  = isNegativeZero . primal
+  isIEEE          = isIEEE . primal
+  exponent = exponent
+  scaleFloat n = unary (scaleFloat n) (scaleFloat n 1)
+  significand x =  unary significand (scaleFloat (- floatDigits x) 1) x
+  atan2 = lift2 atan2 $ \vx vy -> let r = recip (join (*) vx + join (*) vy) in (vy * r, negate vx * r)
+
+instance (Reifies s Tape) => RealFrac (ReverseDouble s) where
+  properFraction a = (w, a `withPrimal` pb) where
+    pa = primal a
+    (w, pb) = properFraction pa
+  truncate = truncate . primal
+  round    = round . primal
+  ceiling  = ceiling . primal
+  floor    = floor . primal
+
+instance (Reifies s Tape) => Erf (ReverseDouble s) where
+  erf = lift1 erf $ \x -> (2 / sqrt pi) * exp (negate x * x)
+  erfc = lift1 erfc $ \x -> ((-2) / sqrt pi) * exp (negate x * x)
+  normcdf = lift1 normcdf $ \x -> (recip $ sqrt (2 * pi)) * exp (- x * x / 2)
+
+instance (Reifies s Tape) => InvErf (ReverseDouble s) where
+  inverf = lift1_ inverf $ \x _ -> sqrt pi / 2 * exp (x * x)
+  inverfc = lift1_ inverfc $ \x _ -> negate (sqrt pi / 2) * exp (x * x)
+  invnormcdf = lift1_ invnormcdf $ \x _ -> sqrt (2 * pi) * exp (x * x / 2)
+
+-- | Helper that extracts the derivative of a chain when the chain was constructed with 1 variable.
+derivativeOf :: (Reifies s Tape) => Proxy s -> ReverseDouble s -> Double
+derivativeOf _ = sum . partials
+{-# INLINE derivativeOf #-}
+
+-- | Helper that extracts both the primal and derivative of a chain when the chain was constructed with 1 variable.
+derivativeOf' :: (Reifies s Tape) => Proxy s -> ReverseDouble s -> (Double, Double)
+derivativeOf' p r = (primal r, derivativeOf p r)
+{-# INLINE derivativeOf' #-}
+
+-- | Used internally to push sensitivities down the chain.
+backPropagate :: Int -> Cells -> STArray s Int Double -> 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.
+partials :: forall s. (Reifies s Tape) => ReverseDouble s -> [Double]
+partials Zero        = []
+partials (Lift _)    = []
+partials (ReverseDouble 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) => Proxy s -> (Int, Int) -> ReverseDouble s -> Array Int Double
+partialArrayOf _ vbounds = accumArray (+) 0 vbounds . zip [0..] . partials
+{-# INLINE partialArrayOf #-}
+
+-- | Return an 'IntMap' of sparse partials
+partialMapOf :: (Reifies s Tape) => Proxy s -> ReverseDouble s-> IntMap Double
+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 #-}
+
+var :: Double -> Int -> ReverseDouble s
+var a v = ReverseDouble v a
+
+varId :: ReverseDouble s -> Int
+varId (ReverseDouble v _) = v
+varId _ = error "varId: not a Var"
+
+bind :: Traversable f => f Double -> (f (ReverseDouble s), (Int,Int))
+bind xs = (r,(0,hi)) where
+  (r,hi) = runState (mapM freshVar xs) 0
+  freshVar a = state $ \s -> let s' = s + 1 in s' `seq` (var a s, s')
+
+unbind :: Functor f => f (ReverseDouble s) -> Array Int Double -> f Double
+unbind xs ys = fmap (\v -> ys ! varId v) xs
+
+unbindWith :: (Functor f) => (Double -> b -> c) -> f (ReverseDouble s) -> Array Int b -> f c
+unbindWith f xs ys = fmap (\v -> f (primal v) (ys ! varId v)) xs
+
+unbindMap :: (Functor f) => f (ReverseDouble s) -> IntMap Double -> f Double
+unbindMap xs ys = fmap (\v -> findWithDefault 0 (varId v) ys) xs
+
+unbindMapWithDefault :: (Functor f) => b -> (Double -> b -> c) -> f (ReverseDouble s) -> IntMap b -> f c
+unbindMapWithDefault z f xs ys = fmap (\v -> f (primal v) $ findWithDefault z (varId v) ys) xs
diff --git a/src/Numeric/AD/Jet.hs b/src/Numeric/AD/Jet.hs
--- a/src/Numeric/AD/Jet.hs
+++ b/src/Numeric/AD/Jet.hs
@@ -36,7 +36,7 @@
 import Data.Typeable
 import Control.Comonad.Cofree
 
-infixl 3 :-
+infixr 3 :-
 
 -- | A 'Jet' is a tower of all (higher order) partial derivatives of a function
 --
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
@@ -150,21 +150,21 @@
 
 -- | Compute the gradient of a function using forward mode AD.
 --
--- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.grad' for @n@ inputs, in exchange for better space utilization.
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.grad' for @n@ inputs, in exchange for better space utilization.
 grad :: (Traversable f, Num a) => (forall s. f (AD s (Forward a)) -> AD s (Forward a)) -> f a -> f a
 grad f = Rank1.grad (runAD.f.fmap AD)
 {-# INLINE grad #-}
 
 -- | Compute the gradient and answer to a function using forward mode AD.
 --
--- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.grad'' for @n@ inputs, in exchange for better space utilization.
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.grad'' for @n@ inputs, in exchange for better space utilization.
 grad' :: (Traversable f, Num a) => (forall s. f (AD s (Forward a)) -> AD s (Forward a)) -> f a -> (a, f a)
 grad' f = Rank1.grad' (runAD.f.fmap AD)
 {-# INLINE grad' #-}
 
 -- | Compute the gradient of a function using forward mode AD and combine the result with the input using a user-specified function.
 --
--- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.gradWith' for @n@ inputs, in exchange for better space utilization.
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.gradWith' for @n@ inputs, in exchange for better space utilization.
 gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. f (AD s (Forward a)) -> AD s (Forward a)) -> f a -> f b
 gradWith g f = Rank1.gradWith g (runAD.f.fmap AD)
 {-# INLINE gradWith #-}
@@ -172,7 +172,7 @@
 -- | 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.Wengert.gradWith'' for @n@ inputs, in exchange for better space utilization.
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.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)])
diff --git a/src/Numeric/AD/Mode/Forward/Double.hs b/src/Numeric/AD/Mode/Forward/Double.hs
--- a/src/Numeric/AD/Mode/Forward/Double.hs
+++ b/src/Numeric/AD/Mode/Forward/Double.hs
@@ -144,21 +144,21 @@
 
 -- | Compute the gradient of a function using forward mode AD.
 --
--- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.grad' for @n@ inputs, in exchange for better space utilization.
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.grad' for @n@ inputs, in exchange for better space utilization.
 grad :: Traversable f => (forall s. f (AD s ForwardDouble) -> AD s ForwardDouble) -> f Double -> f Double
 grad f = Rank1.grad (runAD.f.fmap AD)
 {-# INLINE grad #-}
 
 -- | Compute the gradient and answer to a function using forward mode AD.
 --
--- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.grad'' for @n@ inputs, in exchange for better space utilization.
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.grad'' for @n@ inputs, in exchange for better space utilization.
 grad' :: Traversable f => (forall s. f (AD s ForwardDouble) -> AD s ForwardDouble) -> f Double -> (Double, f Double)
 grad' f = Rank1.grad' (runAD.f.fmap AD)
 {-# INLINE grad' #-}
 
 -- | Compute the gradient of a function using forward mode AD and combine the result with the input using a user-specified function.
 --
--- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.gradWith' for @n@ inputs, in exchange for better space utilization.
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.gradWith' for @n@ inputs, in exchange for better space utilization.
 gradWith :: Traversable f => (Double -> Double -> b) -> (forall s. f (AD s ForwardDouble) -> AD s ForwardDouble) -> f Double -> f b
 gradWith g f = Rank1.gradWith g (runAD.f.fmap AD)
 {-# INLINE gradWith #-}
@@ -166,7 +166,7 @@
 -- | 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.Wengert.gradWith'' for @n@ inputs, in exchange for better space utilization.
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.gradWith'' for @n@ inputs, in exchange for better space utilization.
 --
 -- >>> gradWith' (,) sum [0..4]
 -- (10.0,[(0.0,1.0),(1.0,1.0),(2.0,1.0),(3.0,1.0),(4.0,1.0)])
diff --git a/src/Numeric/AD/Mode/Reverse/Double.hs b/src/Numeric/AD/Mode/Reverse/Double.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/AD/Mode/Reverse/Double.hs
@@ -0,0 +1,216 @@
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE Rank2Types #-}
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+-----------------------------------------------------------------------------
+-- |
+-- Copyright   :  (c) Edward Kmett 2010-2015
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only
+--
+-- Reverse-mode automatic differentiation using Wengert lists and
+-- Data.Reflection
+--
+-- This version is specialized to `Double` enabling the entire
+-- structure
+--
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Mode.Reverse.Double
+  ( ReverseDouble, auto
+  -- * Gradient
+  , grad
+  , grad'
+  , gradWith
+  , gradWith'
+
+  -- * Jacobian
+  , jacobian
+  , jacobian'
+  , jacobianWith
+  , jacobianWith'
+
+  -- * Hessian
+  , hessian
+  , hessianF
+
+  -- * Derivatives
+  , diff
+  , diff'
+  , diffF
+  , diffF'
+  ) where
+
+#if __GLASGOW_HASKELL__ < 710
+import Data.Functor ((<$>))
+import Data.Traversable (Traversable)
+#endif
+import Data.Functor.Compose
+import Data.Reflection (Reifies)
+import Numeric.AD.Internal.On
+import qualified Numeric.AD.Internal.Reverse as R
+import qualified Numeric.AD.Mode.Reverse as M
+import Numeric.AD.Internal.Reverse.Double
+import Numeric.AD.Mode
+
+-- $setup
+--
+-- >>> import Numeric.AD.Internal.Doctest
+
+-- | 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.0,1.0,1.0]
+--
+-- >>> grad (\[x,y] -> x**y) [0,2]
+-- [0.0,NaN]
+grad :: (Traversable f) => (forall s. Reifies s Tape => f (ReverseDouble s) -> ReverseDouble s) -> f Double -> f Double
+grad f as = reifyTape (snd bds) $ \p -> unbind vs $! partialArrayOf p 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] -> x*y+z) [1,2,3]
+-- (5.0,[2.0,1.0,1.0])
+grad' :: (Traversable f) => (forall s. Reifies s Tape => f (ReverseDouble s) -> ReverseDouble s) -> f Double -> (Double, f Double)
+grad' f as = reifyTape (snd bds) $ \p -> case f vs of
+   r -> (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) => (Double -> Double -> b) -> (forall s. Reifies s Tape => f (ReverseDouble s) -> ReverseDouble s) -> f Double -> f b
+gradWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs $! partialArrayOf p 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) => (Double -> Double -> b) -> (forall s. Reifies s Tape => f (ReverseDouble s) -> ReverseDouble s) -> f Double -> (Double, f b)
+gradWith' g f as = reifyTape (snd bds) $ \p -> case f vs of
+   r -> (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.0,1.0],[1.0,0.0],[1.0,2.0]]
+jacobian :: (Traversable f, Functor g) => (forall s. Reifies s Tape => f (ReverseDouble s) -> g (ReverseDouble s)) -> f Double -> g (f Double)
+jacobian f as = reifyTape (snd bds) $ \p -> unbind vs . partialArrayOf p 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''
+--
+-- >>> jacobian' (\[x,y] -> [y,x,x*y]) [2,1]
+-- [(1.0,[0.0,1.0]),(2.0,[1.0,0.0]),(2.0,[1.0,2.0])]
+jacobian' :: (Traversable f, Functor g) => (forall s. Reifies s Tape => f (ReverseDouble s) -> g (ReverseDouble s)) -> f Double -> g (Double, f Double)
+jacobian' f as = reifyTape (snd bds) $ \p ->
+  let row a = (primal a, unbind vs $! partialArrayOf p bds $! a)
+  in row <$> f 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) => (Double -> Double -> b) -> (forall s. Reifies s Tape => f (ReverseDouble s) -> g (ReverseDouble s)) -> f Double -> g (f b)
+jacobianWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs . partialArrayOf p 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) => (Double -> Double -> b) -> (forall s. Reifies s Tape => f (ReverseDouble s) -> g (ReverseDouble s)) -> f Double -> g (Double, 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 vs
+  where (vs, bds) = bind as
+{-# INLINE jacobianWith' #-}
+
+-- | Compute the derivative of a function.
+--
+-- >>> diff sin 0
+-- 1.0
+diff :: (forall s. Reifies s Tape => ReverseDouble s -> ReverseDouble s) -> Double -> Double
+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' :: (forall s. Reifies s Tape => ReverseDouble s -> ReverseDouble s) -> Double -> (Double, Double)
+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) => (forall s. Reifies s Tape => ReverseDouble s -> f (ReverseDouble s)) -> Double -> f Double
+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) => (forall s. Reifies s Tape => ReverseDouble s -> f (ReverseDouble s)) -> Double -> f (Double, Double)
+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.0,1.0],[1.0,0.0]]
+hessian :: (Traversable f) => (forall s s'. (Reifies s R.Tape, Reifies s' Tape) => f (On (R.Reverse s (ReverseDouble s'))) -> (On (R.Reverse s (ReverseDouble s')))) -> f Double -> f (f Double)
+hessian f = jacobian (M.grad (off . f . fmap On))
+{-# INLINE hessian #-}
+
+-- | 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 :: Double]
+-- [[[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) => (forall s s'. (Reifies s R.Tape, Reifies s' Tape) => f (On (R.Reverse s (ReverseDouble s'))) -> g (On (R.Reverse s (ReverseDouble s')))) -> f Double -> g (f (f Double))
+hessianF f = getCompose . jacobian (Compose . M.jacobian (fmap off . f . fmap On))
+{-# INLINE hessianF #-}
diff --git a/src/Numeric/AD/Rank1/Forward.hs b/src/Numeric/AD/Rank1/Forward.hs
--- a/src/Numeric/AD/Rank1/Forward.hs
+++ b/src/Numeric/AD/Rank1/Forward.hs
@@ -160,14 +160,14 @@
 
 -- | Compute the gradient of a function using forward mode AD.
 --
--- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.grad' for @n@ inputs, in exchange for better space utilization.
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.grad' for @n@ inputs, in exchange for better space utilization.
 grad :: (Traversable f, Num a) => (f (Forward a) -> Forward 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.Wengert.grad'' for @n@ inputs, in exchange for better space utilization.
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.grad'' for @n@ inputs, in exchange for better space utilization.
 grad' :: (Traversable f, Num a) => (f (Forward a) -> Forward a) -> f a -> (a, f a)
 grad' f as = (primal b, tangent <$> bs) where
   (b, bs) = bind' f as
@@ -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.Wengert.gradWith' for @n@ inputs, in exchange for better space utilization.
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.gradWith' for @n@ inputs, in exchange for better space utilization.
 gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (f (Forward a) -> Forward a) -> f a -> f b
 gradWith g f = bindWith g (tangent . f)
 {-# INLINE gradWith #-}
@@ -183,7 +183,7 @@
 -- | 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.Wengert.gradWith'' for @n@ inputs, in exchange for better space utilization.
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.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)])
diff --git a/src/Numeric/AD/Rank1/Forward/Double.hs b/src/Numeric/AD/Rank1/Forward/Double.hs
--- a/src/Numeric/AD/Rank1/Forward/Double.hs
+++ b/src/Numeric/AD/Rank1/Forward/Double.hs
@@ -139,14 +139,14 @@
 
 -- | Compute the gradient of a function using forward mode AD.
 --
--- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.grad' for @n@ inputs, in exchange for better space utilization.
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.grad' for @n@ inputs, in exchange for better space utilization.
 grad :: Traversable f => (f ForwardDouble -> ForwardDouble) -> f Double -> f Double
 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.Wengert.grad'' for @n@ inputs, in exchange for better space utilization.
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.grad'' for @n@ inputs, in exchange for better space utilization.
 grad' :: Traversable f => (f ForwardDouble -> ForwardDouble) -> f Double -> (Double, f Double)
 grad' f as = (primal b, tangent <$> bs)
     where
@@ -155,7 +155,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.Wengert.gradWith' for @n@ inputs, in exchange for better space utilization.
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.gradWith' for @n@ inputs, in exchange for better space utilization.
 gradWith :: Traversable f => (Double -> Double -> b) -> (f ForwardDouble -> ForwardDouble) -> f Double -> f b
 gradWith g f = bindWith g (tangent . f)
 {-# INLINE gradWith #-}
@@ -163,7 +163,7 @@
 -- | 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.Wengert.gradWith'' for @n@ inputs, in exchange for better space utilization.
+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.gradWith'' for @n@ inputs, in exchange for better space utilization.
 --
 -- >>> gradWith' (,) sum [0..4]
 -- (10.0,[(0.0,1.0),(1.0,1.0),(2.0,1.0),(3.0,1.0),(4.0,1.0)])