diff --git a/src/Data/AdditiveGroup.hs b/src/Data/AdditiveGroup.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/AdditiveGroup.hs
@@ -0,0 +1,71 @@
+----------------------------------------------------------------------
+-- |
+-- Module      :   Data.AdditiveGroup
+-- Copyright   :  (c) Conal Elliott and Andy J Gill 2008
+-- License     :  BSD3
+-- 
+-- Maintainer  :  conal@conal.net, andygill@ku.edu
+-- Stability   :  experimental
+-- 
+-- Groups: zero, addition, and negation (additive inverse)
+----------------------------------------------------------------------
+
+module Data.AdditiveGroup
+  ( 
+    AdditiveGroup(..), (^-^)
+  ) where
+
+import Control.Applicative
+import Data.Complex hiding (magnitude)
+
+infixl 6 ^+^, ^-^
+
+-- | Additive group @v@.
+class AdditiveGroup v where
+  -- | The zero element: identity for '(^+^)'
+  zeroV :: v
+  -- | Add vectors
+  (^+^) :: v -> v -> v
+  -- | Additive inverse
+  negateV :: v -> v
+
+-- | Group subtraction
+(^-^) :: AdditiveGroup v => v -> v -> v
+v ^-^ v' = v ^+^ negateV v'
+
+instance AdditiveGroup Double where
+  zeroV   = 0.0
+  (^+^)   = (+)
+  negateV = negate
+
+instance AdditiveGroup Float where
+  zeroV   = 0.0
+  (^+^)   = (+)
+  negateV = negate
+
+instance (RealFloat v, AdditiveGroup v) => AdditiveGroup (Complex v) where
+  zeroV   = zeroV :+ zeroV
+  (^+^)   = (+)
+  negateV = negate
+
+-- Hm.  The 'RealFloat' constraint is unfortunate here.  It's due to a
+-- questionable decision to place 'RealFloat' into the definition of the
+-- 'Complex' /type/, rather than in functions and instances as needed.
+
+instance (AdditiveGroup u,AdditiveGroup v) => AdditiveGroup (u,v) where
+  zeroV             = (zeroV,zeroV)
+  (u,v) ^+^ (u',v') = (u^+^u',v^+^v')
+  negateV (u,v)     = (negateV u, negateV v)
+
+instance (AdditiveGroup u,AdditiveGroup v,AdditiveGroup w)
+    => AdditiveGroup (u,v,w) where
+  zeroV                  = (zeroV,zeroV,zeroV)
+  (u,v,w) ^+^ (u',v',w') = (u^+^u',v^+^v',w^+^w')
+  negateV (u,v,w)        = (negateV u, negateV v, negateV w)
+
+
+-- Standard instance for an applicative functor applied to a vector space.
+instance AdditiveGroup v => AdditiveGroup (a->v) where
+  zeroV   = pure   zeroV
+  (^+^)   = liftA2 (^+^)
+  negateV = fmap   negateV
diff --git a/src/Data/Cross.hs b/src/Data/Cross.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Cross.hs
@@ -0,0 +1,110 @@
+{-# LANGUAGE FlexibleInstances, TypeOperators, UndecidableInstances
+           , TypeSynonymInstances
+  #-}
+{-# OPTIONS_GHC -Wall #-}
+----------------------------------------------------------------------
+-- |
+-- Module      :  Data.Cross
+-- Copyright   :  (c) Conal Elliott 2008
+-- License     :  BSD3
+-- 
+-- Maintainer  :  conal@conal.net
+-- Stability   :  experimental
+-- 
+-- Cross products and normals
+----------------------------------------------------------------------
+
+module Data.Cross
+  (
+    HasNormal(..), normal
+  , One, Two, Three
+  , HasCross2(..), HasCross3(..)
+  ) where
+
+import Data.VectorSpace
+import Data.LinearMap
+import Data.Derivative
+
+-- | Thing with a normal vector (not necessarily normalized).
+class HasNormal v where normalVec :: v -> v
+
+-- | Normalized normal vector.  See also 'cross'.
+normal :: (HasNormal v, InnerSpace v s, Floating s) => v -> v
+normal = normalized . normalVec
+
+
+-- | Singleton
+type One   s = s
+
+-- | Homogeneous pair
+type Two   s = (s,s)
+
+-- | Homogeneous triple
+type Three s = (s,s,s)
+
+-- | Cross product of various forms of 2D vectors
+class HasCross2 v where cross2 :: v -> v
+
+instance Num s => HasCross2 (s,s) where
+  cross2 (x,y) = (-y,x)  -- or @(y,-x)@?
+
+-- TODO: Eliminate the 'Num' constraint by using negateV.
+
+-- "Variable occurs more often in a constraint than in the instance
+-- head".  Hence UndecidableInstances.
+
+instance (LMapDom a s, VectorSpace v s, HasCross2 v) => HasCross2 (a:>v) where
+  -- 2d cross-product is linear
+  cross2 = fmapD cross2
+
+instance (Num s, LMapDom s s) => HasNormal (One s :> Two s) where
+  normalVec v = cross2 (derivativeAt v 1)
+
+-- Does this problem come from the choice of 'VectorSpace' instance?
+
+instance (Num s, LMapDom s s, VectorSpace s s)
+    => HasNormal (Two (One s :> s)) where
+  normalVec = unpairD . normalVec . pairD
+
+
+-- | Cross product of various forms of 3D vectors
+class HasCross3 v where cross3 :: v -> v -> v
+
+instance Num s => HasCross3 (s,s,s) where
+  (ax,ay,az) `cross3` (bx,by,bz) = ( ay * bz - az * by
+                                   , az * bx - ax * bz
+                                   , ax * by - ay * bx )
+
+-- TODO: Eliminate the 'Num' constraint by using 'VectorSpace' operations.
+
+instance (LMapDom a s, VectorSpace v s, HasCross3 v) => HasCross3 (a:>v) where
+  -- 3D cross-product is bilinear (curried linear)
+  cross3 = distrib cross3
+
+instance (Num s, LMapDom s s) => HasNormal (Two s :> Three s) where
+  normalVec v = d (1,0) `cross3` d (0,1)
+   where
+     d = derivativeAt v
+
+instance (Num s, VectorSpace s s, LMapDom s s) => HasNormal (Three (Two s :> s)) where
+  normalVec = untripleD . normalVec . tripleD
+
+---- Could go elsewhere
+
+pairD :: (LMapDom a s, VectorSpace b s, VectorSpace c s) =>
+         (a:>b,a:>c) -> a:>(b,c)
+pairD (u,v) = liftD2 (,) u v
+
+tripleD :: (LMapDom a s, VectorSpace b s, VectorSpace c s, VectorSpace d s) =>
+           (a:>b,a:>c,a:>d) -> a:>(b,c,d)
+tripleD (u,v,w) = liftD3 (,,) u v w
+
+unpairD :: (LMapDom a s, VectorSpace a s, VectorSpace b s, VectorSpace c s) =>
+           (a :> (b,c)) -> (a:>b, a:>c)
+unpairD d = (fst <$>> d, snd <$>> d)
+
+untripleD :: ( LMapDom a s , VectorSpace a s, VectorSpace b s
+             , VectorSpace c s, VectorSpace d s) =>
+             (a :> (b,c,d)) -> (a:>b, a:>c, a:>d)
+untripleD d =
+  ((\ (a,_,_) -> a) <$>> d, (\ (_,b,_) -> b) <$>> d, (\ (_,_,c) -> c) <$>> d)
diff --git a/src/Data/Derivative.hs b/src/Data/Derivative.hs
--- a/src/Data/Derivative.hs
+++ b/src/Data/Derivative.hs
@@ -1,5 +1,12 @@
-{-# LANGUAGE TypeOperators, MultiParamTypeClasses, UndecidableInstances #-}
-{-# OPTIONS_GHC -Wall #-}
+{-# LANGUAGE TypeOperators, MultiParamTypeClasses, UndecidableInstances
+           , TypeSynonymInstances, FlexibleInstances, FunctionalDependencies
+           , FlexibleContexts
+           , GeneralizedNewtypeDeriving, StandaloneDeriving
+  #-}
+
+-- TODO: remove FlexibleContexts
+
+-- {-# OPTIONS_GHC -Wall #-}
 ----------------------------------------------------------------------
 -- |
 -- Module      :  Data.Derivative
@@ -9,171 +16,178 @@
 -- Maintainer  :  conal@conal.net
 -- Stability   :  experimental
 -- 
--- Infinite derivative towers via linear maps.  See blog posts
--- <http://conal.net/blog/tag/derivatives/>
+-- This module is a wrapper around Data.Maclaurin or Data.Horner, to
+-- change the 'VectorSpace' instance for '(:>)'.
 ----------------------------------------------------------------------
 
 module Data.Derivative
   (
-    (:>)(..), (:~>), dZero, dConst
+    (:>), powVal, derivative, derivativeAt
+  , (:~>), dZero, pureD
+  , fmapD, (<$>>){-, (<*>>)-}, liftD2, liftD3
   , idD, fstD, sndD
-  , linearD, distribD
+  , linearD, distrib
   , (@.), (>-<)
   ) where
 
-import Control.Applicative
+import Data.Function (on)
 
 import Data.VectorSpace
 import Data.NumInstances ()
+import Data.LinearMap
 
+import qualified Data.Maclaurin as D
+-- import qualified Data.Horner as D
 
-infixr 9 `D`, @.
+-- TODO: sync Data.Horner interface with Data.Maclaurin.  Better: refactor
+-- Maclaurin into a module that depends on the representation (e.g.,
+-- Horner vs Maclaurin) and the rest that doesn't.
+
+infixr 9 @.
+infixl 4 {-<*>>,-} <$>>
 infix  0 >-<
 
 
 -- | Tower of derivatives.
--- 
--- Warning, the 'Applicative' instance is missing its 'pure' (due to a
--- 'VectorSpace' type constraint).  Use 'dConst' instead.
-data a :> b = D { dVal :: b, dDeriv :: a :-* (a :> b) }
+newtype (a :> b) = T { unT :: a D.:> b }
 
--- data a :> b = D b (a :-* (a :> b))
+deriving instance (VectorSpace b b, LMapDom a b, Num b       ) => Num           (a :> b)
+deriving instance (VectorSpace b b, LMapDom a b, Fractional b) => Fractional    (a :> b)
+deriving instance (VectorSpace b b, LMapDom a b, Floating b  ) => Floating      (a :> b)
+deriving instance (VectorSpace b s, LMapDom a s              ) => AdditiveGroup (a :> b)
 
+
+inT :: ((a D.:> b) -> (c D.:> d))
+    -> ((a   :> b) -> (c   :> d))
+inT  = (T .).(. unT)
+
+inT2 :: ((a D.:> b) -> (c D.:> d) -> (e D.:> f))
+     -> ((a   :> b) -> (c   :> d) -> (e   :> f))
+inT2  = (inT .).(. unT)
+
+inT3 :: ((a D.:> b) -> (c D.:> d) -> (e D.:> f) -> (g D.:> h))
+     -> ((a   :> b) -> (c   :> d) -> (e   :> f) -> (g   :> h))
+inT3  = (inT2 .).(. unT)
+
+-- | Extract the value from a derivative tower
+powVal :: (a :> b) -> b
+powVal = D.powVal . unT
+
+-- | Extract the derivative from a derivative tower
+derivative :: (VectorSpace b s, LMapDom a s) =>
+              (a :> b) -> (a :-* (a :> b))
+derivative = fmapL T . D.derivative . unT
+
+-- | Sampled derivative.  For avoiding an awkward typing problem related
+-- to the two required 'VectorSpace' instances.
+derivativeAt :: (LMapDom a s, VectorSpace b s) =>
+                (a :> b) -> a -> (a :> b)
+derivativeAt d = T . D.derivativeAt (unT d)
+
+-- The definition of 'derivativeAt' takes care to share partial
+-- applications of 'D.derivativeAt', which is useful in power series
+-- representations for which 'derivative' is not free (Horner).
+
 -- | Infinitely differentiable functions
 type a :~> b = a -> (a:>b)
 
--- So we could define
--- 
---   data a :> b = D b (a :~> b)
--- 
--- with the restriction that the a :~> b is linear
+-- -- Handy for missing methods.
+-- noOv :: String -> a
+-- noOv op = error (op ++ ": not defined on a :> b")
 
-instance Functor ((:>) a) where
-  fmap f (D b b') = D (f b) ((fmap.fmap) f b')
 
--- I think fmap will be meaningful only with *linear* functions.
+-- | Derivative tower full of 'zeroV'.
+dZero :: (LMapDom a s, VectorSpace b s) => a:>b
+dZero = T D.dZero
 
--- Handy for missing methods.
-noOv :: String -> a
-noOv op = error (op ++ ": not defined on a :> b")
 
-instance Applicative ((:>) a) where
-    -- pure = dConst    -- not!  see below.
-    pure = noOv "pure.  use dConst instead."
-    D f f' <*> D b b' = D (f b) (liftA2 (<*>) f' b')
+-- | Constant derivative tower.
+pureD :: (LMapDom a s, VectorSpace b s) => b -> a:>b
+pureD = fmap T D.pureD
 
--- Why can't we define 'pure' as 'dConst'?  Because of the extra type
--- constraint that @VectorSpace b@ (not @a@).  Oh well.  Be careful not to
--- use 'pure', okay?  Alternatively, I could define the '(<*>)' (naming it
--- something else) and then say @foo <$> p <*^> q <*^> ...@.
+-- | Map a /linear/ function over a derivative tower.
+fmapD, (<$>>) :: (LMapDom a s, VectorSpace b s) =>
+                 (b -> c) -> (a :> b) -> (a :> c)
+fmapD = fmap inT D.fmapD
 
--- | Constant derivative tower.
-dConst :: VectorSpace b s => b -> a:>b
-dConst b = b `D` const dZero
+(<$>>) = fmapD
 
--- | Derivative tower full of 'zeroV'.
-dZero :: VectorSpace b s => a:>b
-dZero = dConst zeroV
+-- -- | Like '(<*>)' for derivative towers.
+-- (<*>>) :: (LMapDom a s, VectorSpace b s, VectorSpace c s) =>
+--           (a :> (b -> c)) -> (a :> b) -> (a :> c)
+-- (<*>>) = inT2 (D.<*>>)
 
+-- | Apply a /linear/ binary function over derivative towers.
+liftD2 :: (VectorSpace b s, LMapDom a s, VectorSpace c s, VectorSpace d s) =>
+          (b -> c -> d) -> (a :> b) -> (a :> c) -> (a :> d)
+liftD2 = fmap inT2 D.liftD2
+
+
+-- | Apply a /linear/ ternary function over derivative towers.
+liftD3 :: ( LMapDom a s
+          , VectorSpace b s, VectorSpace c s
+          , VectorSpace d s, VectorSpace e s ) =>
+          (b -> c -> d -> e)
+       -> (a :> b) -> (a :> c) -> (a :> d) -> (a :> e)
+liftD3 = fmap inT3 D.liftD3
+
 -- | Differentiable identity function.  Sometimes called "the
 -- derivation variable" or similar, but it's not really a variable.
-idD :: VectorSpace u s => u :~> u
-idD = linearD id
-
--- or
---   dId v = D v dConst
+idD :: (LMapDom u s, VectorSpace u s) => u :~> u
+idD = fmap T D.idD
 
 -- | Every linear function has a constant derivative equal to the function
 -- itself (as a linear map).
-linearD :: VectorSpace v s => (u :-* v) -> (u :~> v)
-linearD f u = D (f u) (dConst . f)
-
+linearD :: (LMapDom u s, VectorSpace v s) => (u -> v) -> (u :~> v)
+linearD = (fmap.fmap) T D.linearD
 
 -- Other examples of linear functions
 
 -- | Differentiable version of 'fst'
-fstD :: VectorSpace a s => (a,b) :~> a
-fstD = linearD fst
+fstD :: (VectorSpace a s, LMapDom b s, LMapDom a s) => (a,b) :~> a
+fstD = fmap T D.fstD
 
 -- | Differentiable version of 'snd'
-sndD :: VectorSpace b s => (a,b) :~> b
-sndD = linearD snd
+sndD :: (VectorSpace b s, LMapDom b s, LMapDom a s) => (a,b) :~> b
+sndD = fmap T D.sndD
 
 -- | Derivative tower for applying a binary function that distributes over
--- addition, such as multiplication.
-distribD :: (VectorSpace u s) =>
-             (b -> c -> u) -> ((a :> b) -> (a :> c) -> (a :> u))
-          -> (a :> b) -> (a :> c) -> (a :> u)
-distribD op opD u@(D u0 u') v@(D v0 v') =
-  D (u0 `op` v0) ((u `opD`) . v' ^+^ (`opD` v) . u')
+-- addition, such as multiplication.  A bit weaker assumption than
+-- bilinearity.
+distrib :: (LMapDom a s, VectorSpace b s, VectorSpace c s, VectorSpace u s) =>
+           (b -> c -> u) -> (a :> b) -> (a :> c) -> (a :> u)
+distrib = fmap inT2 D.distrib
 
--- Equivalently,
--- 
---   distribD op opD u@(D u0 u') v@(D v0 v') =
---     D (u0 `op` v0) (\ da -> (u `opD` v' da) ^+^ (u' da `opD` v))
+-- I'm not sure about the next three, which discard information
 
+instance Show b => Show (a :> b) where show    = show     .   unT
+instance Eq   b => Eq   (a :> b) where (==)    = (==)    `on` unT
+instance Ord  b => Ord  (a :> b) where compare = compare `on` unT
 
--- I'm not sure about the next three, which discard information
-instance Show b => Show (a :> b) where show    = noOv "show"
-instance Eq   b => Eq   (a :> b) where (==)    = noOv "(==)"
-instance Ord  b => Ord  (a :> b) where compare = noOv "compare"
+-- These next two instances are the reason for having this module.  The
+-- scalar field differs from the one used in the underlying
+-- representation.  I can't have both instances.  First, there is a
+-- functional dependency on 'VectorSpace', which says that a vector space
+-- type determines its scalar field type.  I experimented with dropping
+-- the fundep, and then managing the resulting ambiguity.  However,
+-- the two 'VectorSpace' instances overlap considerably.
 
-instance VectorSpace u s => VectorSpace (a :> u) (a :> s) where
-  zeroV   = dConst   zeroV    -- or dZero
-  (*^)    = distribD (*^) (*^)
-  negateV = fmap     negateV
-  (^+^)   = liftA2   (^+^)
+instance (LMapDom a s, VectorSpace u s, VectorSpace s s)
+    => VectorSpace (a :> u) (a :> s) where
+  (*^)    = inT2 (D.**^) -- not '(D.*^)'
 
-instance (InnerSpace u s, InnerSpace s s') =>
+instance (InnerSpace u s, InnerSpace s s', VectorSpace s s, LMapDom a s) =>
      InnerSpace (a :> u) (a :> s) where
-  (<.>) = distribD (<.>) (<.>)
+  (<.>)    = inT2 (D.<*.>) -- not '(D.<.>)'
 
--- | Chain rule.
-(@.) :: (b :~> c) -> (a :~> b) -> (a :~> c)
-(h @. g) a0 = D c0 (c' @. b')
-  where
-    D b0 b' = g a0
-    D c0 c' = h b0
 
+-- | Chain rule.
+(@.) :: (LMapDom b s, LMapDom a s, VectorSpace c s) =>
+        (b :~> c) -> (a :~> b) -> (a :~> c)
+h @. g = T . ((unT . h) D.@. (unT . g))
 
 -- | Specialized chain rule.
-(>-<) :: VectorSpace u s => (u -> u) -> ((a :> u) -> (a :> s))
+(>-<) :: (LMapDom a s, VectorSpace s s, VectorSpace u s) =>
+         (u -> u) -> ((a :> u) -> (a :> s))
       -> (a :> u) -> (a :> u)
-
-f >-< f' = \ u@(D u0 u') -> D (f u0) ((f' u *^) . u')
-
--- Equivalently:
--- 
---   f >-< f' = \ u@(D u0 u') -> D (f u0) (\ da -> f' u *^ u' da)
-
-instance (Num b, VectorSpace b b) => Num (a:>b) where
-  fromInteger = dConst . fromInteger
-  (+) = liftA2   (+)
-  (-) = liftA2   (-)
-  (*) = distribD (*) (*)
-  negate = negate >-< -1
-  abs    = abs    >-< signum
-  signum = signum >-< 0  -- derivative wrong at zero
-
-instance (Fractional b, VectorSpace b b) => Fractional (a:>b) where
-  fromRational = dConst . fromRational
-  recip        = recip >-< recip sqr
-
-sqr :: Num a => a -> a
-sqr x = x*x
-
-instance (Floating b, VectorSpace b b) => Floating (a:>b) where
-  pi    = dConst pi
-  exp   = exp   >-< exp
-  log   = log   >-< recip
-  sqrt  = sqrt  >-< recip (2 * sqrt)
-  sin   = sin   >-< cos
-  cos   = cos   >-< - sin
-  sinh  = sinh  >-< cosh
-  cosh  = cosh  >-< sinh
-  asin  = asin  >-< recip (sqrt (1-sqr))
-  acos  = acos  >-< recip (- sqrt (1-sqr))
-  atan  = atan  >-< recip (1+sqr)
-  asinh = asinh >-< recip (sqrt (1+sqr))
-  acosh = acosh >-< recip (- sqrt (sqr-1))
-  atanh = atanh >-< recip (1-sqr)
+f >-< f' = inT (f D.>-< (unT . f' . T))
diff --git a/src/Data/Horner.hs b/src/Data/Horner.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Horner.hs
@@ -0,0 +1,220 @@
+{-# LANGUAGE TypeOperators, MultiParamTypeClasses, UndecidableInstances
+           , TypeSynonymInstances, FlexibleInstances, FunctionalDependencies
+  #-}
+{-# OPTIONS_GHC -Wall #-}
+----------------------------------------------------------------------
+-- |
+-- Module      :  Data.Horner
+-- Copyright   :  (c) Conal Elliott 2008
+-- License     :  BSD3
+-- 
+-- Maintainer  :  conal@conal.net
+-- Stability   :  experimental
+-- 
+-- Infinite derivative towers via linear maps, using the Horner
+-- representation.  See blog posts <http://conal.net/blog/tag/derivatives/>.
+----------------------------------------------------------------------
+
+module Data.Horner
+  (
+    (:>), powVal, derivative, integral
+  , (:~>), dZero, dConst
+  , idD, fstD, sndD
+  , linearD, distrib
+  , (@.), (>-<)
+  -- , HasDeriv(..)
+  )
+  where
+
+import Control.Applicative
+
+import Data.VectorSpace
+import Data.LinearMap
+import Data.NumInstances ()
+
+infixr 9 `H`, @.
+infix  0 >-<
+
+-- | Power series
+-- 
+-- Warning, the 'Applicative' instance is missing its 'pure' (due to a
+-- 'VectorSpace' type constraint).  Use 'dConst' instead.
+data a :> b = H b (a :-* (a :> b))
+
+-- | The plain-old (0th order) value
+powVal :: (a :> b) -> b
+powVal (H b _) = b
+
+-- Apply successive functions to successive values
+apPow :: [b -> c] -> (a :> b) -> (a :> c)
+apPow [] _ = error "apPow: finite function list"
+apPow (f : fs) (b0 `H` bt) = H (f b0) (apPow fs . bt)
+
+-- Count.  Avoids the 'Enum' requirement of [1..]
+from :: Num s => s -> [s]
+from n = n : from (n+1) 
+
+-- | Derivative of a power series
+derivative :: (VectorSpace b s, Num s) =>
+         (a :> b) -> (a :-* (a :> b))
+derivative (H _ bt) = apPow ((*^) <$> from 1) . bt
+
+-- | Integral of a power series
+integral :: (VectorSpace b s, Fractional s) =>
+            b -> (a :-* (a :> b)) -> (a :> b)
+integral b0 bt = H b0 (apPow (((*^).recip) <$> from 1) . bt)
+
+-- | Infinitely differentiable functions
+type a :~> b = a -> (a:>b)
+
+-- So we could define
+-- 
+--   data a :> b = H b (a :~> b)
+-- 
+-- with the restriction that the a :~> b is linear
+
+instance Functor ((:>) a) where
+  fmap f (H b b') = H (f b) ((fmap.fmap) f b')
+
+-- I think fmap will be meaningful only with *linear* functions.
+
+-- Handy for missing methods.
+noOv :: String -> a
+noOv op = error (op ++ ": not defined on a :> b")
+
+instance Applicative ((:>) a) where
+  -- pure = dConst    -- not!  see below.
+  pure = noOv "pure"  -- use dConst instead
+  H f f' <*> H b b' = H (f b) (liftA2 (<*>) f' b')
+
+-- Why can't we define 'pure' as 'dConst'?  Because of the extra type
+-- constraint that @VectorSpace b@ (not @a@).  Oh well.  Be careful not to
+-- use 'pure', okay?  Alternatively, I could define the '(<*>)' (naming it
+-- something else) and then say @foo <$> p <*^> q <*^> ...@.
+
+-- | Constant derivative tower.
+dConst :: VectorSpace b s => b -> a:>b
+dConst b = b `H` const dZero
+
+-- | Derivative tower full of 'zeroV'.
+dZero :: VectorSpace b s => a:>b
+dZero = dConst zeroV
+
+-- | Differentiable identity function.  Sometimes called "the
+-- derivation variable" or similar, but it's not really a variable.
+idD :: VectorSpace u s => u :~> u
+idD = linearD id
+
+-- or
+--   dId v = H v dConst
+
+-- | Every linear function has a constant derivative equal to the function
+-- itself (as a linear map).
+linearD :: VectorSpace v s => (u :-* v) -> (u :~> v)
+linearD f u = H (f u) (dConst . f)
+
+
+-- Other examples of linear functions
+
+-- | Differentiable version of 'fst'
+fstD :: VectorSpace a s => (a,b) :~> a
+fstD = linearD fst
+
+-- | Differentiable version of 'snd'
+sndD :: VectorSpace b s => (a,b) :~> b
+sndD = linearD snd
+
+-- | Derivative tower for applying a binary function that distributes over
+-- addition, such as multiplication.  A bit weaker assumption than
+-- bilinearity.
+distrib :: (VectorSpace u s) =>
+           (b -> c -> u) -> (a :> b) -> (a :> c) -> (a :> u)
+distrib op = opD
+ where
+   opD (H u0 ut) v@(H v0 vt) =
+     H (u0 `op` v0) (fmap (u0 `op`) . vt ^+^ (`opD` v) . ut)
+
+
+-- Equivalently,
+-- 
+--   distrib op = opD
+--    where
+--      opD u@(H u0 u') v@(H v0 v') =
+--        H (u0 `op` v0) (\ da -> ((u0 `op`) <$> v' da) ^+^ (u' da `opD` v))
+
+
+
+-- I'm not sure about the next three, which discard information
+
+instance Show b => Show (a :> b) where show    = noOv "show"
+instance Eq   b => Eq   (a :> b) where (==)    = noOv "(==)"
+instance Ord  b => Ord  (a :> b) where compare = noOv "compare"
+
+instance (LMapDom a s, VectorSpace u s) => AdditiveGroup (a :> u) where
+  zeroV   = pureD  zeroV    -- or dZero
+  negateV = fmapD  negateV
+  (^+^)   = liftD2 (^+^)
+
+instance (LMapDom a s, VectorSpace u s) => VectorSpace (a :> u) s where
+  (*^) s = fmapD  ((*^) s)
+
+(**^) :: (VectorSpace c s, VectorSpace s s, LMapDom a s) =>
+         (a :> s) -> (a :> c) -> (a :> c)
+(**^) = distrib (*^)
+
+-- | Chain rule.
+(@.) :: (VectorSpace b s, VectorSpace c s, Num s) =>
+        (b :~> c) -> (a :~> b) -> (a :~> c)
+(h @. g) a0 = H c0 (derivative c @. derivative b)
+  where
+    b@(H b0 _) = g a0
+    c@(H c0 _) = h b0
+
+
+-- | Specialized chain rule.
+(>-<) :: (VectorSpace u s, Fractional s) => (u -> u) -> ((a :> u) -> (a :> s))
+      -> (a :> u) -> (a :> u)
+
+-- f >-< f' = \ u@(D u0 u') -> D (f u0) ((f' u *^) . u')
+
+f >-< f' = \ u@(H u0 _) -> integral (f u0) ((f' u *^) . derivative u)
+
+-- TODO: consider eliminating @Num s@.  I just need a multiplicative unit.
+
+-- Equivalently:
+-- 
+--   f >-< f' = \ u@(H u0 u') -> H (f u0) (\ da -> f' u *^ u' da)
+
+instance (Fractional b, VectorSpace b b) => Num (a:>b) where
+  fromInteger = dConst . fromInteger
+  (+) = liftA2  (+)
+  (-) = liftA2  (-)
+  (*) = distrib (*)
+  
+  negate = negate >-< -1
+  abs    = abs    >-< signum
+  signum = signum >-< 0  -- derivative wrong at zero
+
+instance (Fractional b, VectorSpace b b) => Fractional (a:>b) where
+  fromRational = dConst . fromRational
+  recip        = recip >-< recip sqr
+
+sqr :: Num a => a -> a
+sqr x = x*x
+
+instance (Floating b, VectorSpace b b) => Floating (a:>b) where
+  pi    = dConst pi
+  exp   = exp   >-< exp
+  log   = log   >-< recip
+  sqrt  = sqrt  >-< recip (2 * sqrt)
+  sin   = sin   >-< cos
+  cos   = cos   >-< - sin
+  sinh  = sinh  >-< cosh
+  cosh  = cosh  >-< sinh
+  asin  = asin  >-< recip (sqrt (1-sqr))
+  acos  = acos  >-< recip (- sqrt (1-sqr))
+  atan  = atan  >-< recip (1+sqr)
+  asinh = asinh >-< recip (sqrt (1+sqr))
+  acosh = acosh >-< recip (- sqrt (sqr-1))
+  atanh = atanh >-< recip (1-sqr)
+
diff --git a/src/Data/LinearMap.hs b/src/Data/LinearMap.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/LinearMap.hs
@@ -0,0 +1,315 @@
+{-# LANGUAGE TypeOperators, TypeFamilies, UndecidableInstances
+           , MultiParamTypeClasses, FlexibleInstances
+           , FunctionalDependencies
+  #-}
+{-# OPTIONS_GHC -Wall -fno-warn-orphans -fglasgow-exts #-}
+----------------------------------------------------------------------
+-- |
+-- Module      :  Data.LinearMap
+-- Copyright   :  (c) Conal Elliott 2008
+-- License     :  BSD3
+-- 
+-- Maintainer  :  conal@conal.net
+-- Stability   :  experimental
+-- 
+-- Linear maps
+----------------------------------------------------------------------
+
+module Data.LinearMap
+  (
+    LMapDom(..), inL, inL2, inL3
+  , linearK, (.*)
+  , pureL --, (<*>*)
+  , fmapL, (<$>*), liftL2, liftL3, idL
+  ) where
+
+-- -fglasgow-exts above enables the RULES pragma
+
+import Control.Applicative
+import Data.Function
+
+import Data.VectorSpace
+
+-- Temporary
+import Graphics.Rendering.OpenGL.GL.CoordTrans
+  (Vector2(..),Vector3(..))
+
+
+infixr 9 :-*
+infixr 9 .*
+infixl 9 `lapply`
+infixl 4 <$>* -- , <*>*
+
+
+-- | Domain of a linear map.
+class VectorSpace a s => LMapDom a s | a -> s where
+  -- | Linear map type
+  data (:-*) a :: * -> *
+  -- | Linear map as function
+  lapply :: VectorSpace b s => (a :-* b) -> (a -> b)
+  -- | Function (assumed linear) as linear map.
+  linear :: (a -> b) -> (a :-* b)
+
+-- Neither 'VectorSpace' nor even 'AdditiveGroup' is really required as a
+-- 'LMapDom' superclass.  Instead, we could have additional constraints in
+-- some 'LMapDom' instances and related functions.
+
+
+{-# RULES
+"linear.lapply"   forall m. linear (lapply m) = m
+"lapply.linear"   forall f. lapply (linear f) = f
+ #-}
+
+
+-- | Transform a linear map by transforming a linear function.
+inL :: (LMapDom c s, VectorSpace b s', LMapDom a s') =>
+        ((a -> b) -> (c -> d)) -> ((a :-* b) -> (c :-* d))
+{-# INLINE inL #-}
+inL h = linear . h . lapply
+
+-- | Transform a linear maps by transforming linear functions.
+inL2 :: ( LMapDom c s, VectorSpace b s', LMapDom a s'
+        , LMapDom e s, VectorSpace d s ) =>
+        ((a -> b) -> (c -> d) -> (e -> f))
+     -> ((a :-* b) -> (c :-* d) -> (e :-* f))
+{-# INLINE inL2 #-}
+inL2 h = inL . h . lapply
+
+-- inL2 h m n
+--   = (inL . h . lapply) m n
+--   = inL (h (lapply m)) n
+--   = (linear . (h (lapply m)) . lapply) n
+--   = linear (h (lapply m) (lapply n)
+--   = linear (h (lapply m) (lapply n))
+
+
+-- | Transform a linear maps by transforming linear functions.
+inL3 :: ( LMapDom a s, VectorSpace b s
+        , VectorSpace f s , LMapDom p s, LMapDom c s'
+        , VectorSpace d s', LMapDom e s ) =>
+        ((a -> b) -> (c -> d) -> (e -> f) -> (p -> q))
+     -> ((a :-* b) -> (c :-* d) -> (e :-* f) -> (p :-* q))
+{-# INLINE inL3 #-}
+inL3 h = inL2 . h . lapply
+
+
+-- TODO: go through this whole module and relax constraints on scalar
+-- fields.  See if it helps with the scalar dilemma in Mac2
+
+
+-- | Constant value as a linear map
+pureL :: LMapDom a s => b -> (a :-* b)
+pureL b = linear (const b)
+
+-- -- | Like '(<*>)' for linear maps.
+-- (<*>*) :: (LMapDom a s, VectorSpace b s, VectorSpace c s) =>
+--         (a :-* (b -> c)) -> (a :-* b) -> (a :-* c)
+-- (<*>*) = inL2 (<*>)
+
+-- | Map a /linear/ function over a linear map.
+fmapL, (<$>*) :: (LMapDom a s, VectorSpace b s) =>
+                 (b -> c) -> (a :-* b) -> (a :-* c)
+{-# INLINE fmapL #-}
+fmapL = inL . fmap
+
+-- fmapL f
+--   = inL (f .)
+--   = linear . (f .) . lapply
+--   = \ m -> linear (fmap f (lapply m))
+
+(<$>*) = fmapL
+
+
+{-# RULES
+"fmapL.fmapL"  forall f g m. fmapL f (fmapL g m) = fmapL (f.g) m
+ #-}
+
+-- | Apply a /linear/ binary function over linear maps.
+liftL2 :: ( LMapDom a s, VectorSpace b s
+           , VectorSpace c s, VectorSpace d s) =>
+           (b -> c -> d) -> (a :-* b) -> (a :-* c) -> (a :-* d)
+liftL2 = inL2 . liftA2
+
+-- liftL2 f a b
+--   = inL2 (liftA2 f) a b
+--   = linear (liftA2 f (lapply a) (lapply b))
+
+-- I expected the following definition to be equivalent, thanks to rewriting:
+-- 
+--   liftL2 f b c = fmapL f b <*>* c
+--     = linear (f . lapply b) <*>* c
+--     = linear (lapply (linear (f . lapply b)) <*> lapply c)
+--     = linear ((f . lapply b) <*> lapply c)
+--     = linear (liftA2 f (lapply b) (lapply c))
+--     = (inL2.liftA2) f b c
+
+-- The rewrite isn't happening, however.  And the '(<*>*)' definition
+-- yields an incredibly slow implementation.
+-- 
+-- Here's that derivation again, in slo-mo:
+
+--   liftL2 f b c
+--     = fmapL f b <*>* c                                  -- inline liftL2
+--     = inL (f .) b <*>* c                                -- inline fmapL
+--     = (linear . (f .) lapply) b <*>* c                  -- inline inL
+--     = linear (f . lapply b) <*>* c                      -- inline (.)
+--     = inL2 (<*>) (linear (f . lapply b)) c              -- inline (<*>*)
+--     = (inL . (<*>) . lapply) (linear (f . lapply b)) c  -- inline inL2
+--     = inL ((<*>) (lapply (linear (f . lapply b)))) c    -- inline (.)
+--     = inL ((<*>) (f . lapply b)) c                      -- RULE lapply.linear
+--     = (linear . ((<*>) (f . lapply b)) .lapply) c       -- inline inL
+--     = linear ((<*>) (f . lapply b) (lapply c))          -- inline (.)
+--     = linear ((f . lapply b) <*> (lapply c))            -- infix <*>
+--     = linear ((f <$> lapply b) <*> (lapply c))          -- <$> on (a ->)
+--     = linear (liftA2 f (lapply b) (lapply c))           -- uninline liftA2
+
+-- When I compile this module, I don't get any firings of lapply.linear
+
+
+-- | Apply a /linear/ ternary function over linear maps.
+liftL3 :: ( LMapDom a s, VectorSpace b s, VectorSpace c s
+           , VectorSpace d s, VectorSpace e s) =>
+           (b -> c -> d -> e)
+        -> (a :-* b) -> (a :-* c) -> (a :-* d) -> (a :-* e)
+liftL3 = inL3 . liftA3
+
+
+-- TODO: Get clear about the linearity requirements of 'apL', 'liftL2',
+-- and 'liftL3'.
+
+-- | Identity linear map
+idL :: LMapDom a s => a :-* a
+idL = linear id
+
+-- | Compose linear maps
+(.*) :: (VectorSpace c s, LMapDom b s, LMapDom a s) =>
+        (b :-* c) -> (a :-* b) -> (a :-* c)
+(.*) = inL2 (.)
+
+
+instance (VectorSpace v s, LMapDom u s) => AdditiveGroup (u :-* v) where
+  zeroV   = pureL  zeroV
+  (^+^)   = liftL2 (^+^)
+  negateV = fmapL  negateV
+
+instance (VectorSpace v s, LMapDom u s) => VectorSpace (u :-* v) s where
+  (*^) s  = fmapL  ((*^) s)
+
+-- Or possibly the following non-standard definition:
+
+-- instance (VectorSpace v s, VectorSpace s s, LMapDom u s)
+--       => VectorSpace (u :-* v) (u :-* s) where
+--   zeroV   = pureL  zeroV
+--   (*^)    = liftL2 (*^)
+--   (^+^)   = liftL2 (^+^)
+--   negateV = fmapL  negateV
+
+-- Alternatively, add some methods to 'LMapDom' and use them as follows.
+-- May be more efficient.
+
+-- -- Linear maps form a vector space
+-- instance (VectorSpace o s, LinearMap a o s) => VectorSpace (a :--> o) s where
+--   zeroV   = zeroL
+--   (^+^)   = addL
+--   (*^)    = scaleL
+--   negateV = negateL
+
+
+
+--- Instances of LMapDom
+
+instance LMapDom Float Float where
+  data Float :-* o   = FloatL o
+  lapply (FloatL o)  = (*^ o)
+  linear f           = FloatL (f 1)
+
+instance LMapDom Double Double where
+  data Double :-* o  = DoubleL o
+  lapply (DoubleL o) = (*^ o)
+  linear f           = DoubleL (f 1)
+
+
+-- | Convenience function for 'linear' definitions.  Both functions are
+-- assumed linear.
+linearK :: (LMapDom a s) => (a -> b) -> (b -> c) -> a :-* c
+linearK k f = linear (f . k)
+
+-- instance LMapDom (Double,Double) Double where
+--   data (Double,Double) :-* o = PairD o o
+--   PairD ao bo `lapply` (a,b) = a *^ ao ^+^ b *^ bo
+--   linear f = PairD (f (0,1)) (f (1,0))
+
+instance (LMapDom a s, LMapDom b s) => LMapDom (a,b) s where
+  data (a,b) :-* o = PairL (a :-* o) (b :-* o)
+  PairL ao bo `lapply` (a,b) = ao `lapply` a ^+^ bo `lapply` b
+  linear f = PairL (linear (\ a -> f (a,zeroV)))
+                   (linear (\ b -> f (zeroV,b)))
+
+--   linear = liftA2 PairL
+--              (linearK (\ a -> (a,zeroV)))
+--              (linearK (\ b -> (zeroV,b)))
+
+instance (LMapDom a s, LMapDom b s, LMapDom c s) => LMapDom (a,b,c) s where
+  data (a,b,c) :-* o = TripleL (a :-* o) (b :-* o) (c :-* o)
+  TripleL ao bo co `lapply` (a,b,c) =
+    ao `lapply` a ^+^ bo `lapply` b ^+^ co `lapply` c
+  linear = liftA3 TripleL
+             (linearK (\ a -> (a,zeroV,zeroV)))
+             (linearK (\ b -> (zeroV,b,zeroV)))
+             (linearK (\ c -> (zeroV,zeroV,c)))
+
+
+
+
+-- TODO: unfst, unsnd, pair, unpair
+
+
+---- OpenGL stuff.
+
+-- I'd rather this code be in a different package.  It's here as a
+-- temporary bug workaround.  In ghc-6.8.2, I get the following error
+-- message if the 'LMapDom' instance (below) is compiled in a separate
+-- module.
+-- 
+--     Type indexes must match class instance head
+--     Found o but expected Vector2 u
+--     In the associated type instance for `:-*'
+--     In the instance declaration for `LMapDom (Vector2 u) s'
+
+
+
+-- TODO: is UndecidableInstances still necessary?
+
+instance AdditiveGroup u => AdditiveGroup (Vector2 u) where
+  zeroV                         = Vector2 zeroV zeroV
+  Vector2 u v ^+^ Vector2 u' v' = Vector2 (u^+^u') (v^+^v')
+  negateV (Vector2 u v)         = Vector2 (negateV u) (negateV v)
+
+instance (VectorSpace u s) => VectorSpace (Vector2 u) s where
+  s *^ Vector2 u v            = Vector2 (s*^u) (s*^v)
+
+instance (InnerSpace u s, AdditiveGroup s)
+    => InnerSpace (Vector2 u) s where
+  Vector2 u v <.> Vector2 u' v' = u<.>u' ^+^ v<.>v'
+
+instance LMapDom u s => LMapDom (Vector2 u) s where
+  data Vector2 u :-* o = VecL (u :-* o) (u :-* o)
+  VecL ao bo `lapply` Vector2 a b = ao `lapply` a ^+^ bo `lapply` b
+  linear = liftA2 VecL
+             (linearK (\ a -> Vector2 a zeroV))
+             (linearK (\ b -> Vector2 zeroV b))
+
+
+instance AdditiveGroup u => AdditiveGroup (Vector3 u) where
+  zeroV                   = Vector3 zeroV zeroV zeroV
+  Vector3 u v w ^+^ Vector3 u' v' w'
+                          = Vector3 (u^+^u') (v^+^v') (w^+^w')
+  negateV (Vector3 u v w) = Vector3 (negateV u) (negateV v) (negateV w)
+
+instance VectorSpace u s => VectorSpace (Vector3 u) s where
+  s *^ Vector3 u v w    = Vector3 (s*^u) (s*^v) (s*^w)
+
+instance (InnerSpace u s, AdditiveGroup s)
+    => InnerSpace (Vector3 u) s where
+  Vector3 u v w <.> Vector3 u' v' w' = u<.>u' ^+^ v<.>v' ^+^ w<.>w'
+
diff --git a/src/Data/Maclaurin.hs b/src/Data/Maclaurin.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Maclaurin.hs
@@ -0,0 +1,248 @@
+{-# LANGUAGE TypeOperators, MultiParamTypeClasses, UndecidableInstances
+           , TypeSynonymInstances, FlexibleInstances, FunctionalDependencies
+           , FlexibleContexts
+  #-}
+
+-- TODO: remove FlexibleContexts
+
+{-# OPTIONS_GHC -Wall #-}
+----------------------------------------------------------------------
+-- |
+-- Module      :  Data.Maclaurin
+-- Copyright   :  (c) Conal Elliott 2008
+-- License     :  BSD3
+-- 
+-- Maintainer  :  conal@conal.net
+-- Stability   :  experimental
+-- 
+-- Infinite derivative towers via linear maps, using the Maclaurin
+-- representation.  See blog posts <http://conal.net/blog/tag/derivatives/>.
+----------------------------------------------------------------------
+
+module Data.Maclaurin
+  (
+    (:>), powVal, derivative, derivativeAt
+  , (:~>), dZero, pureD
+  , fmapD, (<$>>){-, (<*>>)-}, liftD2, liftD3
+  , idD, fstD, sndD
+  , linearD, distrib
+  , (@.), (>-<)
+  ,(**^), (<*.>)
+  -- , HasDeriv(..)
+  -- experimental
+  -- , liftD3
+  ) where
+
+-- import Control.Applicative
+
+import Data.VectorSpace
+import Data.NumInstances ()
+import Data.LinearMap
+
+
+infixr 9 `D`, @.
+infixl 4 {-<*>>,-} <$>>
+infix  0 >-<
+
+
+-- | Tower of derivatives.
+data a :> b = D { powVal :: b, derivative :: a :-* (a :> b) }
+
+-- | Infinitely differentiable functions
+type a :~> b = a -> (a:>b)
+
+-- | Sampled derivative.  For avoiding an awkward typing problem related
+-- to the two required 'VectorSpace' instances.
+derivativeAt :: (VectorSpace b s, LMapDom a s) =>
+                (a :> b) -> a -> (a :> b)
+derivativeAt d = lapply (derivative d)
+
+-- The crucial point here is for '($*)' to be interpreted with respect to
+-- the 'VectorSpace' instance in this module, not Mac.
+
+-- The argument order for 'derivativeAt' allows partial evaluation, which
+-- is useful in power series representations for which 'derivative' is not
+-- free (Horner).
+
+-- Handy for missing methods.
+noOv :: String -> a
+noOv op = error (op ++ ": not defined on a :> b")
+
+-- | Derivative tower full of 'zeroV'.
+dZero :: (LMapDom a s, AdditiveGroup b) => a:>b
+dZero = pureD zeroV
+
+-- | Constant derivative tower.
+pureD :: (LMapDom a s, AdditiveGroup b) => b -> a:>b
+pureD b = b `D` pureL dZero
+
+-- | Map a /linear/ function over a derivative tower.
+fmapD, (<$>>) :: (LMapDom a s, VectorSpace b s) =>
+                 (b -> c) -> (a :> b) -> (a :> c)
+fmapD f (D b0 b') = D (f b0) ((fmapL.fmapD) f b')
+
+(<$>>) = fmapD
+
+-- -- | Like '(<*>)' for derivative towers.
+-- (<*>>) :: (LMapDom a s, VectorSpace b s, VectorSpace c s) =>
+--           (a :> (b -> c)) -> (a :> b) -> (a :> c)
+-- D f0 f' <*>> D x0 x' = D (f0 x0) (liftL2 (<*>>) f' x')
+
+-- | Apply a /linear/ binary function over derivative towers.
+liftD2 :: (VectorSpace b s, LMapDom a s, VectorSpace c s, VectorSpace d s) =>
+          (b -> c -> d) -> (a :> b) -> (a :> c) -> (a :> d)
+liftD2 f (D b0 b') (D c0 c') = D (f b0 c0) (liftL2 (liftD2 f) b' c')
+
+-- | Apply a /linear/ ternary function over derivative towers.
+liftD3 :: ( LMapDom a s
+          , VectorSpace b s, VectorSpace c s
+          , VectorSpace d s, VectorSpace e s ) =>
+          (b -> c -> d -> e)
+       -> (a :> b) -> (a :> c) -> (a :> d) -> (a :> e)
+liftD3 f (D b0 b') (D c0 c') (D d0 d') = D (f b0 c0 d0) (liftL3 (liftD3 f) b' c' d')
+
+-- | Differentiable identity function.  Sometimes called "the
+-- derivation variable" or similar, but it's not really a variable.
+idD :: (LMapDom u s, VectorSpace u s) => u :~> u
+idD = linearD id
+
+-- or
+--   dId v = D v pureD
+
+-- | Every linear function has a constant derivative equal to the function
+-- itself (as a linear map).
+linearD :: (LMapDom u s, VectorSpace v s) => (u -> v) -> (u :~> v)
+linearD f u = D (f u) (linear (pureD . f))
+
+-- Other examples of linear functions
+
+-- | Differentiable version of 'fst'
+fstD :: (VectorSpace a s, LMapDom b s, LMapDom a s) => (a,b) :~> a
+fstD = linearD fst
+
+-- | Differentiable version of 'snd'
+sndD :: (VectorSpace b s, LMapDom b s, LMapDom a s) => (a,b) :~> b
+sndD = linearD snd
+
+-- | Derivative tower for applying a binary function that distributes over
+-- addition, such as multiplication.  A bit weaker assumption than
+-- bilinearity.
+distrib :: (LMapDom a s, VectorSpace b s, VectorSpace c s, VectorSpace u s) =>
+           (b -> c -> u) -> (a :> b) -> (a :> c) -> (a :> u)
+distrib op = opD
+ where
+   opD u@(D u0 u') v@(D v0 v') =
+     D (u0 `op` v0) (linear (\ da -> u `opD` (v' `lapply` da) ^+^
+                                     (u' `lapply` da) `opD` v))
+
+-- Equivalently:
+-- 
+--    opD u@(D u0 u') v@(D v0 v') =
+--      D (u0 `op` v0) (linear ((u `opD`) . lapply v' ^+^ (`opD` v) . lapply u'))
+-- 
+-- or
+-- 
+--    opD u@(D u0 u') v@(D v0 v') =
+--      D (u0 `op` v0) ( linear ((u `opD`) . lapply v') ^+^
+--                       linear ((`opD` v) . lapply u') )
+-- or even
+-- 
+--    opD u@(D u0 u') v@(D v0 v') =
+--      D (u0 `op` v0) ( inL ((u `opD`) .) v' ^+^ inL ((`opD` v) .) u' )
+
+
+
+
+
+-- TODO: look for a simpler definition of distrib.  this definition almost
+-- fits liftLM2.
+
+
+-- I'm not sure about the next three, which discard information
+
+instance Show b => Show (a :> b) where show    = noOv "show"
+instance Eq   b => Eq   (a :> b) where (==)    = noOv "(==)"
+instance Ord  b => Ord  (a :> b) where compare = noOv "compare"
+
+instance (LMapDom a s, VectorSpace u s) => AdditiveGroup (a :> u) where
+  zeroV   = pureD  zeroV    -- or dZero
+  negateV = fmapD  negateV
+  (^+^)   = liftD2 (^+^)
+
+instance (LMapDom a s, VectorSpace u s) => VectorSpace (a :> u) s where
+  (*^) s = fmapD  ((*^) s)
+
+(**^) :: (VectorSpace c s, VectorSpace s s, LMapDom a s) =>
+         (a :> s) -> (a :> c) -> (a :> c)
+(**^) = distrib (*^)
+
+-- ouch!  InnerSpace one won't work at all, for the same reason as for functions.
+                  
+-- instance (InnerSpace u s) => InnerSpace (a :> u) s where
+--   (<.>) = distrib (<.>)
+
+(<*.>) :: (LMapDom a s, InnerSpace b s, VectorSpace s s) =>
+          (a :> b) -> (a :> b) -> (a :> s)
+(<*.>) s = distrib (<.>) s
+
+
+-- The instances below are the one I think we'll want externally.
+-- However, the ones above allow the definition of @a:>b@ to work out.
+-- The module "Data.Mac" rewraps to provide the alternate instances.
+
+-- instance (LMapDom a s, VectorSpace u s, VectorSpace s s)
+--     => VectorSpace (a :> u) (a :> s) where
+--   (*^) = (**^)
+
+-- instance (InnerSpace u s, InnerSpace s s', VectorSpace s s, LMapDom a s) =>
+--      InnerSpace (a :> u) (a :> s) where
+--   (<.>) = (<*.>)
+
+
+-- | Chain rule.  See also '(>-<)'.
+(@.) :: (LMapDom b s, LMapDom a s, VectorSpace c s) =>
+        (b :~> c) -> (a :~> b) -> (a :~> c)
+(h @. g) a0 = D c0 (inL2 (@.) c' b')
+  where
+    D b0 b' = g a0
+    D c0 c' = h b0
+
+-- | Specialized chain rule.  See also '(\@.)'
+(>-<) :: (LMapDom a s, VectorSpace s s, VectorSpace u s) =>
+         (u -> u) -> ((a :> u) -> (a :> s))
+      -> (a :> u) -> (a :> u)
+f >-< f' = \ u@(D u0 u') -> D (f u0) ((f' u **^) <$>* u')
+
+-- TODO: express '(>-<)' in terms of '(@.)'.  If I can't, then understand why not.
+
+instance (LMapDom a b, Num b, VectorSpace b b) => Num (a:>b) where
+  fromInteger = pureD . fromInteger
+  (+) = liftD2  (+)
+  (-) = liftD2  (-)
+  (*) = distrib (*)
+  negate = negate >-< -1
+  abs    = abs    >-< signum
+  signum = signum >-< 0  -- derivative wrong at zero
+
+instance (LMapDom a b, Fractional b, VectorSpace b b) => Fractional (a:>b) where
+  fromRational = pureD . fromRational
+  recip        = recip >-< recip sqr
+
+sqr :: Num a => a -> a
+sqr x = x*x
+
+instance (LMapDom a b, Floating b, VectorSpace b b) => Floating (a:>b) where
+  pi    = pureD pi
+  exp   = exp   >-< exp
+  log   = log   >-< recip
+  sqrt  = sqrt  >-< recip (2 * sqrt)
+  sin   = sin   >-< cos
+  cos   = cos   >-< - sin
+  sinh  = sinh  >-< cosh
+  cosh  = cosh  >-< sinh
+  asin  = asin  >-< recip (sqrt (1-sqr))
+  acos  = acos  >-< recip (- sqrt (1-sqr))
+  atan  = atan  >-< recip (1+sqr)
+  asinh = asinh >-< recip (sqrt (1+sqr))
+  acosh = acosh >-< recip (- sqrt (sqr-1))
+  atanh = atanh >-< recip (1-sqr)
diff --git a/src/Data/VectorSpace.hs b/src/Data/VectorSpace.hs
--- a/src/Data/VectorSpace.hs
+++ b/src/Data/VectorSpace.hs
@@ -1,5 +1,5 @@
 {-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies 
-           , FlexibleInstances, FlexibleContexts, UndecidableInstances
+           , FlexibleInstances, UndecidableInstances
  #-}
 ----------------------------------------------------------------------
 -- |
@@ -15,54 +15,30 @@
 
 module Data.VectorSpace
   ( 
-    VectorSpace(..), (^-^), (^/), (^*)
-  , InnerSpace(..) --, Scalar
+    module Data.AdditiveGroup
+  , VectorSpace(..), (^/), (^*)
+  , InnerSpace(..)
   , lerp, magnitudeSq, magnitude, normalized
-  , (:-*)
   ) where
 
-import Control.Applicative
 import Data.Complex hiding (magnitude)
 
-infixr 9 :-*
+import Data.AdditiveGroup
+
 infixr 7 *^, ^/, <.>
 infixl 7 ^*
-infixl 6 ^+^, ^-^
 
--- | Vector space @v@ over a scalar field @s@
-class VectorSpace v s | v -> s where
-  -- | The zero vector
-  zeroV :: v
+-- | Vector space @v@ over a scalar field @s@.  Extends 'AdditiveGroup'
+-- with scalar multiplication.
+class AdditiveGroup v => VectorSpace v s | v -> s where
   -- | Scale a vector
   (*^)  :: s -> v -> v
-  -- | Add vectors
-  (^+^) :: v -> v -> v
-  -- | Additive inverse
-  negateV :: v -> v
 
--- | Adds inner (dot) products
-class VectorSpace v s => InnerSpace v s | v -> s where
+-- | Adds inner (dot) products.
+class VectorSpace v s => InnerSpace v s where
   -- | Inner/dot product
   (<.>) :: v -> v -> s
 
--- | Convenience.  Maybe add methods later.
--- class VectorSpace s s => Scalar s
-
--- TODO: consider replacing v with a type constructor argument:
--- 
--- class VectorSpace v where
---   zeroV :: v s
---   (*^)  :: s -> v s -> v s
---   (^+^) :: v s -> v s -> v s
---   (<.>)   :: v s -> v s -> s
--- 
--- Perhaps with constraints on s.  We couldn't then define instances for
--- doubles & floats.
-
--- | Vector subtraction
-(^-^) :: VectorSpace v s => v -> v -> v
-v ^-^ v' = v ^+^ negateV v'
-
 -- | Vector divided by scalar
 (^/) :: (Fractional s, VectorSpace v s) => v -> s -> v
 v ^/ s = (1/s) *^ v
@@ -77,7 +53,7 @@
 
 -- | Square of the length of a vector.  Sometimes useful for efficiency.
 -- See also 'magnitude'.
-magnitudeSq :: InnerSpace v s =>  v -> s
+magnitudeSq :: InnerSpace v s => v -> s
 magnitudeSq v = v <.> v
 
 -- | Length of a vector.   See also 'magnitudeSq'.
@@ -87,34 +63,16 @@
 -- | Vector in same direction as given one but with length of one.  If
 -- given the zero vector, then return it.
 normalized :: (InnerSpace v s, Floating s) =>  v -> v
-normalized v | mag /= 0  = v ^/ mag
-             | otherwise = v
-  where
-    mag = magnitude v
-
-instance VectorSpace Double Double where
-  zeroV   = 0.0
-  (*^)    = (*)
-  (^+^)   = (+)
-  negateV = negate
-
-instance InnerSpace Double Double where
- (<.>) = (*)
+normalized v = v ^/ magnitude v
 
-instance VectorSpace Float Float where
-  zeroV   = 0.0
-  (*^)    = (*)
-  (^+^)   = (+)
-  negateV = negate
+instance VectorSpace Double Double where (*^)  = (*)
+instance InnerSpace  Double Double where (<.>) = (*)
 
-instance InnerSpace Float Float where
-  (<.>) = (*)
+instance VectorSpace Float  Float  where (*^)  = (*)
+instance InnerSpace  Float  Float  where (<.>) = (*)
 
 instance (RealFloat v, VectorSpace v s) => VectorSpace (Complex v) s where
-  zeroV       = zeroV :+ zeroV
   s*^(u :+ v) = s*^u :+ s*^v
-  (^+^)       = (+)
-  negateV     = negate
 
 instance (RealFloat v, InnerSpace v s, VectorSpace s s')
      => InnerSpace (Complex v) s where
@@ -129,10 +87,7 @@
 --   Coverage Condition fails for one of the functional dependencies ...)
 
 instance (VectorSpace u s,VectorSpace v s) => VectorSpace (u,v) s where
-  zeroV             = (zeroV,zeroV)
-  s *^ (u,v)        = (s*^u,s*^v)
-  (u,v) ^+^ (u',v') = (u^+^u',v^+^v')
-  negateV (u,v)     = (negateV u, negateV v)
+  s *^ (u,v) = (s*^u,s*^v)
 
 instance (InnerSpace u s,InnerSpace v s, VectorSpace s s')
     => InnerSpace (u,v) s where
@@ -143,48 +98,15 @@
 
 instance (VectorSpace u s,VectorSpace v s,VectorSpace w s)
     => VectorSpace (u,v,w) s where
-  zeroV                  = (zeroV,zeroV,zeroV)
-  s *^ (u,v,w)           = (s*^u,s*^v,s*^w)
-  (u,v,w) ^+^ (u',v',w') = (u^+^u',v^+^v',w^+^w')
-  negateV (u,v,w)        = (negateV u, negateV v, negateV w)
+  s *^ (u,v,w) = (s*^u,s*^v,s*^w)
 
 instance (InnerSpace u s,InnerSpace v s,InnerSpace w s, VectorSpace s s')
     => InnerSpace (u,v,w) s where
   (u,v,w) <.> (u',v',w') = u<.>u' ^+^ v<.>v' ^+^ w<.>w'
 
 
--- Standard instance for an applicative functor applied to a vector space.
-instance VectorSpace v s => VectorSpace (a->v) s where
-  zeroV   = pure   zeroV
-  (*^) s  = fmap   (s *^)
-  (^+^)   = liftA2 (^+^)
-  negateV = fmap   negateV
-
--- I don't know how to make the InnerSpace class work out, because the
--- inner product would have to combine two vector *functions* into a
--- scalar value.
--- 
---   instance InnerSpace v s => InnerSpace (a->v) s where
---     (<.>) = ???
-
--- Alternatively, we could use (a->s) as the scalar field:
--- 
---   -- Standard instance for an applicative functor applied to a vector space.
---   instance VectorSpace v s => VectorSpace (a->v) (a->s) where
---     zeroV   = pure   zeroV
---     (*^)    = liftA2 (*^)
---     (^+^)   = liftA2 (^+^)
---     negateV = fmap negateV
--- 
---   instance InnerSpace v s => InnerSpace (a->v) (a->s) where
---     (<.>) = liftA2 (<.>)
--- 
--- This definition, however, doesn't fit the standard notion of linear
--- maps as vector spaces.
-
-
--- | Linear transformations/maps.  For now, represented as simple
--- functions.  The 'VectorSpace' instance for functions gives the usual
--- meaning for a vector space of linear transformations.
+-- -- Standard instance for an applicative functor applied to a vector space.
+-- instance VectorSpace v s => VectorSpace (a->v) s where
+--   (*^) s = fmap (s *^)
 
-type a :-* b = a -> b
+-- No 'InnerSpace' instance for @(a->v)@.
diff --git a/tests/src/Perf.hs b/tests/src/Perf.hs
new file mode 100644
--- /dev/null
+++ b/tests/src/Perf.hs
@@ -0,0 +1,197 @@
+{-# LANGUAGE TypeOperators, MultiParamTypeClasses, UndecidableInstances, FlexibleInstances
+  #-}
+
+
+-- This module tests *performance* of the vector-space operations, such that it is possible to catch performance regressions.
+
+
+module Main where
+
+import Control.Applicative
+import System.Time
+import Data.List
+
+import Data.NumInstances ()
+import Data.VectorSpace
+import Data.Cross
+import Data.Derivative
+import Data.LinearMap
+
+type Surf s        = (s,s) -> (s,s,s)
+type HeightField s = (s,s) -> s
+type Curve2 s      = s -> (s,s)
+
+type Warp1 s        = s -> s
+type Warp2 s        = (s,s) -> (s,s)
+type Warp3 s        = (s,s,s) -> (s,s,s)
+
+type R = Double
+
+cosU, sinU :: Floating s => s -> s
+cosU = cos . mul2pi
+sinU = sin . mul2pi
+
+mul2pi :: Floating s => s -> s
+mul2pi = (* (2*pi))
+
+torus :: (Floating s, VectorSpace s s) => s -> s -> Surf s
+torus sr cr = revolve (\ s -> (sr,0) ^+^ cr *^ circle s)
+
+-- Try use rules to optimize?
+-- # RULES "sphere" sphere1 = spec_sphere1
+sphere1 :: Floating s => Surf s
+sphere1 = revolve semiCircle
+
+spec_sphere1 :: Surf ((Double,Double) :> Double)
+spec_sphere1 = sphere1
+
+semiCircle :: Floating s => Curve2 s
+semiCircle = circle . (/ 2)
+
+circle :: Floating s => Curve2 s
+circle = liftA2 (,) cosU sinU
+
+revolveG :: Floating s => (s -> Curve2 s) -> Surf s
+revolveG curveF = \ (u,v) -> onXY (rotate (-2*pi*v)) (addY (curveF v) u)
+
+revolve :: Floating s => Curve2 s -> Surf s
+revolve curve = revolveG (const curve)
+
+rotate :: Floating s => s -> Warp2 s
+rotate theta = \ (x,y) -> (x * c - y * s, y * c + x * s)
+ where c = cos theta
+       s = sin theta
+
+addX, addY, addZ :: Num s => (a -> Two s) -> (a -> Three s)
+addX = fmap (\ (y,z) -> (0,y,z))
+addY = fmap (\ (x,z) -> (x,0,z))
+addZ = fmap (\ (x,y) -> (x,y,0))
+
+addYZ,addXZ,addXY :: Num s => (a -> One s) -> (a -> Three s)
+addYZ = fmap (\ x -> (x,0,0))
+addXZ = fmap (\ y -> (0,y,0))
+addXY = fmap (\ z -> (0,0,z))
+
+onX,onY,onZ :: Warp1 s -> Warp3 s
+onX f (x,y,z) = (f x, y, z)
+onY f (x,y,z) = (x, f y, z)
+onZ f (x,y,z) = (x, y, f z)
+
+onXY,onYZ,onXZ :: Warp2 s -> Warp3 s
+onXY f (x,y,z) = (x',y',z ) where (x',y') = f (x,y)
+onXZ f (x,y,z) = (x',y ,z') where (x',z') = f (x,z)
+onYZ f (x,y,z) = (x ,y',z') where (y',z') = f (y,z)
+
+
+onX',onY',onZ' :: Warp1 s -> (a -> Three s) -> (a -> Three s)
+onX' = fmap fmap onX
+onY' = fmap fmap onY
+onZ' = fmap fmap onZ
+
+onXY',onXZ',onYZ' :: Warp2 s -> (a -> Three s) -> (a -> Three s)
+onXY' = fmap fmap onXY
+onXZ' = fmap fmap onXZ
+onYZ' = fmap fmap onYZ
+
+displace :: (InnerSpace v s, Floating s, HasNormal v, Applicative f) =>
+            f v -> f s -> f v
+displace = liftA2 displaceV
+
+displaceV :: (InnerSpace v s, Floating s, HasNormal v) =>
+             v -> s -> v
+displaceV v s = v ^+^ s *^ normal v
+
+------------------------------------------------------------------------------
+
+surfs3 :: [(Surf ((Double,Double) :> Double),String)]
+surfs3 = [ (displace surf hmap,m1 ++ " `displace` " ++ m2) 
+	 | (surf,m1) <- surfs2
+	 , (hmap,m2) <- hmaps
+	 ]
+
+surfs2 :: [(Surf ((Double,Double) :> Double),String)]
+surfs2 = [ (displace surf hmap,m1 ++ " `displace` " ++ m2) 
+	 | (surf,m1) <- surfs
+	 , (hmap,m2) <- hmaps
+	 ]
+
+surfs :: [(Surf ((Double,Double) :> Double),String)]
+surfs =
+  [ (torus 1 (1/2) ,"torus")
+  , (sphere1,"sphere")
+  ]
+
+hmaps :: [(HeightField ((Double,Double) :> Double),String)]
+hmaps = 
+  [ (\ (_,_) -> 0,"flat")
+  , (\ (u,v) -> cosU u * sinU v,"eggcrate")
+  ]
+
+main :: IO ()
+main = do 
+	let loop msg fun t count (points:pss) = do
+		sequence_ [ p1 `seq` p2 `seq` p3 `seq` n1 `seq` n2 `seq` n3 `seq` return ()
+        	          | (x,y) <- points
+		          , let ((p1,p2,p3),(n1,n2,n3)) = vsurf fun (x,y) ]
+		diff <- currRelTime t
+--		print diff
+		if diff > 2
+		  then do let count' = count + length points
+			  putStrLn $ "Sample count rate for " ++ msg ++ " is " ++ show (fromIntegral count' / diff) ++ " (total count = " ++ show count' ++ ")"
+			  return ()
+		  else loop msg fun t (count + length points) pss
+	    loop _ _ _ _ _ = return ()
+
+	sequence_ [ do t <- getClockTime
+		       loop msg fun t 0 (progressive_filter samples_2d)
+		  | (fun,msg) <- concat [ surfs, surfs, surfs, surfs2, surfs3 ]
+	 	  ]
+
+currRelTime :: ClockTime -> IO Double
+currRelTime (TOD sec0 pico0) = fmap delta getClockTime
+ where
+   delta (TOD sec pico) =
+     fromIntegral (sec-sec0) + 1.0e-12 * fromIntegral (pico-pico0)
+
+------------------------------------------------------------------------------
+
+vsurf :: Surf ((R,R) :> R) -> (R,R) -> ((R,R,R),(R,R,R))
+vsurf surf = toVN3 . vector3D . surf . unvector2D . idD
+
+type SurfPt s = (s,s) :> (s,s,s)
+
+toVN3 :: (LMapDom s s,Floating s, InnerSpace s s) => SurfPt s -> ((s,s,s),(s,s,s))
+toVN3 v = ( powVal v
+	  , powVal (normal v)
+	  )
+vector3D :: (LMapDom a s,VectorSpace s s) => (a :> s,a :> s,a :> s) -> (a :> (s,s,s))
+vector3D (u,v,w) = liftD3 (,,) u v w
+unvector2D :: (Data.LinearMap.LMapDom a s,VectorSpace s s) => (a :> (s,s)) -> (a :> s,a :> s) 
+unvector2D d = ( (\ (x,_) -> x) <$>> d
+	       , (\ (_,y) -> y) <$>> d
+	       )
+
+------------------------------------------------------------------------------
+
+between :: [Double] -> [Double]
+between xs = [ (n + m) / 2 | (n,m) <- zip xs (tail xs) ]
+
+samples_1d :: [[Double]]
+samples_1d = fn [0,1]
+     where
+	fn :: [Double] -> [[Double]]
+	fn points = points : fn (sort (points ++ between points))
+
+samples_2d :: [[(Double,Double)]]
+samples_2d =  [ [ (a,b) 
+		| a <- sam
+		, b <- sam
+		]
+  	      | sam <- samples_1d
+	      ]
+
+-- only allows new points through.
+progressive_filter :: (Ord a) => [[a]] -> [[a]]
+progressive_filter xs = head sorted_xs : [ y \\ x | (x,y) <- zip sorted_xs (tail sorted_xs) ]
+  where
+	sorted_xs = map sort xs
diff --git a/vector-space.cabal b/vector-space.cabal
--- a/vector-space.cabal
+++ b/vector-space.cabal
@@ -1,6 +1,7 @@
 Name:                vector-space
-Version:             0.1.3
-Synopsis: 	     Vector & affine spaces, plus derivatives
+Version:             0.2.0
+Cabal-Version:       >= 1.2
+Synopsis:            Vector & affine spaces, plus derivatives
 Category:            math
 Description:
   vector-space provides classes and generic operations for vector
@@ -18,18 +19,35 @@
 Author:              Conal Elliott 
 Maintainer:          conal@conal.net
 Homepage:            http://haskell.org/haskellwiki/vector-space
-Package-Url:	     http://code.haskell.org/vector-space
+Package-Url:         http://code.haskell.org/vector-space
 Copyright:           (c) 2007-2008 by Conal Elliott
 License:             BSD3
 Stability:           experimental
-build-type:	     Simple
-Hs-Source-Dirs:      src
-Extensions:          
-Build-Depends:       base
-Exposed-Modules:     
-		     Data.VectorSpace
-		     Data.Derivative
-		     Data.AffineSpace
-		     Data.NumInstances
-		     
-ghc-options:         -Wall -O2
+build-type:          Simple
+Build-Depends:       base, OpenGL
+
+
+Library
+  Hs-Source-Dirs:      src
+  Extensions:          
+  Exposed-Modules:     
+                     Data.AdditiveGroup
+                     Data.VectorSpace
+                     Data.LinearMap
+                     Data.Maclaurin
+                     Data.Derivative
+                     Data.Cross
+                     Data.AffineSpace
+                     Data.NumInstances
+  ghc-options:         -Wall -O2
+  ghc-prof-options:    -prof -auto-all 
+
+Executable Perf
+  main-is:           Perf.hs
+  build-depends:     base, OpenGL, old-time
+  Hs-Source-Dirs:    src, tests/src
+  ghc-options:       -Wall -O2
+  ghc-prof-options:    -prof -auto-all   
+
+--                   Data.Horner
+-- The OpenGL dependency is a temporary workaround for a ghc-6.8.2 type family bug.