diff --git a/COPYING b/COPYING
new file mode 100644
--- /dev/null
+++ b/COPYING
@@ -0,0 +1,25 @@
+Copyright (c) 2009 Conal Elliott
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+1. Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+3. The names of the authors may not be used to endorse or promote products
+   derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
+IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
diff --git a/Makefile b/Makefile
deleted file mode 100644
--- a/Makefile
+++ /dev/null
@@ -1,9 +0,0 @@
-# For special configuration, especially for docs.  Otherwise see README.
-
-server = code.haskell.org
-server-dir = /srv/code
-server-url-dir =
-
-# extra-configure-args += --enable-library-profiling --enable-executable-profiling
-
-include ../my-cabal-make.inc
diff --git a/src/Data/AdditiveGroup.hs b/src/Data/AdditiveGroup.hs
--- a/src/Data/AdditiveGroup.hs
+++ b/src/Data/AdditiveGroup.hs
@@ -1,47 +1,74 @@
-{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE TypeOperators, CPP #-}
+{-# LANGUAGE FlexibleInstances  #-}
+{-# LANGUAGE FlexibleContexts   #-}
+{-# LANGUAGE DefaultSignatures   #-}
+{-# LANGUAGE ScopedTypeVariables #-}
 ----------------------------------------------------------------------
 -- |
 -- 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(..), (^-^), sumV
+  (
+    AdditiveGroup(..), sumV
   , Sum(..), inSum, inSum2
   ) where
 
+import Prelude hiding (foldr)
+
 import Control.Applicative
+#if !(MIN_VERSION_base(4,8,0))
 import Data.Monoid (Monoid(..))
+import Data.Foldable (Foldable)
+#endif
+import Data.Foldable (foldr)
 import Data.Complex hiding (magnitude)
+import Data.Ratio
+#if !(MIN_VERSION_base(4,11,0))
+import Data.Semigroup (Semigroup(..))
+#endif
+import Foreign.C.Types (CSChar, CInt, CShort, CLong, CLLong, CIntMax, CFloat, CDouble)
 
 import Data.MemoTrie
 
+import Data.VectorSpace.Generic
+import qualified GHC.Generics as Gnrx
+import GHC.Generics (Generic, (:*:)(..))
+
 infixl 6 ^+^, ^-^
 
 -- | Additive group @v@.
 class AdditiveGroup v where
   -- | The zero element: identity for '(^+^)'
   zeroV :: v
+  default zeroV :: (Generic v, AdditiveGroup (VRep v)) => v
+  zeroV = Gnrx.to (zeroV :: VRep v)
+  {-# INLINE zeroV #-}
   -- | Add vectors
   (^+^) :: v -> v -> v
+  default (^+^) :: (Generic v, AdditiveGroup (VRep v)) => v -> v -> v
+  v ^+^ v' = Gnrx.to (Gnrx.from v ^+^ Gnrx.from v' :: VRep v)
+  {-# INLINE (^+^) #-}
   -- | Additive inverse
   negateV :: v -> v
-
--- | Group subtraction
-(^-^) :: AdditiveGroup v => v -> v -> v
-v ^-^ v' = v ^+^ negateV v'
+  default negateV :: (Generic v, AdditiveGroup (VRep v)) => v -> v
+  negateV v = Gnrx.to (negateV $ Gnrx.from v :: VRep v)
+  {-# INLINE negateV #-}
+  -- | Group subtraction
+  (^-^) :: v -> v -> v
+  v ^-^ v' = v ^+^ negateV v'
 
 -- | Sum over several vectors
-sumV :: AdditiveGroup v => [v] -> v
+sumV :: (Foldable f, AdditiveGroup v) => f v -> v
 sumV = foldr (^+^) zeroV
-
+{-# INLINE sumV #-}
 
 instance AdditiveGroup () where
   zeroV     = ()
@@ -49,15 +76,28 @@
   negateV   = id
 
 -- For 'Num' types:
--- 
+--
 -- instance AdditiveGroup n where {zeroV=0; (^+^) = (+); negateV = negate}
 
-instance AdditiveGroup Int     where {zeroV=0; (^+^) = (+); negateV = negate}
-instance AdditiveGroup Integer where {zeroV=0; (^+^) = (+); negateV = negate}
-instance AdditiveGroup Float   where {zeroV=0; (^+^) = (+); negateV = negate}
-instance AdditiveGroup Double  where {zeroV=0; (^+^) = (+); negateV = negate}
+#define ScalarTypeCon(con,t) \
+  instance con => AdditiveGroup (t) where {zeroV=0; (^+^) = (+); negateV = negate}
 
+#define ScalarType(t) ScalarTypeCon((),t)
 
+ScalarType(Int)
+ScalarType(Integer)
+ScalarType(Float)
+ScalarType(Double)
+ScalarType(CSChar)
+ScalarType(CInt)
+ScalarType(CShort)
+ScalarType(CLong)
+ScalarType(CLLong)
+ScalarType(CIntMax)
+ScalarType(CFloat)
+ScalarType(CDouble)
+ScalarTypeCon(Integral a,Ratio a)
+
 instance (RealFloat v, AdditiveGroup v) => AdditiveGroup (Complex v) where
   zeroV   = zeroV :+ zeroV
   (^+^)   = (+)
@@ -100,7 +140,37 @@
   Just a' ^+^ Just b' = Just (a' ^+^ b')
   negateV = fmap negateV
 
+{-
 
+Alexey Khudyakov wrote:
+
+  I looked through vector-space package and found lawless instance. Namely Maybe's AdditiveGroup instance
+
+  It's group so following relation is expected to hold. Otherwise it's not a group.
+  > x ^+^ negateV x == zeroV
+
+  Here is counterexample:
+
+  > let x = Just 2 in x ^+^ negateV x == zeroV
+  False
+
+  I think it's not possible to sensibly define group instance for
+  Maybe a at all.
+
+
+I see that the problem here is in distinguishing 'Just zeroV' from
+Nothing. I could fix the Just + Just line to use Nothing instead of Just
+zeroV when a' ^+^ b' == zeroV, although doing so would require Eq a and
+hence lose some generality. Even so, the abstraction leak would probably
+show up elsewhere.
+
+Hm.
+
+-}
+
+
+
+
 -- Memo tries
 instance (HasTrie u, AdditiveGroup v) => AdditiveGroup (u :->: v) where
   zeroV   = pure   zeroV
@@ -115,6 +185,7 @@
 
 instance Functor Sum where
   fmap f (Sum a) = Sum (f a)
+  {-# INLINE fmap #-}
 
 -- instance Applicative Sum where
 --   pure a = Sum a
@@ -122,32 +193,40 @@
 
 instance Applicative Sum where
   pure  = Sum
+  {-# INLINE pure #-}
   (<*>) = inSum2 ($)
+  {-# INLINE (<*>) #-}
 
+instance AdditiveGroup a => Semigroup (Sum a) where
+  (<>) = liftA2 (^+^)
+  {-# INLINE (<>) #-}
+
 instance AdditiveGroup a => Monoid (Sum a) where
   mempty  = Sum zeroV
-  mappend = liftA2 (^+^)
-
+#if !(MIN_VERSION_base(4,11,0))
+  mappend = (<>)
+#endif
 
 -- | Application a unary function inside a 'Sum'
 inSum :: (a -> b) -> (Sum a -> Sum b)
 inSum = getSum ~> Sum
+{-# INLINE inSum #-}
 
 -- | Application a binary function inside a 'Sum'
 inSum2 :: (a -> b -> c) -> (Sum a -> Sum b -> Sum c)
 inSum2 = getSum ~> inSum
-
+{-# INLINE inSum2 #-}
 
 instance AdditiveGroup a => AdditiveGroup (Sum a) where
-  zeroV   = mempty
-  (^+^)   = mappend
+  zeroV   = Sum zeroV
+  (^+^)   = (<>)
   negateV = inSum negateV
 
-
 ---- to go elsewhere
 
 (~>) :: (a' -> a) -> (b -> b') -> ((a -> b) -> (a' -> b'))
 (i ~> o) f = o . f . i
+{-# INLINE (~>) #-}
 
 -- result :: (b -> b') -> ((a -> b) -> (a -> b'))
 -- result = (.)
@@ -156,3 +235,33 @@
 -- argument = flip (.)
 
 -- g ~> f = result g . argument f
+
+
+
+instance AdditiveGroup a => AdditiveGroup (Gnrx.Rec0 a s) where
+  zeroV = Gnrx.K1 zeroV
+  {-# INLINE zeroV #-}
+  negateV (Gnrx.K1 v) = Gnrx.K1 $ negateV v
+  {-# INLINE negateV #-}
+  Gnrx.K1 v ^+^ Gnrx.K1 w = Gnrx.K1 $ v ^+^ w
+  {-# INLINE (^+^) #-}
+  Gnrx.K1 v ^-^ Gnrx.K1 w = Gnrx.K1 $ v ^-^ w
+  {-# INLINE (^-^) #-}
+instance AdditiveGroup (f p) => AdditiveGroup (Gnrx.M1 i c f p) where
+  zeroV = Gnrx.M1 zeroV
+  {-# INLINE zeroV #-}
+  negateV (Gnrx.M1 v) = Gnrx.M1 $ negateV v
+  {-# INLINE negateV #-}
+  Gnrx.M1 v ^+^ Gnrx.M1 w = Gnrx.M1 $ v ^+^ w
+  {-# INLINE (^+^) #-}
+  Gnrx.M1 v ^-^ Gnrx.M1 w = Gnrx.M1 $ v ^-^ w
+  {-# INLINE (^-^) #-}
+instance (AdditiveGroup (f p), AdditiveGroup (g p)) => AdditiveGroup ((f :*: g) p) where
+  zeroV = zeroV :*: zeroV
+  {-# INLINE zeroV #-}
+  negateV (x:*:y) = negateV x :*: negateV y
+  {-# INLINE negateV #-}
+  (x:*:y) ^+^ (ξ:*:υ) = (x^+^ξ) :*: (y^+^υ)
+  {-# INLINE (^+^) #-}
+  (x:*:y) ^-^ (ξ:*:υ) = (x^-^ξ) :*: (y^-^υ)
+  {-# INLINE (^-^) #-}
diff --git a/src/Data/AffineSpace.hs b/src/Data/AffineSpace.hs
--- a/src/Data/AffineSpace.hs
+++ b/src/Data/AffineSpace.hs
@@ -1,4 +1,10 @@
-{-# LANGUAGE MultiParamTypeClasses, FlexibleContexts, TypeFamilies #-}
+{-# LANGUAGE MultiParamTypeClasses, FlexibleContexts, TypeFamilies, CPP #-}
+{-# LANGUAGE FlexibleInstances  #-}
+{-# LANGUAGE DefaultSignatures   #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeOperators       #-}
+{-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE DeriveGeneric        #-}
 ----------------------------------------------------------------------
 -- |
 -- Module      :  Data.AffineSpace
@@ -13,25 +19,51 @@
 
 module Data.AffineSpace
   (
-    AffineSpace(..), (.-^), distanceSq, distance, alerp
+    AffineSpace(..), (.-^), distanceSq, distance, alerp, affineCombo
   ) where
-
+#if !MIN_VERSION_base(4,10,0)
 import Control.Applicative (liftA2)
+#endif
+import Data.Ratio
+import Foreign.C.Types (CSChar, CInt, CShort, CLong, CLLong, CIntMax, CFloat, CDouble)
+import Control.Arrow(first)
 
 import Data.VectorSpace
+import Data.Basis
 
-infix 4 .+^, .-^, .-.
+import Data.VectorSpace.Generic
+import qualified GHC.Generics as Gnrx
+import GHC.Generics (Generic, (:*:)(..))
 
+-- Through 0.8.4, I used the following fixities.
+-- 
+--   infix 4 .+^, .-^, .-.
+-- 
+-- Changed in 0.8.5 to match precedence of + and -, and to associate usefully.
+-- Thanks to Ben Gamari for suggesting left-associativity.
+
+infixl 6 .+^, .-^
+infix  6 .-.
+
+
 -- TODO: Convert AffineSpace from fundep to associated type, and eliminate
 -- FunctionalDependencies above.
 
 class AdditiveGroup (Diff p) => AffineSpace p where
   -- | Associated vector space
   type Diff p
+  type Diff p = GenericDiff p
   -- | Subtract points
   (.-.)  :: p -> p -> Diff p
+  default (.-.) :: ( Generic p, Diff p ~ GenericDiff p, AffineSpace (VRep p) )
+              => p -> p -> Diff p
+  p .-. q = GenericDiff
+         $ (Gnrx.from p .-. (Gnrx.from q :: VRep p))
   -- | Point plus vector
   (.+^)  :: p -> Diff p -> p
+  default (.+^) :: ( Generic p, Diff p ~ GenericDiff p, AffineSpace (VRep p) )
+              => p -> Diff p -> p
+  p .+^ GenericDiff q = Gnrx.to (Gnrx.from p .+^ q :: VRep p)
 
 -- | Point minus vector
 (.-^) :: AffineSpace p => p -> Diff p -> p
@@ -55,17 +87,46 @@
          p -> p -> Scalar (Diff p) -> p
 alerp p p' s = p .+^ (s *^ (p' .-. p))
 
-instance  AffineSpace Double where
-  type Diff Double = Double
-  (.-.) =  (-)
-  (.+^) =  (+)
+-- | Compute an affine combination (weighted average) of points.
+-- The first element is used as origin and is weighted
+-- such that all coefficients sum to 1. For example,
+--
+-- > affineCombo a [(0.3,b), (0.2,c)]
+--
+-- is equal to
+--
+-- > a .+^ (0.3 *^ (b .-. a) ^+^ 0.2 *^ (c .-. a))
+--
+-- and if @a@, @b@, and @c@ were in a vector space would also be equal to
+--
+-- > 0.5 *^ a ^+^ 0.3 *^ b ^+^ 0.2 *^ c
+--
+-- See also 'linearCombo' (on vector spaces).
+affineCombo :: (AffineSpace p, v ~ Diff p, VectorSpace v) => p -> [(p,Scalar v)] -> p
+affineCombo z l = z .+^ linearCombo (map (first (.-. z)) l)
 
-instance  AffineSpace Float where
-  type Diff Float = Float
-  (.-.) =  (-)
-  (.+^) =  (+)
+#define ScalarTypeCon(con,t) \
+  instance con => AffineSpace (t) where \
+    { type Diff (t) = t \
+    ; (.-.) = (-) \
+    ; (.+^) = (+) }
 
+#define ScalarType(t) ScalarTypeCon((),t)
 
+ScalarType(Int)
+ScalarType(Integer)
+ScalarType(Double)
+ScalarType(Float)
+ScalarType(CSChar)
+ScalarType(CInt)
+ScalarType(CShort)
+ScalarType(CLong)
+ScalarType(CLLong)
+ScalarType(CIntMax)
+ScalarType(CDouble)
+ScalarType(CFloat)
+ScalarTypeCon(Integral a,Ratio a)
+
 instance (AffineSpace p, AffineSpace q) => AffineSpace (p,q) where
   type Diff (p,q)   = (Diff p, Diff q)
   (p,q) .-. (p',q') = (p .-. p', q .-. q')
@@ -81,3 +142,59 @@
   type Diff (a -> p) = a -> Diff p
   (.-.)              = liftA2 (.-.)
   (.+^)              = liftA2 (.+^)
+
+
+
+newtype GenericDiff p = GenericDiff (Diff (VRep p))
+       deriving (Generic)
+
+instance AdditiveGroup (Diff (VRep p)) => AdditiveGroup (GenericDiff p)
+instance VectorSpace (Diff (VRep p)) => VectorSpace (GenericDiff p)
+instance (AdditiveGroup (Scalar (Diff (VRep p))), InnerSpace (Diff (VRep p))) => InnerSpace (GenericDiff p)
+instance HasBasis (Diff (VRep p)) => HasBasis (GenericDiff p)
+
+data AffineDiffProductSpace f g p = AffineDiffProductSpace
+            !(Diff (f p)) !(Diff (g p)) deriving (Generic)
+instance (AffineSpace (f p), AffineSpace (g p))
+    => AdditiveGroup (AffineDiffProductSpace f g p)
+instance ( AffineSpace (f p), AffineSpace (g p)
+         , VectorSpace (Diff (f p)), VectorSpace (Diff (g p))
+         , Scalar (Diff (f p)) ~ Scalar (Diff (g p)) )
+    => VectorSpace (AffineDiffProductSpace f g p)
+instance ( AdditiveGroup (Scalar (Diff (g p)))
+         , AffineSpace (f p), AffineSpace (g p)
+         , InnerSpace (Diff (f p)), InnerSpace (Diff (g p))
+         , Scalar (Diff (f p)) ~ Scalar (Diff (g p))
+         , Num (Scalar (Diff (f p))) )
+    => InnerSpace (AffineDiffProductSpace f g p)
+instance (AffineSpace (f p), AffineSpace (g p))
+    => AffineSpace (AffineDiffProductSpace f g p) where
+  type Diff (AffineDiffProductSpace f g p) = AffineDiffProductSpace f g p
+  (.+^) = (^+^)
+  (.-.) = (^-^)
+instance ( AffineSpace (f p), AffineSpace (g p)
+         , HasBasis (Diff (f p)), HasBasis (Diff (g p))
+         , Scalar (Diff (f p)) ~ Scalar (Diff (g p)) )
+    => HasBasis (AffineDiffProductSpace f g p) where
+  type Basis (AffineDiffProductSpace f g p) = Either (Basis (Diff (f p)))
+                                                     (Basis (Diff (g p)))
+  basisValue (Left bf) = AffineDiffProductSpace (basisValue bf) zeroV
+  basisValue (Right bg) = AffineDiffProductSpace zeroV (basisValue bg)
+  decompose (AffineDiffProductSpace vf vg)
+        = map (first Left) (decompose vf) ++ map (first Right) (decompose vg)
+  decompose' (AffineDiffProductSpace vf _) (Left bf) = decompose' vf bf
+  decompose' (AffineDiffProductSpace _ vg) (Right bg) = decompose' vg bg
+
+
+instance AffineSpace a => AffineSpace (Gnrx.Rec0 a s) where
+  type Diff (Gnrx.Rec0 a s) = Diff a
+  Gnrx.K1 v .+^ w = Gnrx.K1 $ v .+^ w
+  Gnrx.K1 v .-. Gnrx.K1 w = v .-. w
+instance AffineSpace (f p) => AffineSpace (Gnrx.M1 i c f p) where
+  type Diff (Gnrx.M1 i c f p) = Diff (f p)
+  Gnrx.M1 v .+^ w = Gnrx.M1 $ v .+^ w
+  Gnrx.M1 v .-. Gnrx.M1 w = v .-. w
+instance (AffineSpace (f p), AffineSpace (g p)) => AffineSpace ((f :*: g) p) where
+  type Diff ((f:*:g) p) = AffineDiffProductSpace f g p
+  (x:*:y) .+^ AffineDiffProductSpace ξ υ = (x.+^ξ) :*: (y.+^υ)
+  (x:*:y) .-. (ξ:*:υ) = AffineDiffProductSpace (x.-.ξ) (y.-.υ)
diff --git a/src/Data/Basis.hs b/src/Data/Basis.hs
--- a/src/Data/Basis.hs
+++ b/src/Data/Basis.hs
@@ -1,10 +1,8 @@
--- WARNING: this module depends on type families working fairly well, and
--- requires ghc version at least 6.9.  I didn't find a way to specify that
--- dependency in the .cabal.
--- 
 {-# LANGUAGE TypeOperators, TypeFamilies, UndecidableInstances
-  , FlexibleInstances, MultiParamTypeClasses
-  #-}
+  , FlexibleInstances, MultiParamTypeClasses, CPP  #-}
+{-# LANGUAGE DefaultSignatures    #-}
+{-# LANGUAGE FlexibleContexts     #-}
+{-# LANGUAGE ScopedTypeVariables  #-}
 {-# OPTIONS_GHC -Wall -fno-warn-orphans #-}
 ----------------------------------------------------------------------
 -- |
@@ -23,37 +21,51 @@
 
 -- import Control.Applicative ((<$>))
 import Control.Arrow (first)
-import Data.Either
+import Data.Ratio
+import Foreign.C.Types (CFloat, CDouble)
+import Data.Kind
+-- import Data.Either
 
 import Data.VectorSpace
 
+import Data.VectorSpace.Generic
+import qualified GHC.Generics as Gnrx
+import GHC.Generics (Generic, (:*:)(..))
+
 -- using associated data type instead of associated type synonym to work
 -- around ghc bug <http://hackage.haskell.org/trac/ghc/ticket/3038>
 
 class VectorSpace v => HasBasis v where
   -- | Representation of the canonical basis for @v@
-  type Basis v :: *
+  type Basis v :: Type
+  type Basis v = Basis (VRep v)
   -- | Interpret basis rep as a vector
   basisValue   :: Basis v -> v
+  default basisValue :: (Generic v, HasBasis (VRep v), Basis (VRep v) ~ Basis v)
+                    => Basis v -> v
+  basisValue b = Gnrx.to (basisValue b :: VRep v)
   -- | Extract coordinates
   decompose    :: v -> [(Basis v, Scalar v)]
+  default decompose :: ( Generic v, HasBasis (VRep v)
+                       , Scalar (VRep v) ~ Scalar v, Basis (VRep v) ~ Basis v )
+                    => v -> [(Basis v, Scalar v)]
+  decompose v = decompose (Gnrx.from v :: VRep v)
   -- | Experimental version.  More elegant definitions, and friendly to
   -- infinite-dimensional vector spaces.
   decompose'   :: v -> (Basis v -> Scalar v)
+  default decompose' :: ( Generic v, HasBasis (VRep v)
+                        , Scalar (VRep v) ~ Scalar v, Basis (VRep v) ~ Basis v )
+                    => v -> Basis v -> Scalar v
+  decompose' v = decompose' (Gnrx.from v :: VRep v)
 
 -- Defining property: recompose . decompose == id
 
--- | Linear combination
-linearCombo :: VectorSpace v => [(v,Scalar v)] -> v
-linearCombo ps = sumV [s *^ v | (v,s) <- ps]
-
 -- Turn a basis decomposition back into a vector.
 recompose :: HasBasis v => [(Basis v, Scalar v)] -> v
 recompose = linearCombo . fmap (first basisValue)
 
 -- recompose ps = linearCombo (first basisValue <$> ps)
 
-
 -- I don't know how to define
 -- 
 --   recompose' :: HasBasis v => (Basis v -> Scalar v) -> v
@@ -61,18 +73,21 @@
 -- However, I don't seem to use recompose anywhere.
 -- I don't even use basisValue or decompose.
 
-instance HasBasis Float where
-  type Basis Float = ()
-  basisValue ()    = 1
-  decompose s      = [((),s)]
-  decompose' s     = const s
+#define ScalarTypeCon(con,t) \
+  instance con => HasBasis (t) where \
+    { type Basis (t) = () \
+    ; basisValue ()  = 1 \
+    ; decompose s    = [((),s)] \
+    ; decompose' s   = const s }
 
-instance HasBasis Double where
-  type Basis Double = ()
-  basisValue ()     = 1
-  decompose s       = [((),s)]
-  decompose' s      = const s
+#define ScalarType(t) ScalarTypeCon((),t)
 
+ScalarType(Float)
+ScalarType(CFloat)
+ScalarType(Double)
+ScalarType(CDouble)
+ScalarTypeCon(Integral a, Ratio a)
+
 instance ( HasBasis u, s ~ Scalar u
          , HasBasis v, s ~ Scalar v )
       => HasBasis (u,v) where
@@ -136,3 +151,21 @@
 t4 = basisValue (Right (Left ())) :: (Float,Double,Float)
 
 -}
+
+instance HasBasis a => HasBasis (Gnrx.Rec0 a s) where
+  type Basis (Gnrx.Rec0 a s) = Basis a
+  basisValue = Gnrx.K1 . basisValue
+  decompose = decompose . Gnrx.unK1
+  decompose' = decompose' . Gnrx.unK1
+instance HasBasis (f p) => HasBasis (Gnrx.M1 i c f p) where
+  type Basis (Gnrx.M1 i c f p) = Basis (f p)
+  basisValue = Gnrx.M1 . basisValue
+  decompose = decompose . Gnrx.unM1
+  decompose' = decompose' . Gnrx.unM1
+instance (HasBasis (f p), HasBasis (g p), Scalar (f p) ~ Scalar (g p))
+         => HasBasis ((f :*: g) p) where
+  type Basis ((f:*:g) p) = Either (Basis (f p)) (Basis (g p))
+  basisValue (Left bf) = basisValue bf :*: zeroV
+  basisValue (Right bg) = zeroV :*: basisValue bg
+  decompose  (u:*:v)     = decomp2 Left u ++ decomp2 Right v
+  decompose' (u:*:v)     = decompose' u `either` decompose' v
diff --git a/src/Data/Cross.hs b/src/Data/Cross.hs
--- a/src/Data/Cross.hs
+++ b/src/Data/Cross.hs
@@ -1,6 +1,6 @@
 {-# LANGUAGE FlexibleInstances, FlexibleContexts, TypeOperators
-           , TypeFamilies, TypeSynonymInstances
-  #-}
+           , TypeFamilies, TypeSynonymInstances 
+           , UndecidableInstances  #-}
 {-# OPTIONS_GHC -Wall #-}
 ----------------------------------------------------------------------
 -- |
@@ -49,8 +49,7 @@
 instance AdditiveGroup u => HasCross2 (u,u) where
   cross2 (x,y) = (negateV y,x)  -- or @(y,-x)@?
 
-instance ( HasBasis a, HasTrie (Basis a)
-         , VectorSpace v, HasCross2 v) => HasCross2 (a:>v) where
+instance (HasTrie (Basis a), HasCross2 v) => HasCross2 (a:>v) where
   -- 2d cross-product is linear
   cross2 = fmapD cross2
 
@@ -74,8 +73,7 @@
 -- l `atB` b = maybe zeroV (`untrie` b) l
 
 
-instance ( Num s, VectorSpace s
-         , HasBasis s, HasTrie (Basis s), Basis s ~ ())
+instance (VectorSpace s, HasBasis s, HasTrie (Basis s), Basis s ~ ())
     => HasNormal (Two (One s :> s)) where
   normalVec = unpairD . normalVec . pairD
 
@@ -102,7 +100,7 @@
    where
      d = derivAtBasis v
 
-instance ( Num s, VectorSpace s, HasBasis s, HasTrie (Basis s)
-         , HasNormal (Two s :> Three s))
+instance ( VectorSpace s, HasBasis s, HasTrie (Basis s)
+         , HasNormal (Two s :> Three s) )
          => HasNormal (Three (Two s :> s)) where
   normalVec = untripleD . normalVec . tripleD
diff --git a/src/Data/Horner.hs b/src/Data/Horner.hs
deleted file mode 100644
--- a/src/Data/Horner.hs
+++ /dev/null
@@ -1,220 +0,0 @@
-{-# LANGUAGE TypeOperators, MultiParamTypeClasses
-           , TypeSynonymInstances, FlexibleInstances
-  #-}
-{-# 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
--- a/src/Data/LinearMap.hs
+++ b/src/Data/LinearMap.hs
@@ -1,34 +1,40 @@
-{-# LANGUAGE TypeOperators, FlexibleContexts, TypeFamilies #-}
+{-# LANGUAGE TupleSections #-}
+{-# LANGUAGE CPP, TypeOperators, FlexibleContexts, TypeFamilies
+  , GeneralizedNewtypeDeriving, StandaloneDeriving, UndecidableInstances #-}
 {-# OPTIONS_GHC -Wall -fno-warn-orphans #-}
--- {-# OPTIONS_GHC -funbox-strict-fields #-}
--- {-# OPTIONS_GHC -ddump-simpl-stats -ddump-simpl #-}
 ----------------------------------------------------------------------
 -- |
 -- Module      :  Data.LinearMap
--- Copyright   :  (c) Conal Elliott 2008
+-- Copyright   :  (c) Conal Elliott 2008-2016
 -- License     :  BSD3
--- 
+--
 -- Maintainer  :  conal@conal.net
 -- Stability   :  experimental
--- 
+--
 -- Linear maps
 ----------------------------------------------------------------------
 
 module Data.LinearMap
-  ( (:-*) , linear, lapply, atBasis, idL, (*.*)
-  , liftMS, liftMS2, liftMS3
-  , liftL, liftL2, liftL3
-  ) where
+   ( (:-*) , linear, lapply, atBasis, idL, (*.*)
+   , inLMap, inLMap2, inLMap3
+   , liftMS, liftMS2, liftMS3
+   , liftL, liftL2, liftL3
+   , exlL, exrL, forkL, firstL, secondL
+   , inlL, inrL, joinL -- , leftL, rightL
+   )
+  where
 
-import Control.Applicative ((<$>),Applicative,liftA2,liftA3)
-import Control.Arrow       (first)
+#if !(MIN_VERSION_base(4,8,0))
+import Control.Applicative (Applicative, liftA2)
+#endif
+import Control.Applicative (liftA3)
+import Control.Arrow       (first,second)
 
-import Data.MemoTrie      ((:->:)(..))
-import Data.AdditiveGroup (Sum(..),inSum2, AdditiveGroup(..))
+import Data.MemoTrie      (HasTrie(..),(:->:))
+import Data.AdditiveGroup (Sum(..), AdditiveGroup(..))
 import Data.VectorSpace   (VectorSpace(..))
 import Data.Basis         (HasBasis(..), linearCombo)
 
-
 -- Linear maps are almost but not quite a Control.Category.  The type
 -- class constraints interfere.  They're almost an Arrow also, but for the
 -- constraints and the generality of arr.
@@ -36,35 +42,97 @@
 -- | An optional additive value
 type MSum a = Maybe (Sum a)
 
--- nsum :: MSum a
--- nsum = Nothing
-
 jsum :: a -> MSum a
 jsum = Just . Sum
 
+type LMap' u v = MSum (Basis u :->: v)
+
+infixr 1 :-*
 -- | Linear map, represented as an optional memo-trie from basis to
 -- values, where 'Nothing' means the zero map (an optimization).
-type u :-* v = MSum (Basis u :->: v)
+newtype u :-* v = LMap { unLMap :: LMap' u v }
 
--- TODO: Try a partial trie instead, excluding (known) zero elements.
--- Then 'lapply' could be much faster for sparse situations.  Make sure to
--- correctly sum them.  It'd be more like Jason Foutz's formulation
--- <http://metavar.blogspot.com/2008/02/higher-order-multivariate-automatic.html>
--- which uses in @IntMap@.
+deriving instance (HasTrie (Basis u), AdditiveGroup v) => AdditiveGroup (u :-* v)
 
+instance (HasTrie (Basis u), VectorSpace v) =>
+         VectorSpace (u :-* v) where
+  type Scalar (u :-* v) = Scalar v
+  (*^) s = (inLMap.liftMS.fmap) (s *^)
 
--- PROBLEM: u :-* v is a type synonym, and Basis is an associated type synonym, resulting in a subtle
+-- In GHC 7.10:
+-- Constraint is no smaller than the instance head
+-- in the constraint: HasTrie (Basis u)
+-- (Use UndecidableInstances to permit this)
+
+exlL :: ( HasBasis a, HasTrie (Basis a), HasBasis b, HasTrie (Basis b)
+        , Scalar a ~ Scalar b )
+     => (a,b) :-* a
+exlL = linear fst
+
+exrL :: ( HasBasis a, HasTrie (Basis a), HasBasis b, HasTrie (Basis b)
+        , Scalar a ~ Scalar b )
+     => (a,b) :-* b
+exrL = linear snd
+
+forkL :: (HasTrie (Basis a), HasBasis c, HasBasis d)
+      => (a :-* c) -> (a :-* d) -> (a :-* (c,d))
+forkL = (inLMap2.liftL2) (,)
+
+firstL  :: ( HasBasis u, HasBasis u', HasBasis v
+           , HasTrie (Basis u), HasTrie (Basis v) 
+           , Scalar u ~ Scalar v, Scalar u ~ Scalar u'
+           ) =>
+           (u :-* u') -> ((u,v) :-* (u',v))
+firstL  = linear.first.lapply
+
+secondL :: ( HasBasis u, HasBasis v, HasBasis v'
+           , HasTrie (Basis u), HasTrie (Basis v) 
+           , Scalar u ~ Scalar v, Scalar u ~ Scalar v'
+           ) =>
+           (v :-* v') -> ((u,v) :-* (u,v'))
+secondL = linear.second.lapply
+
+-- TODO: more efficient firstL
+
+-- liftMS :: (AdditiveGroup a) => (a -> b) -> (MSum a -> MSum b)
+
+-- (s *^) :: v -> v
+-- fmap (s *^) :: (Basis u :->: v) -> (Basis u :->: v)
+-- (liftMS.fmap) (s *^) :: LMap' u v -> LMap' u v
+-- (inLMap.liftMS.fmap) (s *^) :: (u :-* v) -> (u :-* v)
+
+
+inlL :: (HasBasis a, HasTrie (Basis a), HasBasis b)
+     => a :-* (a,b)
+inlL = linear (,zeroV)
+
+inrL :: (HasBasis a, HasBasis b, HasTrie (Basis b))
+     => b :-* (a,b)
+inrL = linear (zeroV,)
+
+joinL :: ( HasBasis a, HasTrie (Basis a)
+         , HasBasis b, HasTrie (Basis b)
+         , Scalar a ~ Scalar b, Scalar a ~ Scalar c
+         , VectorSpace c )
+      => (a :-* c) -> (b :-* c) -> ((a,b) :-* c)
+f `joinL` g = linear (\ (a,b) -> lapply f a ^+^ lapply g b)
+
+-- Before version 0.7, u :-* v was a type synonym, resulting in a subtle
 -- ambiguity: u:-*v == u':-*v' does not imply that u==u', since Basis
 -- might map different types to the same basis (e.g., Float & Double).
--- See <http://hackage.haskell.org/trac/ghc/ticket/1897>
--- 
--- Work in progress.  See NewLinearMap.hs
+-- See <http://hackage.haskell.org/trac/ghc/ticket/1897>.
+-- See also <http://thread.gmane.org/gmane.comp.lang.haskell.cafe/73271/focus=73332>.
 
+-- TODO: Try a partial trie instead, excluding (known) zero elements.
+-- Then 'lapply' could be much faster for sparse situations.  Make sure to
+-- correctly sum them.  It'd be more like Jason Foutz's formulation
+-- <http://metavar.blogspot.com/2008/02/higher-order-multivariate-automatic.html>
+-- which uses in @IntMap@.
 
 -- | Function (assumed linear) as linear map.
 linear :: (HasBasis u, HasTrie (Basis u)) =>
           (u -> v) -> (u :-* v)
-linear f = jsum (trie (f . basisValue))
+linear f = LMap (jsum (trie (f . basisValue)))
 
 atZ :: AdditiveGroup b => (a -> b) -> (MSum a -> b)
 atZ f = maybe zeroV (f . getSum)
@@ -72,78 +140,88 @@
 -- atZ :: AdditiveGroup b => (a -> b) -> (a -> b)
 -- atZ = id
 
--- | Evaluate a linear map on a basis element.  I've loosened the type to
--- work around a typing problem in 'derivAtBasis'.
--- atBasis :: (AdditiveGroup v, HasTrie (Basis u)) =>
---            (u :-* v) -> Basis u -> v
-atBasis :: (HasTrie a, AdditiveGroup b) => MSum (a :->: b) -> a -> b
-m `atBasis` b = atZ (`untrie` b) m
+inLMap :: (LMap' r s -> LMap' t u) -> ((r :-* s) -> (t :-* u))
+inLMap = unLMap ~> LMap
 
+inLMap2 :: (LMap' r s -> LMap' t u -> LMap' v w)
+        -> ((r :-* s) -> (t :-* u) -> (v :-* w))
+inLMap2 = unLMap ~> inLMap
+
+inLMap3 :: (LMap' r s -> LMap' t u -> LMap' v w -> LMap' x y)
+        -> ((r :-* s) -> (t :-* u) -> (v :-* w) -> (x :-* y))
+inLMap3 = unLMap ~> inLMap2
+
 -- | Apply a linear map to a vector.
 lapply :: ( VectorSpace v, Scalar u ~ Scalar v
           , HasBasis u, HasTrie (Basis u) ) =>
           (u :-* v) -> (u -> v)
-lapply = atZ lapply'
+lapply = atZ lapply' . unLMap
 
--- Handy for 'lapply' and '(*.*)'.
+-- | Evaluate a linear map on a basis element.
+atBasis :: (AdditiveGroup v, HasTrie (Basis u)) =>
+           (u :-* v) -> Basis u -> v
+LMap m `atBasis` b = atZ (`untrie` b) m
+
+-- | Handy for 'lapply' and '(*.*)'.
 lapply' :: ( VectorSpace v, Scalar u ~ Scalar v
            , HasBasis u, HasTrie (Basis u) ) =>
            (Basis u :->: v) -> (u -> v)
 lapply' tr = linearCombo . fmap (first (untrie tr)) . decompose
 
-
-
--- Identity linear map
-idL :: (HasBasis u, HasTrie (Basis u)) => 
+-- | Identity linear map
+idL :: (HasBasis u, HasTrie (Basis u)) =>
        u :-* u
 idL = linear id
 
 
 infixr 9 *.*
 -- | Compose linear maps
-(*.*) :: ( HasBasis u, HasTrie (Basis u)
+(*.*) :: ( HasTrie (Basis u)
          , HasBasis v, HasTrie (Basis v)
          , VectorSpace w
          , Scalar v ~ Scalar w ) =>
          (v :-* w) -> (u :-* v) -> (u :-* w)
 
 -- Simple definition, but only optimizes out uv == zero
--- 
--- (*.*) vw = (fmap.fmap) (lapply vw)
 
+-- vw *.* uv = LMap ((fmap.fmap.fmap) (lapply vw) (unLMap uv))
+
+(*.*) vw = (inLMap.fmap.fmap.fmap) (lapply vw)
+
+-- Eep:
+--     (*.*) = inLMap.fmap.fmap.fmap.lapply
+
+
 -- Instead, use Nothing/zero if /either/ map is zeroV (exploiting linearity
 -- when uv == zeroV.)
 
--- Nothing       *.* _             = Nothing
--- _             *.* Nothing       = Nothing
--- Just (Sum vw) *.* Just (Sum uv) = Just (Sum (lapply' vw <$> uv))
+-- LMap Nothing         *.* _                    = LMap Nothing
+-- _                    *.* LMap Nothing         = LMap Nothing
+-- LMap (Just (Sum vw)) *.* LMap (Just (Sum uv)) = LMap (Just (Sum (lapply' vw <$> uv)))
 
--- (*.*) = liftA2 (\ (Sum vw) (Sum uv) -> Sum (lapply' vw <$> uv))
+-- (*.*) = liftA2 (\ (LMap (Sum vw)) (LMap (Sum uv)) -> LMap (Sum (lapply' vw <$> uv)))
 
--- (*.*) = (liftA2.inSum2) (\ vw uv -> lapply' vw <$> uv)
-(*.*) = (liftA2.inSum2) (\ vw uv -> lapply' vw <$> uv)
+-- (*.*) = (liftA2.inSum2.inLMap2) (\ vw uv -> lapply' vw <$> uv)
 
--- (*.*) = (liftA2.inSum2) (\ vw -> fmap (lapply' vw))
+-- (*.*) = (liftA2.inSum2.inLMap2) (\ vw -> fmap (lapply' vw))
 
--- (*.*) = (liftA2.inSum2) (fmap . lapply')
+-- (*.*) = (liftA2.inSum2.inLMap2) (fmap . lapply')
 
 
 -- It may be helpful that @lapply vw@ is evaluated just once and not
 -- once per uv.  'untrie' can strip off all of its trie constructors.
 
 -- Less efficient definition:
--- 
+--
 --   vw `compL` uv = linear (lapply vw . lapply uv)
--- 
+--
 --   i.e., compL = inL2 (.)
--- 
+--
 -- The problem with these definitions is that basis elements get converted
 -- to values and then decomposed, followed by recombination of the
 -- results.
 
-liftMS :: (AdditiveGroup a) =>
-          (a -> b)
-       -> (MSum a -> MSum b)
+liftMS :: (a -> b) -> (MSum a -> MSum b)
 -- liftMS _ Nothing = Nothing
 -- liftMS h ma = Just (Sum (h (z ma)))
 
@@ -168,8 +246,7 @@
 
 -- | Apply a linear function to each element of a linear map.
 -- @liftL f l == linear f *.* l@, but works more efficiently.
-liftL :: (Functor f, AdditiveGroup (f a)) =>
-         (a -> b) -> MSum (f a) -> MSum (f b)
+liftL :: Functor f => (a -> b) -> MSum (f a) -> MSum (f b)
 liftL = liftMS . fmap
 
 -- | Apply a linear binary function (not to be confused with a bilinear
@@ -186,3 +263,42 @@
           (a -> b -> c -> d)
        -> (MSum (f a) -> MSum (f b) -> MSum (f c) -> MSum (f d))
 liftL3 = liftMS3 . liftA3
+
+{-
+
+
+infixr 9 *.*
+-- | Compose linear maps
+(*.*) :: ( HasBasis u, HasTrie (Basis u)
+         , HasBasis v, HasTrie (Basis v)
+         , VectorSpace w
+         , Scalar v ~ Scalar w ) =>
+         (v :-* w) -> (u :-* v) -> (u :-* w)
+
+-- Simple definition, but only optimizes out uv == zero
+--
+-- (*.*) vw = (fmap.fmap) (lapply vw)
+
+-- Instead, use Nothing/zero if /either/ map is zeroV (exploiting linearity
+-- when uv == zeroV.)
+
+-- Nothing       *.* _             = Nothing
+-- _             *.* Nothing       = Nothing
+-- Just (Sum vw) *.* Just (Sum uv) = Just (Sum (lapply' vw <$> uv))
+
+-- (*.*) = liftA2 (\ (Sum vw) (Sum uv) -> Sum (lapply' vw <$> uv))
+
+-- (*.*) = (liftA2.inSum2) (\ vw uv -> lapply' vw <$> uv)
+(*.*) = (liftA2.inSum2) (\ vw uv -> lapply' vw <$> uv)
+
+-- (*.*) = (liftA2.inSum2) (\ vw -> fmap (lapply' vw))
+
+-- (*.*) = (liftA2.inSum2) (fmap . lapply')
+
+
+-}
+
+-----
+
+(~>) :: (a' -> a) -> (b -> b') -> ((a -> b) -> (a' -> b'))
+(f ~> h) g = h . g . f
diff --git a/src/Data/Maclaurin.hs b/src/Data/Maclaurin.hs
--- a/src/Data/Maclaurin.hs
+++ b/src/Data/Maclaurin.hs
@@ -1,8 +1,7 @@
 {-# LANGUAGE TypeOperators, MultiParamTypeClasses, UndecidableInstances
            , TypeSynonymInstances, FlexibleInstances
            , FlexibleContexts, TypeFamilies
-           , ScopedTypeVariables
-  #-}
+           , ScopedTypeVariables, CPP  #-}
 
 -- The ScopedTypeVariables is there just as a bug work-around.  Without it
 -- I get a bogus error about context mismatch for mutually recursive
@@ -25,7 +24,7 @@
 -- Stability   :  experimental
 -- 
 -- Infinite derivative towers via linear maps, using the Maclaurin
--- representation.  See blog posts <http://conal.net/blog/tag/derivatives/>.
+-- representation.  See blog posts <http://conal.net/blog/tag/derivative/>.
 ----------------------------------------------------------------------
 
 module Data.Maclaurin
@@ -53,6 +52,10 @@
 
 import Data.Boolean
 
+#if MIN_VERSION_base(4,8,0)
+import Prelude hiding ((<*))
+#endif
+
 infixr 9 `D`
 -- | Tower of derivatives.
 data a :> b = D { powVal :: b, derivative :: a :-* (a :> b) }
@@ -71,11 +74,10 @@
 
 infixl 4 <$>>
 -- | Map a /linear/ function over a derivative tower.
-fmapD, (<$>>) :: (HasBasis a, HasTrie (Basis a), AdditiveGroup b) =>
-                 (b -> c) -> (a :> b) -> (a :> c)
+fmapD, (<$>>) :: HasTrie (Basis a) => (b -> c) -> (a :> b) -> (a :> c)
 fmapD f = lf
  where
-   lf (D b0 b') = D (f b0) (liftL lf b')
+   lf (D b0 b') = D (f b0) ((inLMap.liftL) lf b')
 
 (<$>>) = fmapD
 
@@ -84,7 +86,7 @@
           (b -> c -> d) -> (a :> b) -> (a :> c) -> (a :> d)
 liftD2 f = lf
  where
-   lf (D b0 b') (D c0 c') = D (f b0 c0) (liftL2 lf b' c')
+   lf (D b0 b') (D c0 c') = D (f b0 c0) ((inLMap2.liftL2) lf b' c')
 
 
 -- | Apply a /linear/ ternary function over derivative towers.
@@ -95,7 +97,7 @@
 liftD3 f = lf
  where
    lf (D b0 b') (D c0 c') (D d0 d') =
-     D (f b0 c0 d0) (liftL3 lf b' c' d')
+     D (f b0 c0 d0) ((inLMap3.liftL3) lf b' c' d')
 
 
 -- TODO: Can liftD2 and liftD3 be defined in terms of a (<*>>) similar to
@@ -108,9 +110,7 @@
 
 -- | Differentiable identity function.  Sometimes called "the
 -- derivation variable" or similar, but it's not really a variable.
-idD :: ( VectorSpace u, s ~ Scalar u
-       , VectorSpace (u :> u), VectorSpace s
-       , HasBasis u, HasTrie (Basis u)) =>
+idD :: (VectorSpace u , HasBasis u, HasTrie (Basis u)) =>
        u :~> u
 idD = linearD id
 
@@ -165,16 +165,14 @@
 -- | Derivative tower for applying a binary function that distributes over
 -- addition, such as multiplication.  A bit weaker assumption than
 -- bilinearity.  Is bilinearity necessary for correctness here?
-distrib :: forall a b c u.
-           ( HasBasis a, HasTrie (Basis a)
-           , AdditiveGroup b, AdditiveGroup c, AdditiveGroup u) =>
+distrib :: forall a b c u. (HasBasis a, HasTrie (Basis a) , AdditiveGroup u) =>
            (b -> c -> u) -> (a :> b) -> (a :> c) -> (a :> u)
 
 distrib op = (#)
  where
    u@(D u0 u') # v@(D v0 v') =
-     D (u0 `op` v0) ( liftMS (inTrie ((# v) .)) u' ^+^
-                      liftMS (inTrie ((u #) .)) v' )
+     D (u0 `op` v0) ( (inLMap.liftMS) (inTrie ((# v) .)) u' ^+^
+                      (inLMap.liftMS) (inTrie ((u #) .)) v' )
 
 
 -- TODO: I think this distrib is exponential in increasing degree.  Switch
@@ -187,18 +185,19 @@
 instance Show b => Show (a :> b) where
   show (D b0 _) = "D " ++ show b0  ++ " ..."
 
-instance Eq   b => Eq   (a :> b) where (==)    = noOv "(==)"
+instance Eq   (a :> b) where (==)    = noOv "(==)"
 
-instance (AdditiveGroup v, HasBasis u, HasTrie (Basis u), IfB b v) =>
-      IfB b (u :> v) where
+type instance BooleanOf (a :> b) = BooleanOf b
+
+instance (AdditiveGroup v, HasBasis u, HasTrie (Basis u), IfB v) =>
+      IfB (u :> v) where
   ifB = liftD2 . ifB
 
-instance (AdditiveGroup v, HasBasis u, HasTrie (Basis u), OrdB b v) =>
-         OrdB b (u :> v) where
+instance OrdB v => OrdB (u :> v) where
   (<*) = (<*) `on` powVal
 
 instance ( AdditiveGroup b, HasBasis a, HasTrie (Basis a)
-         , OrdB bool b, IfB bool b, Ord  b) => Ord  (a :> b) where
+         , OrdB b, IfB b, Ord  b) => Ord  (a :> b) where
   compare = compare `on` powVal
   min     = minB
   max     = maxB
@@ -213,8 +212,7 @@
   -- Less efficient: adds zero
   -- (^+^)   = liftD2 (^+^)
 
-instance ( HasBasis a, HasTrie (Basis a)
-         , VectorSpace u, AdditiveGroup (Scalar u) )
+instance (HasBasis a, HasTrie (Basis a), VectorSpace u)
       => VectorSpace (a :> u) where
   type Scalar (a :> u) = (a :> Scalar u)
   (*^) = distrib (*^)                     
@@ -236,11 +234,10 @@
 infix  0 >-<
 
 -- | Specialized chain rule.  See also '(\@.)'
-(>-<) :: ( HasBasis a, HasTrie (Basis a), VectorSpace u
-         , AdditiveGroup (Scalar u)) =>
+(>-<) :: (HasBasis a, HasTrie (Basis a), VectorSpace u) =>
          (u -> u) -> ((a :> u) -> (a :> Scalar u))
       -> (a :> u) -> (a :> u)
-f >-< f' = \ u@(D u0 u') -> D (f u0) (liftMS (f' u *^) u')
+f >-< f' = \ u@(D u0 u') -> D (f u0) ((inLMap.liftMS) (f' u *^) u')
 
 
 -- TODO: express '(>-<)' in terms of '(@.)'.  If I can't, then understand why not.
@@ -293,31 +290,21 @@
 
 ---- Misc
 
-pairD :: ( HasBasis a, HasTrie (Basis a)
-         , VectorSpace b, VectorSpace c
-         , Scalar b ~ Scalar c
-         ) => (a:>b,a:>c) -> a:>(b,c)
+pairD :: (HasBasis a, HasTrie (Basis a), VectorSpace b, VectorSpace c)
+      => (a:>b,a:>c) -> a:>(b,c)
 
 pairD (u,v) = liftD2 (,) u v
 
-unpairD :: ( HasBasis a, HasTrie (Basis a)
-           , VectorSpace a, VectorSpace b, VectorSpace c
-           , Scalar b ~ Scalar c
-           ) => (a :> (b,c)) -> (a:>b, a:>c)
+unpairD :: HasTrie (Basis a) => (a :> (b,c)) -> (a:>b, a:>c)
 unpairD d = (fst <$>> d, snd <$>> d)
 
 
 tripleD :: ( HasBasis a, HasTrie (Basis a)
            , VectorSpace b, VectorSpace c, VectorSpace d
-           , Scalar b ~ Scalar c, Scalar c ~ Scalar d
            ) => (a:>b,a:>c,a:>d) -> a:>(b,c,d)
 tripleD (u,v,w) = liftD3 (,,) u v w
 
-untripleD :: ( HasBasis a, HasTrie (Basis a)
-             , VectorSpace a, VectorSpace b, VectorSpace c, VectorSpace d
-             , Scalar b ~ Scalar c, Scalar c ~ Scalar d
-             ) =>
-             (a :> (b,c,d)) -> (a:>b, a:>c, a:>d)
+untripleD :: HasTrie (Basis a) => (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/NumInstances.hs b/src/Data/NumInstances.hs
deleted file mode 100644
--- a/src/Data/NumInstances.hs
+++ /dev/null
@@ -1,174 +0,0 @@
-{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-}
-{-# OPTIONS_GHC -fno-warn-orphans #-}
-----------------------------------------------------------------------
--- |
--- Module      :  Data.NumInstances
--- Copyright   :  (c) Conal Elliott 2008
--- License     :  BSD3
--- 
--- Maintainer  :  conal@conal.net
--- Stability   :  experimental
--- 
--- Number class instances for functions and tuples
-----------------------------------------------------------------------
-
-module Data.NumInstances () where
-
-import Control.Applicative
-
-noOv :: String -> String -> a
-noOv ty meth = error $ meth ++ ": No overloading for " ++ ty
-
-noFun :: String -> a
-noFun = noOv "function"
-
--- Eq & Show are prerequisites for Num, so they need to be faked here
-instance Eq (a->b) where
-  (==) = noFun "(==)"
-  (/=) = noFun "(/=)"
-
-instance Ord b => Ord (a->b) where
-  min = liftA2 min
-  max = liftA2 max
-
-instance Show (a->b) where
-  show      = noFun "show"
-  showsPrec = noFun "showsPrec"
-  showList  = noFun "showList"
-
-instance Num b => Num (a->b) where
-  negate      = fmap negate
-  (+)         = liftA2 (+)
-  (*)         = liftA2 (*)
-  fromInteger = pure . fromInteger
-  abs         = fmap abs
-  signum      = fmap signum
-
-instance Fractional b => Fractional (a->b) where
-  recip        = fmap recip
-  fromRational = pure . fromRational
-
-instance Floating b => Floating (a->b) where
-  pi    = pure pi
-  sqrt  = fmap sqrt
-  exp   = fmap exp
-  log   = fmap log
-  sin   = fmap sin
-  cos   = fmap cos
-  asin  = fmap asin
-  atan  = fmap atan
-  acos  = fmap acos
-  sinh  = fmap sinh
-  cosh  = fmap cosh
-  asinh = fmap asinh
-  atanh = fmap atanh
-  acosh = fmap acosh
-
-
------ Tuples
-
-lift2 :: (a->u) -> (b->v) -> (a,b) -> (u,v)
-lift2 f g (a,b) = (f a, g b)
-
--- Equivalently, lift2 = (***)
-
-instance (Num a, Num b) => Num (a,b) where
-  fromInteger n   = (fromInteger n, fromInteger n)
-  (a,b) + (a',b') = (a+a',b+b')
-  (a,b) - (a',b') = (a-a',b-b')
-  (a,b) * (a',b') = (a*a',b*b')
-  negate = lift2 negate negate
-  abs    = lift2 abs abs
-  signum = lift2 signum signum
-
-instance (Fractional a, Fractional b) => Fractional (a,b) where
-  fromRational x = (fromRational x, fromRational x)
-  recip = lift2 recip recip
-
-instance (Floating a, Floating b) => Floating (a,b) where
-  pi    = (pi,pi)
-  exp   = lift2 exp exp
-  log   = lift2 log log
-  sqrt  = lift2 sqrt sqrt
-  sin   = lift2 sin sin
-  cos   = lift2 cos cos
-  sinh  = lift2 sinh sinh
-  cosh  = lift2 cosh cosh
-  asin  = lift2 asin asin
-  acos  = lift2 acos acos
-  atan  = lift2 atan atan
-  asinh = lift2 asinh asinh
-  acosh = lift2 acosh acosh
-  atanh = lift2 atanh atanh
-
-instance (Num a, Num b, Num c) => Num (a,b,c) where
-  fromInteger n = (fromInteger n, fromInteger n, fromInteger n)
-  (a,b,c) + (a',b',c') = (a+a',b+b',c+c')
-  (a,b,c) - (a',b',c') = (a-a',b-b',c-c')
-  (a,b,c) * (a',b',c') = (a*a',b*b',c*c')
-  negate = lift3 negate negate negate
-  abs    = lift3 abs abs abs
-  signum = lift3 signum signum signum
-
-instance (Fractional a, Fractional b, Fractional c)
-    => Fractional (a,b,c) where
-  fromRational x = (fromRational x, fromRational x, fromRational x)
-  recip = lift3 recip recip recip
-
-
-lift3 :: (a->u) -> (b->v) -> (c->w) -> (a,b,c) -> (u,v,w)
-lift3 f g h (a,b,c) = (f a, g b, h c)
-
-instance (Floating a, Floating b, Floating c)
-    => Floating (a,b,c) where
-  pi    = (pi,pi,pi)
-  exp   = lift3 exp exp exp
-  log   = lift3 log log log
-  sqrt  = lift3 sqrt sqrt sqrt
-  sin   = lift3 sin sin sin
-  cos   = lift3 cos cos cos
-  sinh  = lift3 sinh sinh sinh
-  cosh  = lift3 cosh cosh cosh
-  asin  = lift3 asin asin asin
-  acos  = lift3 acos acos acos
-  atan  = lift3 atan atan atan
-  asinh = lift3 asinh asinh asinh
-  acosh = lift3 acosh acosh acosh
-  atanh = lift3 atanh atanh atanh
-
-
-
-lift4 :: (a->u) -> (b->v) -> (c->w) -> (d->x)
-      -> (a,b,c,d) -> (u,v,w,x)
-lift4 f g h k (a,b,c,d) = (f a, g b, h c, k d)
-
-instance (Num a, Num b, Num c, Num d) => Num (a,b,c,d) where
-  fromInteger n = (fromInteger n, fromInteger n, fromInteger n, fromInteger n)
-  (a,b,c,d) + (a',b',c',d') = (a+a',b+b',c+c',d+d')
-  (a,b,c,d) - (a',b',c',d') = (a-a',b-b',c-c',d-d')
-  (a,b,c,d) * (a',b',c',d') = (a*a',b*b',c*c',d*d')
-  negate = lift4 negate negate negate negate
-  abs    = lift4 abs abs abs abs
-  signum = lift4 signum signum signum signum
-
-instance (Fractional a, Fractional b, Fractional c, Fractional d)
-    => Fractional (a,b,c,d) where
-  fromRational x = (fromRational x, fromRational x, fromRational x, fromRational x)
-  recip = lift4 recip recip recip recip
-
-instance (Floating a, Floating b, Floating c, Floating d)
-    => Floating (a,b,c,d) where
-  pi    = (pi,pi,pi,pi)
-  exp   = lift4 exp exp exp exp
-  log   = lift4 log log log log
-  sqrt  = lift4 sqrt sqrt sqrt sqrt
-  sin   = lift4 sin sin sin sin
-  cos   = lift4 cos cos cos cos
-  sinh  = lift4 sinh sinh sinh sinh
-  cosh  = lift4 cosh cosh cosh cosh
-  asin  = lift4 asin asin asin asin
-  acos  = lift4 acos acos acos acos
-  atan  = lift4 atan atan atan atan
-  asinh = lift4 asinh asinh asinh asinh
-  acosh = lift4 acosh acosh acosh acosh
-  atanh = lift4 atanh atanh atanh atanh
diff --git a/src/Data/VectorSpace.hs b/src/Data/VectorSpace.hs
--- a/src/Data/VectorSpace.hs
+++ b/src/Data/VectorSpace.hs
@@ -1,18 +1,21 @@
 {-# LANGUAGE MultiParamTypeClasses, TypeOperators
-           , TypeFamilies, UndecidableInstances
- #-}
+           , TypeFamilies, UndecidableInstances, CPP
+           , FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances  #-}
+{-# LANGUAGE DefaultSignatures   #-}
+{-# LANGUAGE ScopedTypeVariables #-}
 {-# OPTIONS_GHC -Wall #-}
 ----------------------------------------------------------------------
 -- |
 -- Module      :   Data.VectorSpace
 -- Copyright   :  (c) Conal Elliott and Andy J Gill 2008
 -- License     :  BSD3
--- 
+--
 -- Maintainer  :  conal@conal.net, andygill@ku.edu
 -- Stability   :  experimental
--- 
+--
 -- Vector spaces
--- 
+--
 -- This version uses associated types instead of fundeps and
 -- requires ghc-6.10 or later
 ----------------------------------------------------------------------
@@ -20,83 +23,130 @@
 -- NB: I'm attempting to replace fundeps with associated types.  See
 -- NewVectorSpace.hs.  Ran into trouble with type equality constraints not
 -- getting propagated.  Manuel Ch is looking into it.
--- 
+--
 -- Blocking bug: http://hackage.haskell.org/trac/ghc/ticket/2448
 
 module Data.VectorSpace
   ( module Data.AdditiveGroup
   , VectorSpace(..), (^/), (^*)
   , InnerSpace(..)
-  , lerp, magnitudeSq, magnitude, normalized
+  , lerp, linearCombo, magnitudeSq, magnitude, normalized, project
   ) where
-
+#if !(MIN_VERSION_base(4,8,0))
 import Control.Applicative (liftA2)
+#endif
 import Data.Complex hiding (magnitude)
+import Foreign.C.Types (CSChar, CInt, CShort, CLong, CLLong, CIntMax, CFloat, CDouble)
+import Data.Ratio
+import Data.Kind
 
 import Data.AdditiveGroup
 import Data.MemoTrie
 
+import Data.VectorSpace.Generic
+import qualified GHC.Generics as Gnrx
+import GHC.Generics (Generic, (:*:)(..))
+
 infixr 7 *^
 
 -- | Vector space @v@.
 class AdditiveGroup v => VectorSpace v where
-  type Scalar v :: *
+  type Scalar v :: Type
+  type Scalar v = Scalar (VRep v)
   -- | Scale a vector
   (*^) :: Scalar v -> v -> v
+  default (*^) :: (Generic v, VectorSpace (VRep v), Scalar (VRep v) ~ Scalar v)
+                    => Scalar v -> v -> v
+  μ *^ v = Gnrx.to (μ *^ Gnrx.from v :: VRep v)
+  {-# INLINE (*^) #-}
 
 infixr 7 <.>
 
 -- | Adds inner (dot) products.
-class VectorSpace v => InnerSpace v where
+class (VectorSpace v, AdditiveGroup (Scalar v)) => InnerSpace v where
   -- | Inner/dot product
   (<.>) :: v -> v -> Scalar v
+  default (<.>) :: (Generic v, InnerSpace (VRep v), Scalar (VRep v) ~ Scalar v)
+                    => v -> v -> Scalar v
+  v<.>w = (Gnrx.from v :: VRep v) <.> Gnrx.from w
+  {-# INLINE (<.>) #-}
 
 infixr 7 ^/
 infixl 7 ^*
 
 -- | Vector divided by scalar
 (^/) :: (VectorSpace v, s ~ Scalar v, Fractional s) => v -> s -> v
-v ^/ s = (1/s) *^ v
+v ^/ s = recip s *^ v
+{-# INLINE (^/) #-}
 
 -- | Vector multiplied by scalar
 (^*) :: (VectorSpace v, s ~ Scalar v) => v -> s -> v
 (^*) = flip (*^)
+{-# INLINE (^*) #-}
 
 -- | Linear interpolation between @a@ (when @t==0@) and @b@ (when @t==1@).
 
 -- lerp :: (VectorSpace v, s ~ Scalar v, Num s) => v -> v -> s -> v
 lerp :: VectorSpace v => v -> v -> Scalar v -> v
 lerp a b t = a ^+^ t *^ (b ^-^ a)
+{-# INLINE lerp #-}
 
+-- | Linear combination of vectors
+linearCombo :: VectorSpace v => [(v,Scalar v)] -> v
+linearCombo ps = sumV [v ^* s | (v,s) <- ps]
+{-# INLINE linearCombo #-}
+
 -- | Square of the length of a vector.  Sometimes useful for efficiency.
 -- See also 'magnitude'.
 magnitudeSq :: (InnerSpace v, s ~ Scalar v) => v -> s
 magnitudeSq v = v <.> v
+{-# INLINE magnitudeSq #-}
 
 -- | Length of a vector.   See also 'magnitudeSq'.
 magnitude :: (InnerSpace v, s ~ Scalar v, Floating s) =>  v -> s
 magnitude = sqrt . magnitudeSq
+{-# INLINE magnitude #-}
 
--- | Vector in same direction as given one but with length of one.  If
--- given the zero vector, then return it.
+-- | Vector in same direction as given one but with length of one.
+-- Divides by zero for the zero vector.
 normalized :: (InnerSpace v, s ~ Scalar v, Floating s) =>  v -> v
 normalized v = v ^/ magnitude v
+{-# INLINE normalized #-}
 
-instance VectorSpace Double where
-  type Scalar Double = Double
-  (*^) = (*)
-instance InnerSpace  Double where (<.>) = (*)
+-- | @project u v@ computes the projection of @v@ onto @u@.
+project :: (InnerSpace v, s ~ Scalar v, Fractional s) => v -> v -> v
+project u v = ((v <.> u) / magnitudeSq u) *^ u
+{-# INLINE project #-}
 
-instance VectorSpace Float  where
-  type Scalar Float = Float
-  (*^)  = (*)
-instance InnerSpace  Float  where (<.>) = (*)
+#define ScalarType(t) \
+  instance VectorSpace (t) where \
+    { type Scalar (t) = (t) \
+    ; (*^) = (*) } ; \
+  instance InnerSpace  (t) where (<.>) = (*)
 
+ScalarType(Int)
+ScalarType(Integer)
+ScalarType(Double)
+ScalarType(Float)
+ScalarType(CSChar)
+ScalarType(CInt)
+ScalarType(CShort)
+ScalarType(CLong)
+ScalarType(CLLong)
+ScalarType(CIntMax)
+ScalarType(CDouble)
+ScalarType(CFloat)
+
+instance Integral a => VectorSpace (Ratio a) where
+  type Scalar (Ratio a) = Ratio a
+  (*^) = (*)
+instance Integral a => InnerSpace  (Ratio a) where (<.>) = (*)
+
 instance (RealFloat v, VectorSpace v) => VectorSpace (Complex v) where
   type Scalar (Complex v) = Scalar v
   s*^(u :+ v) = s*^u :+ s*^v
 
-instance (RealFloat v, InnerSpace v, s ~ Scalar v, AdditiveGroup s)
+instance (RealFloat v, InnerSpace v)
      => InnerSpace (Complex v) where
   (u :+ v) <.> (u' :+ v') = (u <.> u') ^+^ (v <.> v')
 
@@ -116,8 +166,7 @@
   s *^ (u,v) = (s*^u,s*^v)
 
 instance ( InnerSpace u, s ~ Scalar u
-         , InnerSpace v, s ~ Scalar v
-         , AdditiveGroup (Scalar v) )
+         , InnerSpace v, s ~ Scalar v )
     => InnerSpace (u,v) where
   (u,v) <.> (u',v') = (u <.> u') ^+^ (v <.> v')
 
@@ -130,8 +179,7 @@
 
 instance ( InnerSpace u, s ~ Scalar u
          , InnerSpace v, s ~ Scalar v
-         , InnerSpace w, s ~ Scalar w
-         , AdditiveGroup s )
+         , InnerSpace w, s ~ Scalar w )
     => InnerSpace (u,v,w) where
   (u,v,w) <.> (u',v',w') = u<.>u' ^+^ v<.>v' ^+^ w<.>w'
 
@@ -146,8 +194,7 @@
 instance ( InnerSpace u, s ~ Scalar u
          , InnerSpace v, s ~ Scalar v
          , InnerSpace w, s ~ Scalar w
-         , InnerSpace x, s ~ Scalar x
-         , AdditiveGroup s )
+         , InnerSpace x, s ~ Scalar x )
     => InnerSpace (u,v,w,x) where
   (u,v,w,x) <.> (u',v',w',x') = u<.>u' ^+^ v<.>v' ^+^ w<.>w' ^+^ x<.>x'
 
@@ -185,7 +232,7 @@
 
 
 
-instance (InnerSpace a, AdditiveGroup (Scalar a)) => InnerSpace (Maybe a) where
+instance InnerSpace a => InnerSpace (Maybe a) where
   -- dotting with zero (vector) yields zero (scalar)
   Nothing <.> _     = zeroV
   _ <.> Nothing     = zeroV
@@ -194,3 +241,30 @@
 --   mu <.> mv = fromMaybe zeroV (liftA2 (<.>) mu mv)
 
 --   (<.>) = (fmap.fmap) (fromMaybe zeroV) (liftA2 (<.>))
+
+
+instance VectorSpace a => VectorSpace (Gnrx.Rec0 a s) where
+  type Scalar (Gnrx.Rec0 a s) = Scalar a
+  μ *^ Gnrx.K1 v = Gnrx.K1 $ μ*^v
+  {-# INLINE (*^) #-}
+instance VectorSpace (f p) => VectorSpace (Gnrx.M1 i c f p) where
+  type Scalar (Gnrx.M1 i c f p) = Scalar (f p)
+  μ *^ Gnrx.M1 v = Gnrx.M1 $ μ*^v
+  {-# INLINE (*^) #-}
+instance (VectorSpace (f p), VectorSpace (g p), Scalar (f p) ~ Scalar (g p))
+         => VectorSpace ((f :*: g) p) where
+  type Scalar ((f:*:g) p) = Scalar (f p)
+  μ *^ (x:*:y) = μ*^x :*: μ*^y
+  {-# INLINE (*^) #-}
+
+instance InnerSpace a => InnerSpace (Gnrx.Rec0 a s) where
+  Gnrx.K1 v <.> Gnrx.K1 w = v<.>w
+  {-# INLINE (<.>) #-}
+instance InnerSpace (f p) => InnerSpace (Gnrx.M1 i c f p) where
+  Gnrx.M1 v <.> Gnrx.M1 w = v<.>w
+  {-# INLINE (<.>) #-}
+instance ( InnerSpace (f p), InnerSpace (g p)
+         , Scalar (f p) ~ Scalar (g p), Num (Scalar (f p)) )
+         => InnerSpace ((f :*: g) p) where
+  (x:*:y) <.> (ξ:*:υ) = x<.>ξ + y<.>υ
+  {-# INLINE (<.>) #-}
diff --git a/src/Data/VectorSpace/Generic.hs b/src/Data/VectorSpace/Generic.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/VectorSpace/Generic.hs
@@ -0,0 +1,20 @@
+-- |
+-- Module      :   Data.VectorSpace.Generic
+-- Copyright   :  (c) Conal Elliott and Justus Sagemüller 2017
+-- License     :  BSD3
+-- 
+-- Maintainer  :  conal@conal.net, (@) jsagemue $ uni-koeln.de
+-- Stability   :  experimental
+-- 
+-- Underpinnings of the type that represents vector / affine / etc. spaces
+-- with GHC generics
+
+module Data.VectorSpace.Generic where
+
+
+import qualified GHC.Generics as Gnrx
+
+import Data.Void
+
+
+type VRep v = Gnrx.Rep v Void
diff --git a/tests/src/Perf.hs b/tests/src/Perf.hs
deleted file mode 100644
--- a/tests/src/Perf.hs
+++ /dev/null
@@ -1,203 +0,0 @@
-{-# LANGUAGE TypeOperators, MultiParamTypeClasses, UndecidableInstances, FlexibleInstances
-           , TypeFamilies, FlexibleContexts
-  #-}
-
-
--- 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.Basis
-import Data.MemoTrie
-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 ()
-
-	let samples = samples_2d
-
-	sequence_ [ do t <- getClockTime
-		       loop msg fun t 0 samples
-		  | (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 :: (HasBasis s s, Basis s ~ (), Floating s, InnerSpace s s)
-         => SurfPt s -> ((s,s,s),(s,s,s))
-toVN3 v = ( powVal v
-	  , powVal (normal v)
-	  )
-vector3D :: (HasBasis a s, HasTrie (Basis a), VectorSpace s s) => (a :> s,a :> s,a :> s) -> (a :> (s,s,s))
-vector3D (u,v,w) = liftD3 (,,) u v w
-unvector2D :: (HasBasis a s, HasTrie (Basis a), 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,7 +1,7 @@
 Name:                vector-space
-Version:             0.5.9
-Cabal-Version:       >= 1.2
-Synopsis:            Vector & affine spaces, linear maps, and derivatives (requires ghc 6.9 or better)
+Version:             0.19
+Cabal-Version:       >= 1.10
+Synopsis:            Vector & affine spaces, linear maps, and derivatives
 Category:            math
 Description:
   /vector-space/ provides classes and generic operations for vector
@@ -11,24 +11,30 @@
   (scalars, vectors, matrices, ...).
   .
   /Warning/: this package depends on type families working fairly well,
-  and requires ghc version at least 6.9.
+  requiring GHC version at least 6.9.
   .
   Project wiki page: <http://haskell.org/haskellwiki/vector-space>
   .
-  &#169; 2008 by Conal Elliott; BSD3 license.
+  &#169; 2008-2012 by Conal Elliott; BSD3 license.
 Author:              Conal Elliott 
 Maintainer:          conal@conal.net
-Homepage:            http://haskell.org/haskellwiki/vector-space
-Package-Url:         http://code.haskell.org/vector-space
-Copyright:           (c) 2008 by Conal Elliott
+Copyright:           (c) 2008-2012 by Conal Elliott
 License:             BSD3
+License-File:        COPYING
 Stability:           experimental
 build-type:          Simple
 
+source-repository head
+  type:     git
+  location: git://github.com/conal/vector-space.git
+
 Library
+  default-language:  Haskell2010
   hs-Source-Dirs:      src
   Extensions:          
-  Build-Depends:       base, MemoTrie >= 0.4.2, Boolean
+  Build-Depends:       base<5, MemoTrie >= 0.5
+                     , Boolean >= 0.1.0
+                     , NumInstances >= 1.0
   Exposed-Modules:     
                      Data.AdditiveGroup
                      Data.VectorSpace
@@ -39,14 +45,11 @@
                      Data.Derivative
                      Data.Cross
                      Data.AffineSpace
-                     Data.NumInstances
-
+  Other-Modules:     
+                     Data.VectorSpace.Generic
 
   -- This library relies on type families working as well as in 6.10.
-  if impl(ghc < 6.10) {
-    buildable: False
-  }
-  ghc-options:         -Wall -O2
-  ghc-prof-options:    -prof -auto-all 
-
--- For ghc-options: -ddump-simpl-stats -ddump-rules -ddump-simpl -ddump-simpl-phases
+  if  impl(ghc < 6.10) { buildable: False }
+  if !impl(ghc >= 7.6) { Build-Depends: ghc-prim >= 0.2 }
+  if !impl(ghc >= 7.9) { Build-Depends: void >= 0.4 }
+  if !impl(ghc >= 8.0) { Build-Depends: semigroups >= 0.16 }