diff --git a/ad.cabal b/ad.cabal
--- a/ad.cabal
+++ b/ad.cabal
@@ -1,5 +1,5 @@
 name:         ad
-version:      4.2
+version:      4.2.0.1
 license:      BSD3
 license-File: LICENSE
 copyright:    (c) Edward Kmett 2010-2014,
@@ -142,10 +142,11 @@
     Numeric.AD.Newton
     Numeric.AD.Rank1.Forward
     Numeric.AD.Rank1.Forward.Double
-    Numeric.AD.Rank1.Sparse
-    Numeric.AD.Rank1.Tower
+    Numeric.AD.Rank1.Halley
     Numeric.AD.Rank1.Kahn
     Numeric.AD.Rank1.Newton
+    Numeric.AD.Rank1.Sparse
+    Numeric.AD.Rank1.Tower
 
   other-modules:
     Numeric.AD.Internal.Combinators
diff --git a/bench/BlackScholes.hs b/bench/BlackScholes.hs
--- a/bench/BlackScholes.hs
+++ b/bench/BlackScholes.hs
@@ -3,6 +3,7 @@
 import Data.Number.Erf
 import qualified Numeric.AD as Mixed
 import qualified Numeric.AD.Mode.Forward as Forward
+import qualified Numeric.AD.Mode.Forward.Double as ForwardDouble
 import qualified Numeric.AD.Mode.Kahn as Kahn
 import qualified Numeric.AD.Mode.Reverse as Reverse
 import qualified Numeric.AD.Mode.Sparse as Sparse
@@ -39,6 +40,9 @@
         [ bench "greeks Double" $ nf (runDouble $ \r s v t k -> Forward.jacobian (fromPair . bs) [r, s, v, t, k]) 2
         , bench "greeks Float" $ nf (runFloat $ \r s v t k -> Forward.jacobian (fromPair . bs) [r, s, v, t, k]) 2
         ]
+    , bgroup "ForwardDouble"
+        [ bench "greeks Double" $ nf (runDouble $ \r s v t k -> ForwardDouble.jacobian (fromPair . bs) [r, s, v, t, k]) 2
+        ]
     , bgroup "Kahn"
         [ bench "greeks Double" $ nf (runDouble $ \r s v t k -> Kahn.jacobian (fromPair . bs) [r, s, v, t, k]) 2
         , bench "higherGreeks Double" $ nf (runDouble $ \r s v t k -> (Kahn.hessian (fst . bs) [r, s, v, t, k], Kahn.hessian (snd . bs) [r, s, v, t, k])) 2
@@ -63,12 +67,12 @@
         , bench "higherGreeks Float" $ nf (runFloat $ \r s v t k -> (Sparse.hessian (fst . bs) [r, s, v, t, k], Sparse.hessian (snd . bs) [r, s, v, t, k])) 2
         , bench "highererGreeks Float" $ nf (runFloat $ \r s v t k -> Sparse.hessianF (fromPair . bs) [r, s, v, t, k]) 2
         ]
---    , bgroup "Mixed"
---        [ bench "greeks Double" $ nf (runDouble $ \r s v t k -> Mixed.jacobian (fromPair . bs) [r, s, v, t, k]) 2
---        , bench "higherGreeks Double" $ nf (runDouble $ \r s v t k -> (Mixed.hessian (fst . bs) [r, s, v, t, k], Mixed.hessian (snd . bs) [r, s, v, t, k])) 2
---        , bench "highererGreeks Double" $ nf (runDouble $ \r s v t k -> Mixed.hessianF (fromPair . bs) [r, s, v, t, k]) 2
---        , bench "greeks Float" $ nf (runFloat $ \r s v t k -> Mixed.jacobian (fromPair . bs) [r, s, v, t, k]) 2
---        , bench "higherGreeks Float" $ nf (runFloat $ \r s v t k -> (Mixed.hessian (fst . bs) [r, s, v, t, k], Mixed.hessian (snd . bs) [r, s, v, t, k])) 2
---        , bench "highererGreeks Float" $ nf (runFloat $ \r s v t k -> Mixed.hessianF (fromPair . bs) [r, s, v, t, k]) 2
---        ]
+   , bgroup "Mixed"
+       [ bench "greeks Double" $ nf (runDouble $ \r s v t k -> Mixed.jacobian (fromPair . bs) [r, s, v, t, k]) 2
+       , bench "higherGreeks Double" $ nf (runDouble $ \r s v t k -> (Mixed.hessian (fst . bs) [r, s, v, t, k], Mixed.hessian (snd . bs) [r, s, v, t, k])) 2
+       , bench "highererGreeks Double" $ nf (runDouble $ \r s v t k -> Mixed.hessianF (fromPair . bs) [r, s, v, t, k]) 2
+       , bench "greeks Float" $ nf (runFloat $ \r s v t k -> Mixed.jacobian (fromPair . bs) [r, s, v, t, k]) 2
+       , bench "higherGreeks Float" $ nf (runFloat $ \r s v t k -> (Mixed.hessian (fst . bs) [r, s, v, t, k], Mixed.hessian (snd . bs) [r, s, v, t, k])) 2
+       , bench "highererGreeks Float" $ nf (runFloat $ \r s v t k -> Mixed.hessianF (fromPair . bs) [r, s, v, t, k]) 2
+       ]
     ]
