diff --git a/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -1,4 +1,6 @@
 Copyright (c) 2010, Edward Kmett
+          (c) 2008-2009 Barak A. Pearlmutter and Jeffrey Mark Siskind
+
 All rights reserved.
 
 Redistribution and use in source and binary forms, with or without
diff --git a/Numeric/AD.hs b/Numeric/AD.hs
--- a/Numeric/AD.hs
+++ b/Numeric/AD.hs
@@ -17,9 +17,11 @@
     (
     -- * Gradients
       grad, grad2
+    , gradWith, gradWith2
 
     -- * Jacobians
     , jacobian, jacobian2
+    , jacobianWith, jacobianWith2
 
     -- * Synonyms
     , diff
@@ -61,7 +63,7 @@
 import Numeric.AD.Internal (AD(..), probed, unprobe)
 import Numeric.AD.Forward  (diff, diffUU, diff2, diff2UU, diffUF, diff2UF)
 import Numeric.AD.Tower    (diffsUU, diffs0UU , diffsUF, diffs0UF , diffs, diffs0, taylor, taylor0)
-import Numeric.AD.Reverse  (diffFU, diff2FU, grad, grad2)
+import Numeric.AD.Reverse  (diffFU, diff2FU, grad, grad2, gradWith, gradWith2)
 
 import qualified Numeric.AD.Forward as Forward
 import qualified Numeric.AD.Reverse as Reverse
@@ -83,3 +85,26 @@
         size :: Foldable f => f a -> Int
         size = foldr' (\_ b -> 1 + b) 0
 {-# INLINE jacobian2 #-}
+
+-- | @'jacobianWith' g f@ calculates the Jacobian of a non-scalar-to-non-scalar function, automatically choosing between forward and reverse mode AD based on the number of inputs and outputs.
+-- 
+-- The resulting Jacobian matrix is then recombined element-wise with the input using @g@. 
+jacobianWith :: (Traversable f, Traversable 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 bs = snd <$> jacobianWith2 g f bs
+{-# INLINE jacobianWith #-}
+
+-- | @'jacobianWith2' g f@ calculates the answer and Jacobian of a non-scalar-to-non-scalar function, automatically choosing between forward and reverse mode AD based on the number of inputs and outputs.
+--
+-- The resulting Jacobian matrix is then recombined element-wise with the input using @g@. 
+jacobianWith2 :: (Traversable f, Traversable g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f b)
+jacobianWith2 g f bs 
+    | n == 0    = fmap (\x -> (unprobe x, undefined <$> bs)) as
+    | n > m     = Reverse.jacobianWith2 g f bs
+    | otherwise = Forward.jacobianWith2 g f bs
+    where
+        as = f (probed bs)
+        n = size bs
+        m = size as
+        size :: Foldable f => f a -> Int
+        size = foldr' (\_ b -> 1 + b) 0
+{-# INLINE jacobianWith2 #-}
diff --git a/Numeric/AD/Forward.hs b/Numeric/AD/Forward.hs
--- a/Numeric/AD/Forward.hs
+++ b/Numeric/AD/Forward.hs
@@ -17,15 +17,23 @@
     -- * Gradient
       grad
     , grad2
+    , gradWith
+    , gradWith2
     -- * Jacobian
     , jacobian
     , jacobian2
     , jacobianT
+    , jacobianWith
+    , jacobianWith2
+    , jacobianWithT
     -- * Derivatives
     , diffUU
     , diff2UU
     , diffUF
     , diff2UF
+    -- * Directional Derivatives
+    , diffMU 
+    , diff2MU
     -- * Synonyms
     , diff
     , diff2
@@ -40,6 +48,14 @@
 import Numeric.AD.Internal
 import Numeric.AD.Internal.Forward
 
+diffMU :: (Functor f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> a
+diffMU f = tangent . f . fmap (uncurry bundle)
+{-# INLINE diffMU #-}
+
+diff2MU :: (Functor f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> (a, a)
+diff2MU f = unbundle . f . fmap (uncurry bundle)
+{-# INLINE diff2MU #-}
+
 -- | The 'diff2' function calculates the first derivative of scalar-to-scalar function by 'Forward' 'AD'
 diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a
 diff = diffUU
@@ -70,18 +86,31 @@
 diff2UF f a = unbundle <$> apply f a
 {-# INLINE diff2UF #-}
 
--- A fast, simple transposed forward jacobian
+-- | A fast, simple transposed Jacobian computed with forward-mode AD.
 jacobianT :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> f (g a)
 jacobianT f = bind (fmap tangent . f)
 -- jacobianT f as = fmap tangent <$> bind f as
 {-# INLINE jacobianT #-}
 
+-- | A fast, simple transposed Jacobian computed with forward-mode AD.
+jacobianWithT :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> f (g b)
+jacobianWithT g f = bindWith g' f 
+    where g' a ga = g a . tangent <$> ga 
+{-# INLINE jacobianWithT #-}
+
 jacobian :: (Traversable f, Traversable g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)
 jacobian f as = transposeWith (const id) t p
     where
         (p, t) = bind2 (fmap tangent . f) as
 {-# INLINE jacobian #-}
 
+jacobianWith :: (Traversable f, Traversable 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 = transposeWith (const id) t p
+    where
+        (p, t) = bindWith2 g' f as
+        g' a ga = g a . tangent <$> ga 
+{-# INLINE jacobianWith #-}
+
 jacobian2 :: (Traversable f, Traversable g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)
 jacobian2 f as = transposeWith row t p
     where
@@ -89,12 +118,29 @@
         row x as' = (primal x, tangent <$> as')
 {-# INLINE jacobian2 #-}
 
+jacobianWith2 :: (Traversable f, Traversable g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f b)
+jacobianWith2 g f as = transposeWith row t p
+    where
+        (p, t) = bindWith2 g' f as
+        row x as' = (primal x, as')
+        g' a ga = g a . tangent <$> ga 
+{-# INLINE jacobianWith2 #-}
+
 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 #-}
 
+
 grad2 :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)
 grad2 f as = (primal b, tangent <$> bs)
     where
         (b, bs) = bind2 f as
 {-# INLINE grad2 #-}
+
+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 #-}
+
+gradWith2 :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f b)
+gradWith2 g f = bindWith2 g (tangent . f)
+{-# INLINE gradWith2 #-}
diff --git a/Numeric/AD/Internal.hs b/Numeric/AD/Internal.hs
--- a/Numeric/AD/Internal.hs
+++ b/Numeric/AD/Internal.hs
@@ -11,6 +11,7 @@
 -----------------------------------------------------------------------------
 module Numeric.AD.Internal
     ( zipWithT
+    , zipWithDefaultT
     , AD(..)
     , Id(..)
     , probe
@@ -28,6 +29,9 @@
 
 zipWithT :: (Foldable f, Traversable g) => (a -> b -> c) -> f a -> g b -> g c
 zipWithT f as = snd . mapAccumL (\(a:as') b -> (as', f a b)) (toList as)
+
+zipWithDefaultT :: (Foldable f, Traversable g) => a -> (a -> b -> c) -> f a -> g b -> g c
+zipWithDefaultT z f as = zipWithT f (toList as ++ repeat z)
 
 class Iso a b where
     iso :: f a -> f b
diff --git a/Numeric/AD/Internal/Forward.hs b/Numeric/AD/Internal/Forward.hs
--- a/Numeric/AD/Internal/Forward.hs
+++ b/Numeric/AD/Internal/Forward.hs
@@ -21,6 +21,8 @@
     , apply
     , bind
     , bind2
+    , bindWith
+    , bindWith2
     , transposeWith
     ) where
 
@@ -98,6 +100,20 @@
 bind2 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, bundle a $ if i == j then 1 else 0)
+        b0 = f (lift <$> 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, bundle a $ if i == j then 1 else 0)
+
+bindWith2 :: (Traversable f, Num a) => (a -> b -> c) -> (f (AD Forward a) -> b) -> f a -> (b, f c)
+bindWith2 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, bundle a $ if i == j then 1 else 0)
         b0 = f (lift <$> as)
         dropIx ((_,b),bs) = (b,bs)
diff --git a/Numeric/AD/Internal/Reverse.hs b/Numeric/AD/Internal/Reverse.hs
--- a/Numeric/AD/Internal/Reverse.hs
+++ b/Numeric/AD/Internal/Reverse.hs
@@ -26,6 +26,11 @@
     , derivative
     , derivative2
     , Var(..)
+    , bind
+    , unbind
+    , unbindMap
+    , unbindWith
+    , unbindMapWithDefault
     ) where
 
 import Prelude hiding (mapM)
@@ -176,27 +181,33 @@
     S g >>= f = S (\s -> let (a,s') = g s in runS (f a) s')
 
 -- | Used to mark variables for inspection during the reverse pass
-class Var t a | t -> a where
-    var     :: a -> Int -> t
-    fromVar :: t -> Int
+class Primal v => Var v where
+    var   :: a -> Int -> v a
+    varId :: v a -> Int
 
-    bind :: Traversable f => f a -> (f t, (Int,Int))
-    unbind :: Functor f => f t -> Array Int a -> f a
-    unbindMap :: (Functor f, Num a) => f t -> IntMap a -> f a
+instance Var Reverse where
+    var a v = Reverse (Var a v)
+    varId (Reverse (Var _ v)) = v
+    varId _ = error "varId: not a Var"
 
-    -- TODO: tweak bounds
-    bind xs = (r,(0,hi))
-        where
+instance Var (AD Reverse) where
+    var a v = AD (var a v)
+    varId (AD v) = varId v
+
+bind :: (Traversable f, Var v) => f a -> (f (v a), (Int,Int))
+bind xs = (r,(0,hi))
+    where
         (r,hi) = runS (mapM freshVar xs) 0
         freshVar a = S (\s -> let s' = s + 1 in s' `seq` (var a s, s'))
-    unbind xs ys = fmap (\v -> ys ! fromVar v) xs
-    unbindMap xs ys = fmap (\v -> findWithDefault 0 (fromVar v) ys) xs
 
-instance Var (Reverse a) a where
-    var a v = Reverse (Var a v)
-    fromVar (Reverse (Var _ v)) = v
-    fromVar _ = error "fromVar: not a Var"
+unbind :: (Functor f, Var v)  => f (v a) -> Array Int a -> f a
+unbind xs ys = fmap (\v -> ys ! varId v) xs
 
-instance Var (AD Reverse a) a where
-    var a v = AD (var a v)
-    fromVar (AD v) = fromVar v
+unbindWith :: (Functor f, Var v, Num a) => (a -> b -> c) -> f (v a) -> Array Int b -> f c
+unbindWith f xs ys = fmap (\v -> f (primal v) (ys ! varId v)) xs 
+
+unbindMap :: (Functor f, Var v, Num a) => f (v a) -> IntMap a -> f a
+unbindMap xs ys = fmap (\v -> findWithDefault 0 (varId v) ys) xs
+
+unbindMapWithDefault :: (Functor f, Var v, Num a) => b -> (a -> b -> c) -> f (v a) -> IntMap b -> f c
+unbindMapWithDefault z f xs ys = fmap (\v -> f (primal v) $ findWithDefault z (varId v) ys) xs 
diff --git a/Numeric/AD/Newton.hs b/Numeric/AD/Newton.hs
--- a/Numeric/AD/Newton.hs
+++ b/Numeric/AD/Newton.hs
@@ -30,7 +30,7 @@
 import Data.Foldable (all)
 import Data.Traversable (Traversable)
 import Numeric.AD.Forward (diff, diff2)
-import Numeric.AD.Reverse (grad2)
+import Numeric.AD.Reverse (gradWith2)
 
 -- | The 'findZero' function finds a zero of a scalar function using
 -- Newton's method; its output is a stream of increasingly accurate
@@ -81,18 +81,19 @@
 --
 -- It uses reverse mode automatic differentiation to compute the gradient.
 gradientDescent :: (Traversable f, Fractional a, Ord a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> [f a]
-gradientDescent f x0 = go x0 fx0 gx0 0.1 (0 :: Int)
+gradientDescent f x0 = go x0 fx0 xgx0 0.1 (0 :: Int)
     where
-        (fx0, gx0) = grad2 f x0
-        go x fx gx !eta !i
+        (fx0, xgx0) = gradWith2 (,) f x0
+        go x fx xgx !eta !i
             | eta == 0     = [] -- step size is 0
-            | fx1 > fx     = go x fx gx (eta/2) 0 -- we stepped too far
-            | all (==0) gx = [] -- gradient is 0
+            | fx1 > fx     = go x fx xgx (eta/2) 0 -- we stepped too far
+            | zeroGrad xgx = [] -- gradient is 0
             | otherwise    = x1 : if i == 10
-                                  then go x1 fx1 gx1 (eta*2) 0
-                                  else go x1 fx1 gx1 eta (i+1)
+                                  then go x1 fx1 xgx1 (eta*2) 0
+                                  else go x1 fx1 xgx1 eta (i+1)
             where
-                -- should check gx = 0 here
-                x1 = zipWithT (\xi gxi -> xi - eta * gxi) x gx
-                (fx1, gx1) = grad2 f x1
+                zeroGrad = all (\(_,g) -> g == 0)
+                x1 = fmap (\(xi,gxi) -> xi - eta * gxi) xgx
+                (fx1, xgx1) = gradWith2 (,) f x1
+                
 {-# INLINE gradientDescent #-}
diff --git a/Numeric/AD/Reverse.hs b/Numeric/AD/Reverse.hs
--- a/Numeric/AD/Reverse.hs
+++ b/Numeric/AD/Reverse.hs
@@ -22,9 +22,13 @@
     -- * Gradient
       grad
     , grad2
+    , gradWith
+    , gradWith2
     -- * Jacobian
     , jacobian
     , jacobian2
+    , jacobianWith
+    , jacobianWith2
     -- * Derivatives
     , diffUU
     , diff2UU
@@ -60,6 +64,26 @@
           r = f vs
 {-# INLINE grad2 #-}
 
+-- | @'grad' g f@ function calculates the gradient of a non-scalar-to-scalar function @f@ with reverse-mode AD in a single pass.
+-- The gradient is combined element-wise with the argument using the function @g@.
+--
+-- > grad == gradWith (\_ dx -> dx) 
+-- > id == gradWith const
+gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f b
+gradWith g f as = unbindWith g vs (partialArray bds $ f vs)
+    where (vs,bds) = bind as
+{-# INLINE gradWith #-}
+
+-- | @'grad2' g f@ calculates the result and gradient of a non-scalar-to-scalar function @f@ with 'Reverse' AD in a single pass
+-- the gradient is combined element-wise with the argument using the function @g@.
+-- 
+-- > grad2 == gradWith2 (\_ dx -> dx)
+gradWith2 :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f b)
+gradWith2 g f as = (primal r, unbindWith g vs $ partialArray bds r)
+    where (vs, bds) = bind as
+          r = f vs
+{-# INLINE gradWith2 #-}
+
 -- | The 'jacobian' function calculates the jacobian of a non-scalar-to-non-scalar function with reverse AD lazily in @m@ passes for @m@ outputs.
 jacobian :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)
 jacobian f as = unbind vs . partialArray bds <$> f vs where
@@ -73,6 +97,31 @@
     (vs, bds) = bind as
     row a = (primal a, unbind vs (partialArray bds a))
 {-# INLINE jacobian2 #-}
+
+-- | 'jacobianWith g f' calculates the jacobian of a non-scalar-to-non-scalar function @f@ with reverse AD lazily in @m@ passes for @m@ outputs.
+--
+-- Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@.
+-- 
+-- > jacobian == jacobianWith (\_ dx -> dx)
+-- > jacobianWith const == (\f x -> const x <$> f x)
+--
+jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f b)
+jacobianWith g f as = unbindWith g vs . partialArray bds <$> f vs where
+    (vs, bds) = bind as
+{-# INLINE jacobianWith #-}
+
+-- | 'jacobianWith2 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@.
+--
+-- > jacobian2 == jacobianWith2 (\_ dx -> dx)
+--
+jacobianWith2 :: (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)
+jacobianWith2 g f as = row <$> f vs where
+    (vs, bds) = bind as
+    row a = (primal a, unbindWith g vs (partialArray bds a))
+{-# INLINE jacobianWith2 #-}
 
 diffUU :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a
 diffUU f a = derivative $ f (var a 0)
diff --git a/ad.cabal b/ad.cabal
--- a/ad.cabal
+++ b/ad.cabal
@@ -1,9 +1,10 @@
 Name:         ad
-Version:      0.13
+Version:      0.15
 License:      BSD3
 License-File: LICENSE
 Copyright:    Edward Kmett 2010
-Author:       Edward Kmett 2010
+              Barak Pearlmutter and Jeffrey Mark Siskind 2008-2009
+Author:       Edward Kmett
 Maintainer:   ekmett@gmail.com
 Stability:    Experimental
 Category:     Math
@@ -36,5 +37,4 @@
     Numeric.AD.Internal.Tower
 
 Extra-Source-Files: TODO
-
 GHC-Options: -Wall