diff --git a/Numeric/AD/Internal.hs b/Numeric/AD/Internal.hs
--- a/Numeric/AD/Internal.hs
+++ b/Numeric/AD/Internal.hs
@@ -1,4 +1,4 @@
-{-# LANGUAGE GeneralizedNewtypeDeriving, TemplateHaskell, FlexibleContexts, FlexibleInstances, MultiParamTypeClasses #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving, TemplateHaskell, FlexibleContexts, FlexibleInstances, MultiParamTypeClasses, DeriveFunctor, DeriveFoldable, DeriveTraversable #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.AD.Internal
@@ -18,6 +18,7 @@
     , unprobe
     , probed
     , unprobed
+    , Pair(..)
     ) where
 
 import Control.Applicative
@@ -26,6 +27,8 @@
 import Data.Monoid
 import Data.Traversable (Traversable, mapAccumL)
 import Data.Foldable (Foldable, toList)
+
+data Pair a b = Pair a b deriving (Eq, Ord, Show, Read, Functor, Foldable, Traversable)
 
 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)
diff --git a/Numeric/AD/Newton.hs b/Numeric/AD/Newton.hs
--- a/Numeric/AD/Newton.hs
+++ b/Numeric/AD/Newton.hs
@@ -1,4 +1,4 @@
-{-# LANGUAGE Rank2Types, BangPatterns #-}
+{-# LANGUAGE Rank2Types, BangPatterns, ScopedTypeVariables #-}
 -----------------------------------------------------------------------------
 -- |
 -- Module      :  Numeric.AD.Newton
@@ -14,23 +14,32 @@
     (
     -- * Newton's Method (Forward AD)
       findZero
+    , findZeroM
     , inverse
+    , inverseM
     , fixedPoint
+    , fixedPointM
     , extremum
-    -- * Gradient Descent (Reverse AD)
+    , extremumM
+    -- * Gradient Ascent/Descent (Reverse AD)
     , gradientDescent
+    , gradientDescentM
+    , gradientAscent
+    , gradientAscentM
     -- * Exposed Types
     , AD(..)
     , Mode(..)
     ) where
 
 import Prelude hiding (all)
+import Control.Monad (liftM)
+import Data.MList
 import Numeric.AD.Classes
 import Numeric.AD.Internal
 import Data.Foldable (all)
 import Data.Traversable (Traversable)
-import Numeric.AD.Forward (diff, diff')
-import Numeric.AD.Reverse (gradWith')
+import Numeric.AD.Forward (diff, diff', diffM, diffM')
+import Numeric.AD.Reverse (gradWith', gradWithM')
 import Numeric.AD.Internal.Composition
 
 -- | The 'findZero' function finds a zero of a scalar function using
@@ -45,35 +54,65 @@
 --  > take 10 $ findZero ((+1).(^2)) (1 :+ 1)  -- converge to (0 :+ 1)@
 --
 findZero :: Fractional a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
-findZero f x0 = iterate (\x -> let (y,y') = diff' f x in x - y/y') x0
+findZero f = go
+    where
+        go x = x : go (x - y/y') 
+            where
+                (y,y') = diff' f x
 {-# INLINE findZero #-}
 
+findZeroM :: (Monad m, Fractional a) => (forall s. Mode s => AD s a -> m (AD s a)) -> a -> MList m a
+findZeroM f x0 = MList (go x0)
+    where
+        go x = return $ 
+               MCons x $ 
+               MList $ do
+                (y,y') <- diffM' f x
+                go (x - y/y')
+{-# INLINE findZeroM #-}
+
 -- | The 'inverseNewton' function inverts a scalar function using
 -- Newton's method; its output is a stream of increasingly accurate
 -- results.  (Modulo the usual caveats.)
 --
 -- Example:
 --
--- > take 10 $ inverseNewton sqrt 1 (sqrt 10)  -- converge to 10
+-- > take 10 $ inverseNewton sqrt 1 (sqrt 10)  -- converges to 10
 --
 inverse :: Fractional a => (forall s. Mode s => AD s a -> AD s a) -> a -> a -> [a]
 inverse f x0 y = findZero (\x -> f x - lift y) x0
 {-# INLINE inverse  #-}
 
+inverseM :: (Monad m, Fractional a) => (forall s. Mode s => AD s a -> m (AD s a)) -> a -> a -> MList m a
+inverseM f x0 y = findZeroM (\x -> subtract (lift y) `liftM` f x) x0
+{-# INLINE inverseM  #-}
+
 -- | The 'fixedPoint' function find a fixedpoint of a scalar
 -- function using Newton's method; its output is a stream of
 -- increasingly accurate results.  (Modulo the usual caveats.)
+-- 
+-- > take 10 $ fixedPoint cos 1 -- converges to 0.7390851332151607
 fixedPoint :: Fractional a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
 fixedPoint f = findZero (\x -> f x - x)
 {-# INLINE fixedPoint #-}
 
+fixedPointM :: (Monad m, Fractional a) => (forall s. Mode s => AD s a -> m (AD s a)) -> a -> MList m a
+fixedPointM f = findZeroM (\x -> subtract x `liftM` f x)
+{-# INLINE fixedPointM #-}
+
 -- | The 'extremum' function finds an extremum of a scalar
 -- function using Newton's method; produces a stream of increasingly
 -- accurate results.  (Modulo the usual caveats.)
+--
+-- > take 10 $ extremum cos 1 -- convert to 0 
 extremum :: Fractional a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
-extremum f x0 = findZero (diff (decompose . f . compose)) x0
+extremum f = findZero (diff (decompose . f . compose))
 {-# INLINE extremum #-}
 
+extremumM :: (Monad m, Fractional a) => (forall s. Mode s => AD s a -> m (AD s a)) -> a -> MList m a
+extremumM f = findZeroM (diffM (liftM decompose . f . compose))
+{-# INLINE extremumM #-}
+
 -- | The 'gradientDescent' function performs a multivariate
 -- optimization, based on the naive-gradient-descent in the file
 -- @stalingrad\/examples\/flow-tests\/pre-saddle-1a.vlad@ from the
@@ -96,5 +135,36 @@
                 zeroGrad = all (\(_,g) -> g == 0)
                 x1 = fmap (\(xi,gxi) -> xi - eta * gxi) xgx
                 (fx1, xgx1) = gradWith' (,) f x1
-
 {-# INLINE gradientDescent #-}
+
+gradientAscent :: (Traversable f, Fractional a, Ord a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> [f a]
+gradientAscent f = gradientDescent (negate . f)
+{-# INLINE gradientAscent #-}
+
+-- monadic gradient descent
+gradientDescentM :: (Traversable f, Monad m, Fractional a, Ord a) => (forall s. Mode s => f (AD s a) -> m (AD s a)) -> f a -> MList m (f a)
+gradientDescentM f x0 = MList $ do
+        (fx0, xgx0) <- gradWithM' (,) f x0
+        go x0 fx0 xgx0 0.1 (0 :: Int)
+    where
+        go x fx xgx !eta !i
+            | eta == 0  = return MNil -- step size is 0
+            | otherwise = do
+                (fx1, xgx1) <- gradWithM' (,) f x1
+                case () of
+                 _ | fx1 > fx     -> go x fx xgx (eta/2) 0 -- we stepped too far
+                   | zeroGrad xgx -> return MNil -- gradient is 0
+                   | otherwise    -> return $ 
+                                     MCons x1 $ 
+                                     MList $
+                                     if i == 10
+                                     then go x1 fx1 xgx1 (eta*2) 0
+                                     else go x1 fx1 xgx1 eta (i+1)
+            where
+                x1 = fmap (\(xi,gxi) -> xi - eta * gxi) xgx
+                zeroGrad = all (\(_,g) -> g == 0)
+{-# INLINE gradientDescentM #-}
+
+gradientAscentM :: (Traversable f, Monad m, Fractional a, Ord a) => (forall s. Mode s => f (AD s a) -> m (AD s a)) -> f a -> MList m (f a)
+gradientAscentM f = gradientDescentM (liftM negate . f)
+{-# INLINE gradientAscentM #-}
diff --git a/ad.cabal b/ad.cabal
--- a/ad.cabal
+++ b/ad.cabal
@@ -1,5 +1,5 @@
 Name:         ad
-Version:      0.19
+Version:      0.20
 License:      BSD3
 License-File: LICENSE
 Copyright:    Edward Kmett 2010
@@ -21,6 +21,7 @@
     data-reify >= 0.5 && < 0.6, 
     containers >= 0.2 && < 0.4,
     template-haskell >= 2.4 && < 2.5,
+    mlist >= 0.0.2 && <= 0.1,
     array >= 0.2 && < 0.4
 
 Exposed-Modules: