diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown
--- a/CHANGELOG.markdown
+++ b/CHANGELOG.markdown
@@ -1,3 +1,9 @@
+4.3
+---
+* Made drastic improvements in the performance of `Tower` and `Sparse` modes thanks to the help of Björn von Sydow.
+* Added constrained convex optimization.
+* Incorporated some suggestions from [herbie](http://herbie.uwplse.org/z) for improving floating point accuracy.
+
 4.2.4
 -----
 * Added `Newton.Double` modules for performance.
diff --git a/README.markdown b/README.markdown
--- a/README.markdown
+++ b/README.markdown
@@ -31,10 +31,10 @@
     Prelude Numeric.AD> diff' (exp . log) 2
     (2.0,1.0)
 
-You can compute the derivative of a function with a constant parameter using `auto` from Numeric.AD.Types:
+You can compute the derivative of a function with a constant parameter using `auto`:
 
-    Prelude Numeric.AD Numeric.AD.Types> let t = 2.0 :: Double
-    Prelude Numeric.AD Numeric.AD.Types> diff (\ x -> auto t * sin x) 0
+    Prelude Numeric.AD> let t = 2.0 :: Double
+    Prelude Numeric.AD> diff (\ x -> auto t * sin x) 0
     2.0
 
 You can use a symbolic numeric type, like the one from `simple-reflect` to obtain symbolic derivatives:
@@ -65,22 +65,22 @@
 
 The answer:
 
-    Prelude Numeric.AD Numeric.AD.Types> headJet $ jet $  grads (\[x,y] -> exp (x * y)) [1,2]
+    Prelude Numeric.AD> headJet $ jet $  grads (\[x,y] -> exp (x * y)) [1,2]
     7.38905609893065
 
 The gradient:
 
-    Prelude Numeric.AD Numeric.AD.Types> headJet $ tailJet $ jet $  grads (\[x,y] -> exp (x * y)) [1,2]
+    Prelude Numeric.AD> headJet $ tailJet $ jet $  grads (\[x,y] -> exp (x * y)) [1,2]
     [14.7781121978613,7.38905609893065]
 
 The hessian (n * n matrix of 2nd derivatives)
 
-    Prelude Numeric.AD Numeric.AD.Types> headJet $ tailJet $ tailJet $ jet $  grads (\[x,y] -> exp (x * y)) [1,2]
+    Prelude Numeric.AD> headJet $ tailJet $ tailJet $ jet $  grads (\[x,y] -> exp (x * y)) [1,2]
     [[29.5562243957226,22.16716829679195],[22.16716829679195,7.38905609893065]]
 
 Or even higher order tensors of derivatives as a jet.
 
-    Prelude Numeric.AD Numeric.AD.Types> headJet $ tailJet $ tailJet $ tailJet $ jet $  grads (\[x,y] -> exp (x * y)) [1,2]
+    Prelude Numeric.AD> headJet $ tailJet $ tailJet $ tailJet $ jet $  grads (\[x,y] -> exp (x * y)) [1,2]
     [[[59.1124487914452,44.3343365935839],[44.3343365935839,14.7781121978613]],[[44.3343365935839,14.7781121978613],[14.7781121978613,7.38905609893065]]]
 
 Note the redundant values caused by the various symmetries in the tensors. The `ad` library is careful to compute
diff --git a/ad.cabal b/ad.cabal
--- a/ad.cabal
+++ b/ad.cabal
@@ -1,5 +1,5 @@
 name:          ad
-version:       4.2.4
+version:       4.3
 license:       BSD3
 license-File:  LICENSE
 copyright:     (c) Edward Kmett 2010-2015,
@@ -75,6 +75,10 @@
   type: git
   location: git://github.com/ekmett/ad.git
 
+flag herbie
+  default: False
+  manual: True
+
 library
   hs-source-dirs: src
   include-dirs: include
@@ -110,6 +114,11 @@
 
   if impl(ghc < 7.8)
     build-depends: tagged >= 0.7 && < 1
+
+  if flag(herbie)
+    build-depends: HerbiePlugin >= 0.1 && < 0.2
+    cpp-options: -DHERBIE
+    ghc-options: -fplugin=Herbie
 
   exposed-modules:
     Numeric.AD
diff --git a/include/instances.h b/include/instances.h
--- a/include/instances.h
+++ b/include/instances.h
@@ -21,7 +21,7 @@
   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))
