diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown
--- a/CHANGELOG.markdown
+++ b/CHANGELOG.markdown
@@ -1,3 +1,9 @@
+4.3.1
+-----
+* Further improvements have been made in the performance of `Sparse` mode, at least asymptotically, when used on functions with many variables.
+  Since this is the target use-case for `Sparse` in the first place, this seems like a good trade-off. Note: this results in an API change, but
+  only in the API of an `Internal` module, so this is treated as a minor version bump.
+
 4.3
 ---
 * Made drastic improvements in the performance of `Tower` and `Sparse` modes thanks to the help of Björn von Sydow.
diff --git a/ad.cabal b/ad.cabal
--- a/ad.cabal
+++ b/ad.cabal
@@ -1,5 +1,5 @@
 name:          ad
-version:       4.3
+version:       4.3.1
 license:       BSD3
 license-File:  LICENSE
 copyright:     (c) Edward Kmett 2010-2015,
@@ -12,7 +12,7 @@
 bug-reports:   http://github.com/ekmett/ad/issues
 build-type:    Custom
 cabal-version: >= 1.10
-tested-with:   GHC==7.0.1, GHC == 7.0.4, GHC == 7.2.2, GHC == 7.4.2, GHC == 7.6.3, GHC == 7.8.4, GHC == 7.10.1
+tested-with:   GHC==7.0.1, GHC == 7.0.4, GHC == 7.2.2, GHC == 7.4.2, GHC == 7.6.3, GHC == 7.8.4, GHC == 7.10.1, GHC == 7.10.2
 synopsis:      Automatic Differentiation
 extra-source-files:
   .ghci
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
@@ -23,9 +23,9 @@
 -- Handle with care.
 -----------------------------------------------------------------------------
 module Numeric.AD.Internal.Sparse
-  ( Index(..)
-  , emptyIndex
-  , addToIndex
+  ( Monomial(..)
+  , emptyMonomial
+  , addToMonomial
   , indices
   , Sparse(..)
   , apply
@@ -40,7 +40,6 @@
   , Grad(..)
   , Grads(..)
   , terms
-  , deriv
   , primal
   ) where
 
@@ -60,18 +59,18 @@
 import Numeric.AD.Jacobian
 import Numeric.AD.Mode
 
-newtype Index = Index (IntMap Int)
+newtype Monomial = Monomial (IntMap Int)
 
-emptyIndex :: Index
-emptyIndex = Index IntMap.empty
-{-# INLINE emptyIndex #-}
+emptyMonomial :: Monomial
+emptyMonomial = Monomial IntMap.empty
+{-# INLINE emptyMonomial #-}
 
-addToIndex :: Int -> Index -> Index
-addToIndex k (Index m) = Index (insertWith (+) k 1 m)
-{-# INLINE addToIndex #-}
+addToMonomial :: Int -> Monomial -> Monomial
+addToMonomial k (Monomial m) = Monomial (insertWith (+) k 1 m)
+{-# INLINE addToMonomial #-}
 
-indices :: Index -> [Int]
-indices (Index as) = uncurry (flip replicate) `concatMap` toAscList as
+indices :: Monomial -> [Int]
+indices (Monomial as) = uncurry (flip replicate) `concatMap` toAscList as
 {-# INLINE indices #-}
 
 -- | We only store partials in sorted order, so the map contained in a partial
@@ -84,24 +83,6 @@
   | 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 #-}
-
-times :: Num a => Sparse a -> Int -> Sparse a -> Sparse a
-times Zero _ _ = Zero
-times _ _ Zero = Zero
-times a@(Sparse pa da) n b@(Sparse pb db) = Sparse (pa * pb) $
-  unionWith (+)
-    (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
   var !n a = (n + 1, Sparse a $ singleton n $ auto 1)
@@ -127,13 +108,19 @@
 
 ds :: (Traversable f, Num a) => f b -> Sparse a -> Cofree f a
 ds fs Zero = r where r = 0 :< (r <$ fs)
-ds fs (as@(Sparse a _)) = a :< (go emptyIndex <$> fns) where
+ds fs (as@(Sparse a _)) = a :< (go emptyMonomial <$> fns) where
   fns = skeleton fs
-  -- go :: Index -> Int -> Cofree f a
+  -- go :: Monomial -> Int -> Cofree f a
   go ix i = partial (indices ix') as :< (go ix' <$> fns) where
-    ix' = addToIndex i ix
+    ix' = addToMonomial i ix
 {-# INLINE ds #-}
 
+partialS :: Num a => [Int] -> Sparse a -> Sparse a
+partialS []     a             = a
+partialS (n:ns) (Sparse _ da) = partialS ns $ findWithDefault Zero n da
+partialS _      Zero          = Zero
+{-# INLINE partialS #-}
+
 partial :: Num a => [Int] -> Sparse a -> a
 partial []     (Sparse a _)  = a
 partial (n:ns) (Sparse _ da) = partial ns $ findWithDefault (auto 0) n da
@@ -265,56 +252,40 @@
 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
+--
+-- @
+-- 'sum' [a * 'primal' ('partialS' ('indices' b) f) * 'primal' ('partialS' ('indices' c) g) | (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 'Zero's 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
+terms :: Monomial -> [(Integer,Monomial,Monomial)]
+terms (Monomial m) = t (toAscList m) where
+  t [] = [(1,emptyMonomial,emptyMonomial)]
+  t ((k,a):ts) = concatMap (f (t ts)) (zip (bins!!a) [0..a]) where
+    f ps (b,i) = map (\(w,Monomial mf,Monomial mg) -> (w*b,Monomial (IntMap.insert k i mf), Monomial (IntMap.insert k (a-i) mg))) ps
   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
+mul f@(Sparse _ am) g@(Sparse _ bm) = Sparse (primal f * primal g) (derivs 0 emptyMonomial) 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
+        mi' = addToMonomial 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))
+          fder = partialS (indices mif) f
+          gder = partialS (indices mig) g
+  kMax = maybe (-1) (fst.fst) (IntMap.maxViewWithKey am) `max` maybe (-1) (fst.fst) (IntMap.maxViewWithKey bm)