diff --git a/src/Numeric/AD/Rank1/Halley.hs b/src/Numeric/AD/Rank1/Halley.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/AD/Rank1/Halley.hs
@@ -0,0 +1,89 @@
+-----------------------------------------------------------------------------
+-- |
+-- Copyright   :  (c) Edward Kmett 2010-2014
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only
+--
+-- Root finding using Halley's rational method (the second in
+-- the class of Householder methods). Assumes the function is three
+-- times continuously differentiable and converges cubically when
+-- progress can be made.
+--
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Rank1.Halley
+  (
+  -- * Halley's Method (Tower AD)
+    findZero
+  , inverse
+  , fixedPoint
+  , extremum
+  ) where
+
+import Prelude hiding (all)
+import Numeric.AD.Internal.Forward (Forward)
+import Numeric.AD.Internal.On
+import Numeric.AD.Internal.Tower (Tower)
+import Numeric.AD.Mode
+import Numeric.AD.Rank1.Tower (diffs0)
+import Numeric.AD.Rank1.Forward (diff)
+
+-- $setup
+-- >>> import Data.Complex
+
+-- | The 'findZero' function finds a zero of a scalar function using
+-- Halley's method; its output is a stream of increasingly accurate
+-- results.  (Modulo the usual caveats.) If the stream becomes constant
+-- ("it converges"), no further elements are returned.
+--
+-- Examples:
+--
+-- >>> take 10 $ findZero (\x->x^2-4) 1
+-- [1.0,1.8571428571428572,1.9997967892704736,1.9999999999994755,2.0]
+--
+-- >>> last $ take 10 $ findZero ((+1).(^2)) (1 :+ 1)
+-- 0.0 :+ 1.0
+findZero :: (Fractional a, Eq a) => (Tower a -> Tower a) -> a -> [a]
+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'')
+{-# INLINE findZero #-}
+
+-- | The 'inverse' function inverts a scalar function using
+-- Halley's method; its output is a stream of increasingly accurate
+-- results.  (Modulo the usual caveats.) If the stream becomes constant
+-- ("it converges"), no further elements are returned.
+--
+-- 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) => (Tower a -> Tower a) -> a -> a -> [a]
+inverse f x0 y = findZero (\x -> f x - auto y) x0
+{-# INLINE inverse  #-}
+
+-- | The 'fixedPoint' function find a fixedpoint of a scalar
+-- function using Halley's method; its output is a stream of
+-- increasingly accurate results.  (Modulo the usual caveats.)
+--
+-- If the stream becomes constant ("it converges"), no further
+-- elements are returned.
+--
+-- >>> last $ take 10 $ fixedPoint cos 1
+-- 0.7390851332151607
+fixedPoint :: (Fractional a, Eq a) => (Tower a -> Tower a) -> a -> [a]
+fixedPoint f = findZero (\x -> f x - x)
+{-# INLINE fixedPoint #-}
+
+
+-- | The 'extremum' function finds an extremum of a scalar
+-- function using Halley's method; produces a stream of increasingly
+-- accurate results.  (Modulo the usual caveats.) If the stream becomes
+-- constant ("it converges"), no further elements are returned.
+--
+-- >>> take 10 $ extremum cos 1
+-- [1.0,0.29616942658570555,4.59979519460002e-3,1.6220740159042513e-8,0.0]
+extremum :: (Fractional a, Eq a) => (On (Forward (Tower a)) -> On (Forward (Tower a))) -> a -> [a]
+extremum f = findZero (diff (off . f . On))
+{-# INLINE extremum #-}