+  (*)          = mul -- lift2 (*) (\x y -> (y, x))
   negate       = lift1 negate (const (auto (-1)))
   abs          = lift1 abs signum
   signum a     = lift1 signum (const zero) 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
@@ -185,6 +185,9 @@
     (Id dadb, Id dadc) = df (Id a) (Id b) (Id c)
     productRule dbi dci = dadb * dbi + dci * dadc
 
+mul :: (Traversable f, Num a) => Dense f a -> Dense f a -> Dense f a
+mul = lift2 (*) (\x y -> (y, x))
+
 #define BODY1(x)   (Traversable f, x)
 #define BODY2(x,y) (Traversable f, x, y)
 #define HEAD Dense f a
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
@@ -214,3 +214,6 @@
 transposeWith f as = snd . mapAccumL go xss0 where
   go xss b = (tail <$> xss, f b (head <$> xss))
   xss0 = toList <$> as
+
+mul :: (Num a) => Forward a -> Forward a -> Forward a
+mul = lift2 (*) (\x y -> (y, x))
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
@@ -158,6 +158,10 @@
     a = f pb pc
     (dadb, dadc) = df (Id a) (Id pb) (Id pc)
 
+
+mul :: Num a => Kahn a -> Kahn a -> Kahn a
+mul = lift2 (*) (\x y -> (y, x))
+
 #define HEAD Kahn a
 #include <instances.h>
 
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
@@ -199,6 +199,9 @@
     a = f pb pc
     (dadb, dadc) = df (Id a) (Id pb) (Id pc)
 
+mul :: (Reifies s Tape, Num a) => Reverse s a -> Reverse s a -> Reverse s a
+mul = lift2 (*) (\x y -> (y, x))
+
 #define BODY1(x) (Reifies s Tape,x)
 #define BODY2(x,y) (Reifies s Tape,x,y)
 #define HEAD Reverse s a
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
@@ -39,6 +39,9 @@
   , vgrads
   , Grad(..)
   , Grads(..)
+  , terms
+  , deriv
+  , primal
   ) where
 
 import Prelude hiding (lookup)
@@ -48,7 +51,7 @@
 import Control.Comonad.Cofree
 import Control.Monad (join)
 import Data.Data
-import Data.IntMap (IntMap, mapWithKey, unionWith, findWithDefault, toAscList, singleton, insertWith, lookup)
+import Data.IntMap (IntMap, unionWith, findWithDefault, toAscList, singleton, insertWith, lookup)
 import qualified Data.IntMap as IntMap
 import Data.Number.Erf
 import Data.Traversable
@@ -74,13 +77,17 @@
 -- | We only store partials in sorted order, so the map contained in a partial
 -- will only contain partials with equal or greater keys to that of the map in
 -- which it was found. This should be key for efficiently computing sparse hessians.
--- there are only (n + k - 1) choose k distinct nth partial derivatives of a
+-- there are only (n + k - 1) choose (k - 1) distinct nth partial derivatives of a
 -- function with k inputs.
 data Sparse a
   = Sparse !a (IntMap (Sparse a))
   | Zero
   deriving (Show, Data, Typeable)
 
+{-
+
+These functions are now unused.
+
 dropMap :: Int -> IntMap a -> IntMap a
 dropMap n = snd . IntMap.split (n - 1)
 {-# INLINE dropMap #-}
@@ -93,6 +100,7 @@
     (fmap (* b) (dropMap n da))
     (fmap (a *) (dropMap n db))
 {-# INLINE times #-}
+-}
 
 vars :: (Traversable f, Num a) => f a -> f (Sparse a)
 vars = snd . mapAccumL var 0 where
@@ -169,45 +177,42 @@
 a <+> Zero = a
 Sparse a as <+> Sparse b bs = Sparse (a + b) $ unionWith (<+>) as bs
 
+-- The instances for Jacobian for Sparse and Tower are almost identical;
+-- could easily be made exactly equal by small changes.
 instance Num a => Jacobian (Sparse a) where
   type D (Sparse a) = Sparse a
   unary f _ Zero = auto (f 0)
-  unary f dadb (Sparse pb bs) = Sparse (f pb) $ mapWithKey (times dadb) bs
+  unary f dadb (Sparse pb bs) = Sparse (f pb) $ IntMap.map (* dadb) bs
 
   lift1 f _ Zero = auto (f 0)
-  lift1 f df b@(Sparse pb bs) = Sparse (f pb) $ mapWithKey (times (df b)) bs
+  lift1 f df b@(Sparse pb bs) = Sparse (f pb) $ IntMap.map (* df b) bs
 
   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
+    a = Sparse (f pb) $ IntMap.map ((df a b) *) bs
 
   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 _    dadc Zero           (Sparse pc dc) = Sparse (f 0  pc) $ IntMap.map (dadc *) dc
+  binary f dadb _    (Sparse pb db) Zero           = Sparse (f pb 0 ) $ IntMap.map (dadb *) db
   binary f dadb dadc (Sparse pb db) (Sparse pc dc) = Sparse (f pb pc) $
-    unionWith (<+>)
-      (mapWithKey (times dadb) db)
-      (mapWithKey (times dadc) dc)
+    unionWith (<+>)  (IntMap.map (dadb *) db) (IntMap.map (dadc *) dc)
 
   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 Zero c@(Sparse pc dc) = Sparse (f 0 pc) $ IntMap.map (dadc *) dc where dadc = snd (df zero c)
+  lift2 f df b@(Sparse pb db) Zero = Sparse (f pb 0) $ IntMap.map (* 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
     (dadb, dadc) = df b c
-    da = unionWith (<+>)
-      (mapWithKey (times dadb) db)
-      (mapWithKey (times dadc) dc)
+    da = unionWith (<+>) (IntMap.map (dadb *) db) (IntMap.map (dadc *) dc)
 
   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) Zero = a where a = Sparse (f pb 0) (IntMap.map (fst (df a b zero) *) db)
+  lift2_ f df Zero c@(Sparse pc dc) = a where a = Sparse (f 0 pc) (IntMap.map (* snd (df a zero c)) dc)
   lift2_ f df b@(Sparse pb db) c@(Sparse pc dc) = a where
     (dadb, dadc) = df a b c
     a = Sparse (f pb pc) da
-    da = unionWith (<+>)
-      (mapWithKey (times dadb) db)
-      (mapWithKey (times dadc) dc)
+    da = unionWith (<+>) (IntMap.map (dadb *) db) (IntMap.map (dadc *) dc)
 
+
 #define HEAD Sparse a
 #include "instances.h"
 
@@ -254,3 +259,62 @@
 vgrads i = unpacks (unsafeGrads (packs i)) where
   unsafeGrads f as = ds as $ apply f as
 {-# INLINE vgrads #-}
+
+isZero :: Sparse a -> Bool
+isZero Zero = True
+isZero _ = False
+
+-- |
+-- A monomial is used to indicate order of differentiation.
+-- For a k-ary function, it represented as a list of k non-negative Ints.
+-- MI [n_0,n_1...n_{k-1}] denotes differentiation n_0 times with respect
+-- to variable 0, n_1 times to variable 1, etc.
+-- Trailing zeros omitted for efficiency.
+--
+-- Add 1 to variable k (i.e.differentiate once more wrt variable k).
+incMonomial :: Int -> [Int] -> [Int]
+incMonomial k [] = replicate k 0 ++ [1]
+incMonomial 0 (a:as) = a+1:as
+incMonomial k (a:as) = a:incMonomial (k-1) as
+
+-- deriv f mi is the derivative of f of order mi (including higher derivatives).
+deriv :: Sparse a -> [Int] -> Sparse a
+deriv f mi = indx 0 mi f where
+  indx _ [] f = f
+  indx _ _ Zero = Zero
+  indx v (0:as) f = indx (v+1) as f
+  indx v (a:as) (Sparse _ df) = maybe Zero (indx v (a-1 : as)) (lookup v df)
+
+-- The value of the derivative of (f*g) of order mi is
+--       sum [a*primal (deriv f b)*primal (deriv g c) | (a,b,c) <- terms mi ]
+-- It is a bit more complicated in mul' below, since we build the whole tree of
+-- derivatives and want to prune the tree with Zeros as much as possible.
+-- The number of terms in the sum for order MI as of differentiation has
+-- sum (map (+1) as) terms, so this is *much* more efficient
+-- than the naive recursive differentiation with 2^(sum as) terms.
+-- The coefficients a, which collect equivalent derivatives, are suitable products
+-- of binomial coefficients.
+terms :: [Int]-> [(Integer,[Int],[Int])]
+terms [] = [(1,[],[])]
+terms (a:as) = concatMap (f ps) (zip (bins!!a) [0..a]) where
+  ps = terms as
+  bins = iterate next [1]
+  next xs@(_:ts) = 1 : zipWith (+) xs ts ++ [1]
+  next [] = error "impossible"
+  f ps (b,k) = map (\(w,ks,is) -> (w*b,(k:ks),(a-k:is))) ps
+
+mul :: Num a => Sparse a -> Sparse a -> Sparse a
+mul Zero _ = Zero
+mul _ Zero = Zero
+mul f@(Sparse _ am) g@(Sparse _ bm) = Sparse (primal f * primal g) (derivs 0 []) where
+  derivs v mi = IntMap.unions (map fn [v..kMax]) where
+    fn w
+      | and zs = IntMap.empty
+      | otherwise = IntMap.singleton w (Sparse (sum ds) (derivs w mi'))
+      where
+        mi' = incMonomial w mi
+        (zs,ds) = unzip (map derVal (terms mi'))
+        derVal (bin,mif,mig) = (isZero fder || isZero gder, fromIntegral bin * primal fder * primal gder) where
+          fder = deriv f mif
+          gder = deriv g mig
+  kMax = max (maximum (-1:IntMap.keys am)) (maximum (-1:IntMap.keys bm))
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
@@ -34,7 +34,7 @@
   , tower
   ) where
 
-import Prelude hiding (all)
+import Prelude hiding (all, sum)
 #if __GLASGOW_HASKELL__ < 710
 import Control.Applicative hiding ((<**>))
 #endif
@@ -163,6 +163,23 @@
 _        <**> 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 * log xi)) x y
+
+-- mul xs ys = [ sum [xs!!j * ys!!(k-j)*bin k j | j <- [0..k]] | k <- [0..] ]
+-- adapted for efficiency and to handle finite lists xs, ys
+mul:: Num a => Tower a -> Tower a -> Tower a
+mul (Tower []) _ = Tower []
+mul (Tower (a:as)) (Tower bs) = Tower (convs' [1] [a] as bs)
+  where convs' _ _ _ [] = []
+        convs' ps ars as bs = sumProd3 ps ars bs :
+              case as of
+                 [] -> convs'' (next' ps) ars bs
+                 a:as -> convs' (next ps) (a:ars) as bs
+        convs'' _ _ [] = undefined -- convs'' never called with last argument empty
+        convs'' _ _ [_] = []
+        convs'' ps ars (_:bs) = sumProd3 ps ars bs : convs'' (next' ps) ars bs
+        next xs = 1 : zipWith (+) xs (tail xs) ++ [1] -- next row in Pascal's triangle
+        next' xs = zipWith (+) xs (tail xs) ++ [1] -- end part of next row in Pascal's triangle
+        sumProd3 as bs cs = sum (zipWith3 (\x y z -> x*y*z) as bs cs)
 
 #define HEAD Tower a
 #include <instances.h>
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
@@ -1,8 +1,14 @@
-{-# LANGUAGE Rank2Types #-}
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE GADTs #-}
 {-# LANGUAGE BangPatterns #-}
 {-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE DeriveFunctor #-}
+{-# LANGUAGE DeriveFoldable #-}
+{-# LANGUAGE DeriveTraversable #-}
+{-# LANGUAGE ParallelListComp #-}
 -----------------------------------------------------------------------------
 -- |
 -- Copyright   :  (c) Edward Kmett 2010-2015
@@ -21,14 +27,18 @@
   , fixedPoint
   , extremum
   -- * Gradient Ascent/Descent (Reverse AD)
-  , gradientDescent
+  , gradientDescent, constrainedDescent, CC(..), eval
   , gradientAscent
   , conjugateGradientDescent
   , conjugateGradientAscent
   , stochasticGradientDescent
   ) where
 
+#if __GLASGOW_HASKELL__ < 710
+import Data.Foldable (Foldable, all, sum)
+#else
 import Data.Foldable (all, sum)
+#endif
 import Data.Reflection (Reifies)
 import Data.Traversable
 import Numeric.AD.Internal.Combinators
@@ -38,7 +48,7 @@
 import Numeric.AD.Internal.Reverse (Reverse, Tape)
 import Numeric.AD.Internal.Type (AD(..))
 import Numeric.AD.Mode
-import Numeric.AD.Mode.Reverse as Reverse (gradWith, gradWith')
+import Numeric.AD.Mode.Reverse as Reverse (gradWith, gradWith', grad')
 import Numeric.AD.Rank1.Kahn as Kahn (Kahn, grad)
 import qualified Numeric.AD.Rank1.Newton as Rank1
 import Prelude hiding (all, mapM, sum)
@@ -122,6 +132,94 @@
         x1 = fmap (\(xi,gxi) -> xi - eta * gxi) xgx
         (fx1, xgx1) = Reverse.gradWith' (,) f x1
 {-# INLINE gradientDescent #-}
+
+data SEnv (f :: * -> *) a = SEnv { sValue :: a, origEnv :: f a }
+  deriving (Functor, Foldable, Traversable)
+
+-- | Convex constraint, CC, is a GADT wrapper that hides the existential
+-- ('s') which is so prevalent in the rest of the API.  This is an
+-- engineering convenience for managing the skolems.
+data CC f a where
+  CC :: forall f a. (forall s. Reifies s Tape => f (Reverse s a) -> Reverse s a) -> CC f a
+
+-- |@constrainedDescent obj fs env@ optimizes the convex function @obj@
+-- subject to the convex constraints @f <= 0@ where @f `elem` fs@. This is
+-- done using a log barrier to model constraints (i.e. Boyd, Chapter 11.3).
+-- The returned optimal point for the objective function must satisfy @fs@,
+-- but the initial environment, @env@, needn't be feasible.
+constrainedDescent :: forall f a. (Traversable f, RealFloat a, Floating a, Ord a)
+                   => (forall s. Reifies s Tape => f (Reverse s a)
+                                                -> Reverse s a)
+                   -> [CC f a]
+                   -> f a
+                   -> [(a,f a)]
+constrainedDescent objF [] env =
+  map (\x -> (eval objF x, x)) (gradientDescent objF env)
+constrainedDescent objF cs env =
+    let s0       = 1 + maximum [eval c env | CC c <- cs]
+        -- ^ s0 = max ( f_i(0) )
+        cs'      = [CC (\(SEnv sVal rest) -> c rest - sVal) | CC c <- cs]
+        -- ^ f_i' = f_i - s0  and thus f_i' <= 0
+        envS     = SEnv s0 env
+        -- feasible point for f_i', use gd to find feasiblity for f_i
+        cc       = constrainedConvex' (CC sValue) cs' envS ((<=0) . sValue)
+    in case dropWhile ((0 <) . fst) (take (2^(20::Int)) cc) of
+        []                  -> []
+        (_,envFeasible) : _ ->
+            constrainedConvex' (CC objF) cs (origEnv envFeasible) (const True)
+{-# INLINE constrainedDescent #-}
+
+eval :: (Traversable f, Fractional a, Ord a) => (forall s. Reifies s Tape => f (Reverse s a) -> Reverse s a) -> f a -> a
+eval f e = fst (grad' f e)
+{-# INLINE eval #-}
+
+-- | Like 'constrainedDescent' except the initial point must be feasible.
+constrainedConvex' :: forall f a. (Traversable f, RealFloat a, Floating a, Ord a)
+                   => CC f a
+                   -> [CC f a]
+                   -> f a
+                   -> (f a -> Bool)
+                   -> [(a,f a)]
+constrainedConvex' objF cs env term =
+  -- 1. Transform cs using a log barrier with increasing t values.
+  let os   = map (mkOpt objF cs) tValues
+  -- 2. Iteratively run gradientDescent on each os.
+      envs =  [(undefined,env)] :
+              [gD (snd $ last e) o
+                          | o  <- os
+                          | e  <- limEnvs
+                          ]
+      -- Obtain a finite number of elements from the initial len tValues - 1 lists.
+      limEnvs = zipWith id nrSteps envs
+  in dropWhile (not . term . snd) (concat $ drop 1 limEnvs)
+ where
+  tValues = map realToFrac $ take 64 $ iterate (*2) (2 :: a)
+  nrSteps = [take 20 | _ <- [1..length tValues]] ++ [id]
+  -- | `gD f e` is gradient descent with the evaulated result
+  gD e (CC f)  = (eval f e, e) :
+                 map (\x -> (eval f x, x)) (gradientDescent f e)
+{-# INLINE constrainedConvex' #-}
+
+-- @mkOpt u fs t@ converts an inequality convex problem (@u,fs@) into an
+-- unconstrained convex problem using log barrier @u + -(1/t)log(-f_i)@.
+-- As @t@ increases the approximation is more accurate but the gradient
+-- decreases, making the gradient descent more expensive.
+mkOpt :: forall f a. (Traversable f, RealFloat a, Floating a, Ord a)
+      => CC f a -> [CC f a]
+      -> a -> CC f a
+mkOpt (CC o) xs t = CC (\e -> o e + sum (map (\(CC c) -> iHat t c e) xs))
+{-# INLINE mkOpt #-}
+
+iHat :: forall a f. (Traversable f, RealFloat a, Floating a, Ord a)
+     => a
+     -> (forall s. Reifies s Tape => f (Reverse s a) -> Reverse s a)
+     -> (forall s. Reifies s Tape => f (Reverse s a) -> Reverse s a)
+iHat t c e =
+   let r = c e
+   in if r >= 0 || isNaN r
+        then 1  / 0
+        else (-1 / auto t) * log( - (c  e))
+{-# INLINE iHat #-}
 
 -- | The 'stochasticGradientDescent' function approximates
 -- the true gradient of the constFunction by a gradient at
diff --git a/src/Numeric/AD/Rank1/Halley.hs b/src/Numeric/AD/Rank1/Halley.hs
--- a/src/Numeric/AD/Rank1/Halley.hs
+++ b/src/Numeric/AD/Rank1/Halley.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE CPP #-}
 -----------------------------------------------------------------------------
 -- |
 -- Copyright   :  (c) Edward Kmett 2010-2015
@@ -49,7 +50,12 @@
 findZero f = go where
   go x = x : if x == xn then [] else go xn where
     (y:y':y'':_) = diffs0 f x
-    xn = x - 2*y*y'/(2*y'*y'-y*y'')
+    xn = x - 2*y*y'/(2*y'*y'-y*y'') -- 9.606671960457536 bits error
+       -- = x - recip (y'/y - y''/ y') -- "improved error" = 6.640625e-2 bits
+       -- = x - y' / (y'/y/y' - y''/2) -- "improved error" = 1.4
+#ifdef HERBIE
+{-# ANN findZero "NoHerbie" #-}
+#endif
 {-# INLINE findZero #-}
 
 -- | The 'inverse' function inverts a scalar function using