diff --git a/clifford.cabal b/clifford.cabal
--- a/clifford.cabal
+++ b/clifford.cabal
@@ -10,7 +10,7 @@
 -- PVP summary:      +-+------- breaking API changes
 --                   | | +----- non-breaking API additions
 --                   | | | +--- code changes with no API change
-version:             0.1.0.1
+version:             0.1.0.3
 
 -- A short (one-line) description of the package.
 synopsis:            A Clifford algebra library
@@ -48,6 +48,7 @@
 -- Constraint on the version of Cabal needed to build this package.
 cabal-version:       >=1.10
 
+
 library
   -- Modules exported by the library.
   exposed-modules:     Numeric.Clifford.Blade, Numeric.Clifford.Multivector, Numeric.Clifford.NumericIntegration, Numeric.Clifford.NumericIntegration.DefaultIntegrators, Numeric.Clifford.ClassicalMechanics
@@ -56,8 +57,7 @@
   -- other-modules:       
   
   -- LANGUAGE extensions used by modules in this package.
-  -- other-extensions:
-  
+  other-extensions:    NoImplicitPrelude, ScopedTypeVariables, GADTs, DataKinds,KindSignatures,UnicodeSyntax,FlexibleContexts, RankNTypes,TemplateHaskell,NoMonomorphismRestriction,MultiParamTypeClasses,FlexibleInstances, TypeOperators
   -- Other library packages from which modules are imported.
   build-depends:       base >=4.6 && <4.9, numeric-prelude >= 0.4.0.1 && < 0.5.0, permutation >= 0.4.1 && < 0.5, 
                        data-ordlist >= 0.4.5 && < 0.5,  converge >= 0.1.0.1 && < 0.2, lens >= 4.0.3 && < 4.1, 
diff --git a/src/Numeric/Clifford/Blade.lhs b/src/Numeric/Clifford/Blade.lhs
--- a/src/Numeric/Clifford/Blade.lhs
+++ b/src/Numeric/Clifford/Blade.lhs
@@ -49,8 +49,8 @@
 import Test.QuickCheck
 import Control.Lens hiding (indices)
 import Data.DeriveTH
-import GHC.TypeLits hiding (isEven)
-import GHC.Real (fromIntegral)
+import GHC.TypeLits hiding (isEven, isOdd)
+import GHC.Real (fromIntegral, toInteger)
 import Algebra.Field
 import Debug.Trace
 --trace _ a = a
@@ -63,19 +63,19 @@
 \texttt{bScale} is the amplitude of the blade. \texttt{bIndices} are the indices for the basis. 
 \begin{code}
 
-data Blade (n :: Nat) f where
-    Blade :: forall n f . (SingI n, Algebra.Field.C f) => {_scale :: f, _indices :: [Natural]} -> Blade n f
+data Blade (p :: Nat) (q :: Nat) f where
+    Blade :: forall p q f . (SingI p, SingI q, Algebra.Field.C f) => {_scale :: f, _indices :: [Natural]} -> Blade p q f
 
 -- makeLenses ''Blade
-scale :: Lens' (Blade n f) f
+scale :: Lens' (Blade p q f) f
 scale = lens _scale (\blade v -> blade {_scale = v})
-indices :: Lens' (Blade n f) [Natural]
+indices :: Lens' (Blade p q f) [Natural]
 indices = lens _indices (\blade v -> blade {_indices = v})
-dimension :: forall (n::Nat) f. SingI n => Blade n f ->  Natural
-dimension _ = toNatural  ((GHC.Real.fromIntegral $ fromSing (sing :: Sing n))::Word)
+dimension :: forall (p::Nat) (q::Nat) f. (SingI p, SingI q) => Blade p q f ->  (Natural,Natural)
+dimension _ = (toNatural  ((GHC.Real.fromIntegral $ fromSing (sing :: Sing p))::Word),toNatural((GHC.Real.fromIntegral $ fromSing (sing :: Sing q))::Word))
 bScale b =  b^.scale
 bIndices b = b^.indices
-instance(Show f) =>  Show (Blade n f) where
+instance(Show f) =>  Show (Blade p q f) where
     --TODO: Do this with HaTeX
     show  (Blade scale indices) = pref ++  if null indices then "" else basis where
                         pref = show scale
@@ -84,7 +84,7 @@
                         vecced index = "\\vec{e_{" ++ show index ++ "}}"
                                                 
                         
-instance (Algebra.Additive.C f, Eq f) => Eq (Blade n f) where
+instance (Algebra.Additive.C f, Eq f) => Eq (Blade p q f) where
    a == b = aScale == bScale && aIndices == bIndices where
                  (Blade aScale aIndices) = bladeNormalForm a
                  (Blade bScale bIndices) = bladeNormalForm b
@@ -93,10 +93,10 @@
 
 For example, a scalar could be constructed like so: \texttt{Blade s []}
 \begin{code}
-scalarBlade :: (Algebra.Field.C f, SingI n) => f -> Blade n f
+scalarBlade :: (Algebra.Field.C f, SingI p, SingI q) => f -> Blade p q f
 scalarBlade d = Blade d []
 
-zeroBlade :: (Algebra.Field.C f, SingI n) => Blade n f
+zeroBlade :: (Algebra.Field.C f, SingI p, SingI q) => Blade p q f
 zeroBlade = scalarBlade Algebra.Additive.zero
 
 bladeNonZero b = b^.scale /= Algebra.Additive.zero
@@ -117,32 +117,38 @@
 
 
 \begin{code}
-bladeNormalForm :: forall (n::Nat) f. Blade n f -> Blade n f
+bladeNormalForm :: forall (p::Nat) (q::Nat) f. Blade p q f -> Blade p q f
 bladeNormalForm (Blade scale indices)  = result 
         where
-             result = if (any (\i -> (GHC.Real.fromIntegral i) >n') indices) then trace "Blade contains vector with i > d" zeroBlade else Blade scale' uniqueSorted
-             n' = fromSing (sing :: Sing n)
+             result = if (any (\i -> (GHC.Real.toInteger i) > d) indices) then trace "Blade contains vector with i > d" zeroBlade else Blade scale' uniqueSorted
+             p' = (fromSing (sing :: Sing p)) :: Integer
+             q' = (fromSing (sing :: Sing q)) :: Integer
+             d = p' + q' 
              numOfIndices = length indices
              (sorted, perm) = Data.Permute.sort numOfIndices indices
-             scale' = if isEven perm then scale else negate scale
-             uniqueSorted = removeDupPairs sorted
+             scale' = if (isEven perm) /= (negated)  then scale else negate scale
+             (uniqueSorted,negated) = removeDupPairs [] sorted False
                             where
-                              removeDupPairs [] = []
-                              removeDupPairs [x] = [x]
-                              removeDupPairs (x:y:rest) | x == y = removeDupPairs rest
-                                                        | otherwise = x : removeDupPairs (y:rest)
+                              removeDupPairs :: [Natural] -> [Natural] -> Bool -> ([Natural],Bool)
+                              removeDupPairs accum [] negated = (accum,negated)
+                              removeDupPairs accum [x] negated = (accum++[x],negated)
+                              removeDupPairs accum (x:y:rest) negated  | x == y  = 
+                                                                            if  GHC.Real.toInteger x >  p' 
+                                                                            then removeDupPairs accum rest (not negated)
+                                                                            else removeDupPairs accum rest negated
+                                                        | otherwise = removeDupPairs (accum++[x]) (y:rest) negated
 \end{code}
 
 What is the grade of a blade? Just the number of indices.
 
 \begin{code}
-grade :: Blade n f -> Integer
-grade = toInteger . length . bIndices 
+grade :: Blade p q f -> Integer
+grade = GHC.Real.toInteger . length . bIndices 
 
-bladeIsOfGrade :: Blade n f -> Integer -> Bool
+bladeIsOfGrade :: Blade p q f -> Integer -> Bool
 blade `bladeIsOfGrade` k = grade blade == k
 
-bladeGetGrade :: Integer -> Blade n f -> Blade n f
+bladeGetGrade :: Integer -> Blade p q f -> Blade p q f
 bladeGetGrade k blade@(Blade _ _) =
     if blade `bladeIsOfGrade` k then blade else zeroBlade
 \end{code}
@@ -152,9 +158,9 @@
 First up for operations: Blade multiplication. This is no more than assembling orthogonal vectors into k-vectors. 
 
 \begin{code}
-bladeMul ::  Blade n f -> Blade n f-> Blade n f
-bladeMul x@(Blade _ _) y@(Blade _ _)= bladeNormalForm $ Blade (bScale x Algebra.Ring.* bScale y) (bIndices x ++ bIndices y)
-multiplyBladeList :: (SingI n, Algebra.Field.C f) => [Blade n f] -> Blade n f
+bladeMul ::  Blade p q f -> Blade p q f-> Blade p q f
+bladeMul x@(Blade _ _) y@(Blade _ _)= bladeNormalForm $ Blade (bScale x Algebra.Ring.* bScale y) (bIndices x ++ bIndices y) -- HAVE TO MAKE Q INDICES SQUARE NEGATIVE
+multiplyBladeList :: (SingI p, SingI q, Algebra.Field.C f) => [Blade p q f] -> Blade p q f
 multiplyBladeList [] = error "Empty blade list!"
 multiplyBladeList (a:[]) = a
 multiplyBladeList a = bladeNormalForm $ Blade scale indices where
@@ -167,7 +173,7 @@
 Next up: The outer (wedge) product, denoted by $∧$ :3
 
 \begin{code}
-bWedge :: Blade n f -> Blade n f -> Blade n f
+bWedge :: Blade p q f -> Blade p q f -> Blade p q f
 bWedge x y = bladeNormalForm $ bladeGetGrade k xy
              where
                k = grade x + grade y
@@ -179,7 +185,7 @@
 
 
 \begin{code}
-bDot ::Blade n f -> Blade n f -> Blade n f
+bDot ::Blade p q f -> Blade p q f -> Blade p q f
 bDot x y = bladeNormalForm $ bladeGetGrade k xy
           where
             k = Algebra.Absolute.abs $ grade x - grade y
@@ -194,7 +200,7 @@
 Now for linear combinations of (possibly different basis) blades. To start with, let's order basis blades:
 
 \begin{code}
-instance (Algebra.Additive.C f, Ord f) => Ord (Blade n f) where
+instance (Algebra.Additive.C f, Ord f) => Ord (Blade p q f) where
     compare a b | bIndices a == bIndices b = compare (bScale a) (bScale b)
                 | otherwise =  compare (bIndices a) (bIndices b)
 
diff --git a/src/Numeric/Clifford/ClassicalMechanics.lhs b/src/Numeric/Clifford/ClassicalMechanics.lhs
--- a/src/Numeric/Clifford/ClassicalMechanics.lhs
+++ b/src/Numeric/Clifford/ClassicalMechanics.lhs
@@ -16,7 +16,7 @@
 \begin{code}
 {-# LANGUAGE NoImplicitPrelude, FlexibleContexts, RankNTypes, ScopedTypeVariables, DeriveDataTypeable #-}
 {-# LANGUAGE NoMonomorphismRestriction, UnicodeSyntax, GADTs, KindSignatures, DataKinds #-}
-{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleInstances, TypeOperators #-}
 {-# LANGUAGE TemplateHaskell #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 \end{code}
@@ -59,9 +59,9 @@
 import Debug.Trace
 --trace _ a = a
 
-data EnergyMethod (d::Nat) f = Hamiltonian{ _dqs :: [DynamicSystem d f -> Multivector d f], _dps :: [DynamicSystem d f -> Multivector d f]}
+data EnergyMethod (p::Nat) (q::Nat) f = Hamiltonian{ _dqs :: [DynamicSystem p q f -> Multivector p q f], _dps :: [DynamicSystem p q f -> Multivector p q f]}
 
-data DynamicSystem (d::Nat) f = DynamicSystem {_time :: f, coordinates :: [Multivector d f], _momenta :: [Multivector d f], _energyFunction :: EnergyMethod d f, _projector :: DynamicSystem d f -> DynamicSystem d f}
+data DynamicSystem (p::Nat) (q::Nat) f = DynamicSystem {_time :: f, coordinates :: [Multivector p q f], _momenta :: [Multivector p q f], _energyFunction :: EnergyMethod p q f, _projector :: DynamicSystem p q f -> DynamicSystem p q f}
 
 makeLenses ''EnergyMethod
 makeLenses ''DynamicSystem
@@ -80,14 +80,14 @@
 
 Now to make a physical object.
 \begin{code}
-data ReferenceFrame (d::Nat) t = ReferenceFrame {basisVectors :: [Multivector d t]}
-psuedoScalar' :: forall f (d::Nat). (Ord f, Algebra.Field.C f, SingI d) => ReferenceFrame d f -> Multivector d f
+data ReferenceFrame (p::Nat) (q::Nat) t = ReferenceFrame {basisVectors :: [Multivector p q t]}
+psuedoScalar' :: forall f (p::Nat) (q::Nat). (Ord f, Algebra.Field.C f, SingI p, SingI q) => ReferenceFrame p q f -> Multivector p q f
 psuedoScalar'  = multiplyList . basisVectors
-psuedoScalar :: forall (d::Nat) f. (Ord f, Algebra.Field.C f, SingI d) =>  Multivector d f
-psuedoScalar = one `e` [1..(toNatural  ((fromIntegral $ fromSing (sing :: Sing d))::Word))] 
 
+
+
 a `cross` b = (negate $ one)`e`[1,2,3] * (a ∧ b)
-data PhysicalVector (d::Nat) t = PhysicalVector {dimension :: Natural, r :: Multivector d t, referenceFrame :: ReferenceFrame d t}
+data PhysicalVector (p::Nat) (q::Nat) t = PhysicalVector {dimension :: Natural, r :: Multivector p q t, referenceFrame :: ReferenceFrame p q t}
 {-squishToDimension (PhysicalVector d (BladeSum terms) f) = PhysicalVector d r' f where
     r' = BladeSum terms' where
         terms' = terms & filter (\(Blade _ ind) -> all (\k -> k <= d) ind)
@@ -95,14 +95,14 @@
     r' = BladeSum terms' where
         terms' = terms & filter (\(Blade _ ind) -> all (\k -> k <= d) ind)-}
 
-data RigidBody (d::Nat) f where
- RigidBody:: (Algebra.Field.C f, Algebra.Module.C f (Multivector d f)) =>  {position :: PhysicalVector d f,
-                              _momentum :: PhysicalVector d f,
+data RigidBody (p::Nat) (q::Nat) f where
+ RigidBody:: (Algebra.Field.C f, Algebra.Module.C f (Multivector p q f)) =>  {position :: PhysicalVector p q f,
+                              _momentum :: PhysicalVector p q f,
                               _mass :: f,
-                              _attitude :: PhysicalVector d f,
-                              _angularMomentum :: PhysicalVector d f,
-                              _inertia :: PhysicalVector d f
-                             } -> RigidBody d f
+                              _attitude :: PhysicalVector p q f,
+                              _angularMomentum :: PhysicalVector p q f,
+                              _inertia :: PhysicalVector p q f
+                             } -> RigidBody p q f
 
 --makeLenses ''RigidBody doesn't actually work
 {- Things to do: 
diff --git a/src/Numeric/Clifford/Multivector.lhs b/src/Numeric/Clifford/Multivector.lhs
--- a/src/Numeric/Clifford/Multivector.lhs
+++ b/src/Numeric/Clifford/Multivector.lhs
@@ -21,7 +21,7 @@
 {-# LANGUAGE NoImplicitPrelude, FlexibleContexts, RankNTypes, ScopedTypeVariables, DeriveDataTypeable #-}
 {-# LANGUAGE NoMonomorphismRestriction, UnicodeSyntax, GADTs#-}
 {-# LANGUAGE FlexibleInstances, StandaloneDeriving, KindSignatures, DataKinds #-}
-{-# LANGUAGE TemplateHaskell #-}
+{-# LANGUAGE TemplateHaskell, TypeOperators #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 \end{code}
 %if False
@@ -83,17 +83,17 @@
 A multivector is nothing but a linear combination of basis blades.
 
 \begin{code}
-data Multivector (n::Nat) f where
-    BladeSum :: forall n f . (SingI n, Algebra.Field.C f, Ord f) => { _terms :: [Blade n f]} -> Multivector n f
+data Multivector (p::Nat) (q::Nat) f where
+    BladeSum :: forall p q f . (Ord f, Algebra.Field.C f, SingI p, SingI q) => { _terms :: [Blade p q f]} -> Multivector p q f
 
-deriving instance Eq (Multivector n f)
-deriving instance Ord (Multivector n f)
-deriving instance (Show f) => Show (Multivector n f)
+deriving instance Eq (Multivector p q f)
+deriving instance Ord (Multivector p q f)
+deriving instance (Show f) => Show (Multivector p q f)
 
-dimension :: forall (n::Nat) f. SingI n => Multivector n f ->  Natural
-dimension _ = toNatural  ((fromIntegral $ fromSing (sing :: Sing n))::Word)
+dimension :: forall (p::Nat) (q::Nat) f. (SingI p, SingI q) => Multivector p q f ->  (Natural,Natural)
+dimension _ = (toNatural  ((fromIntegral $ fromSing (sing :: Sing p))::Word),toNatural  ((fromIntegral $ fromSing (sing :: Sing q))::Word))
 
-terms :: Lens' (Multivector n f) [Blade n f]
+terms :: Lens' (Multivector p q f) [Blade p q f]
 terms = lens _terms (\bladeSum v -> bladeSum {_terms = v})
 
 mvNormalForm (BladeSum terms) = BladeSum $ if null resultant then [scalarBlade Algebra.Additive.zero] else resultant  where
@@ -102,7 +102,7 @@
 
 addLikeTerms' = sumLikeTerms . groupLikeTerms
 
-groupLikeTerms ::Eq f =>  [Blade n f] -> [[Blade n f]]
+groupLikeTerms ::Eq f =>  [Blade p q f] -> [[Blade p q f]]
 groupLikeTerms = groupBy (\a b -> a^.indices == b^.indices)
 
 compensatedSum' :: (Algebra.Additive.C f) => [f] -> f
@@ -135,7 +135,7 @@
 --      ei = multiplyAdd eim1 ai pii
 
 
-multiplyOutBlades :: (SingI n, Algebra.Ring.C a) => [Blade n a] -> [Blade n a] -> [Blade n a]
+multiplyOutBlades :: (SingI p, SingI q, Algebra.Ring.C a) => [Blade p q a] -> [Blade p q a] -> [Blade p q a]
 multiplyOutBlades x y = [bladeMul l r | l <-x, r <- y]
 
 --multiplyList :: Algebra.Ring.C t => [Multivector t] -> Multivector t
@@ -150,29 +150,29 @@
 
 sumList xs = mvNormalForm $ BladeSum $ concat $ map mvTerms xs
 
-sumLikeTerms :: (Algebra.Field.C f, SingI n) => [[Blade n f]] -> [Blade n f]
+sumLikeTerms :: (Algebra.Field.C f, SingI p, SingI q) => [[Blade p q f]] -> [Blade p q f]
 sumLikeTerms blades = map (\sameIxs -> map bScale sameIxs & compensatedSum' & (\result -> Blade result ((head sameIxs) & bIndices))) blades
 
 
-instance (Algebra.Field.C f, SingI n, Ord f) => Data.Monoid.Monoid (Sum (Multivector n f)) where
+instance (Algebra.Field.C f, SingI p, SingI q, Ord f) => Data.Monoid.Monoid (Sum (Multivector p q f)) where
     mempty = Data.Monoid.Sum Algebra.Additive.zero
     mappend (Data.Monoid.Sum x) (Data.Monoid.Sum y)= Data.Monoid.Sum  (x + y)
     mconcat = Data.Monoid.Sum . sumList . map getSum
 
-instance (Algebra.Field.C f, SingI n, Ord f) => Data.Monoid.Monoid (Product (Multivector n f)) where
+instance (Algebra.Field.C f, SingI p, SingI q, Ord f) => Data.Monoid.Monoid (Product (Multivector p q f)) where
     mempty = Product one
     mappend (Product x) (Product y) = Product (x * y)
     mconcat = Product . foldl (*) one . map getProduct
 
 --Constructs a multivector from a scaled blade.
-e :: (Algebra.Field.C f, Ord f, SingI n) => f -> [Natural] -> Multivector n f
+e :: (Algebra.Field.C f, Ord f, SingI p, SingI q) => f -> [Natural] -> Multivector p q f
 s `e` indices = mvNormalForm $ BladeSum [Blade s indices]
 scalar s = s `e` []
 
 
-instance (Control.DeepSeq.NFData f) => Control.DeepSeq.NFData (Multivector n f)
-instance (Control.DeepSeq.NFData f) => Control.DeepSeq.NFData (Blade n f)
-instance (Algebra.Field.C f, Ord f, SingI n) => Algebra.Additive.C (Multivector n f) where
+instance (Control.DeepSeq.NFData f) => Control.DeepSeq.NFData (Multivector p q f)
+instance (Control.DeepSeq.NFData f) => Control.DeepSeq.NFData (Blade p q f)
+instance (Algebra.Field.C f, Ord f, SingI p, SingI q) => Algebra.Additive.C (Multivector p q f) where
     a + b =  mvNormalForm $ BladeSum (mvTerms a ++ mvTerms b)
     a - b =  mvNormalForm $ BladeSum (mvTerms a ++ map bladeNegate (mvTerms b))
     zero = BladeSum [scalarBlade Algebra.Additive.zero]
@@ -184,7 +184,7 @@
 
 \begin{code}
 
-instance (Algebra.Field.C f, Ord f, SingI n) => Algebra.Ring.C (Multivector n f) where
+instance (Algebra.Field.C f, Ord f, SingI p, SingI q) => Algebra.Ring.C (Multivector p q f) where
     BladeSum [Blade s []] * b = BladeSum $ map (bladeScaleLeft s) $ mvTerms b
     a * BladeSum [Blade s []] = BladeSum $ map (bladeScaleRight s) $ mvTerms a 
     a * b = mvNormalForm $ BladeSum [bladeMul x y | x <- mvTerms a, y <- mvTerms b]
@@ -199,6 +199,11 @@
 
 two = fromInteger 2
 mul = (Algebra.Ring.*)
+
+psuedoScalar :: forall (p::Nat) (q::Nat) f. (Ord f, Algebra.Field.C f, SingI p, SingI q, SingI (p+q)) =>  Multivector p q f
+psuedoScalar = one `e` [1..(toNatural d)] where
+    d = fromIntegral (fromSing (sing :: Sing (p+q)) )::Word
+
 \end{code}
 
 Clifford numbers have a magnitude and absolute value:
@@ -208,18 +213,18 @@
 --magnitude :: (Algebra.Algebraic.C f) => Multivector f -> f
 magnitude = sqrt . compensatedSum' . map (\b -> (bScale b)^ 2) . mvTerms
 
-instance (Algebra.Absolute.C f, Algebra.Algebraic.C f, Ord f, SingI n) => Algebra.Absolute.C (Multivector n f) where
+instance (Algebra.Absolute.C f, Algebra.Algebraic.C f, Ord f, SingI p, SingI q) => Algebra.Absolute.C (Multivector p q f) where
     abs v =  magnitude v `e` []
     signum (BladeSum [Blade scale []]) = scalar $ signum scale 
     signum (BladeSum []) = scalar Algebra.Additive.zero
 
-instance (Algebra.Field.C f, Ord f, SingI n) => Algebra.Module.C f (Multivector n f) where
+instance (Algebra.Field.C f, Ord f, SingI p, SingI q) => Algebra.Module.C f (Multivector p q f) where
 --    (*>) zero v = Algebra.Additive.zero
     (*>) s v = v & mvTerms & map (bladeScaleLeft s) & BladeSum
 
 
 
-(/) :: (Algebra.Field.C f, Ord f, SingI n) => Multivector n f -> f -> Multivector n f
+(/) :: (Algebra.Field.C f, Ord f, SingI p, SingI q) => Multivector p q f -> f -> Multivector p q f
 (/) v d = BladeSum $ map (bladeScaleLeft (NPN.recip d)) $ mvTerms v --Algebra.Field.recip d *> v
 
 (</) n d = Numeric.Clifford.Multivector.inverse d * n
@@ -306,7 +311,7 @@
 inverse a = assert (a /= zero) $ reverseMultivector a Numeric.Clifford.Multivector./ bScale (head $ mvTerms (a * reverseMultivector a))
 recip=Numeric.Clifford.Multivector.inverse
 
-instance (Algebra.Field.C f, Ord f, SingI n) => Algebra.OccasionallyScalar.C f (Multivector n f) where
+instance (Algebra.Field.C f, Ord f, SingI p, SingI q) => Algebra.OccasionallyScalar.C f (Multivector p q f) where
     toScalar = bScale . bladeGetGrade 0 . head . mvTerms
     toMaybeScalar (BladeSum [Blade s []]) = Just s
     toMaybeScalar (BladeSum []) = Just Algebra.Additive.zero
@@ -317,7 +322,7 @@
 Also, we may as well implement the standard prelude Num interface.
 
 \begin{code}
-instance (Algebra.Algebraic.C f, SingI n,  Ord f) => PNum.Num (Multivector n f) where
+instance (Algebra.Algebraic.C f, SingI p, SingI q,  Ord f) => PNum.Num (Multivector p q f) where
     (+) = (Algebra.Additive.+)
     (-) = (Algebra.Additive.-)
     (*) = (Algebra.Ring.*)
@@ -332,23 +337,23 @@
 Let's use Newton or Halley iteration to find the principal n-th root :3
 
 \begin{code}
-root :: (Show f, Ord f, Algebra.Algebraic.C f, SingI d) => NPN.Integer -> Multivector d f -> Multivector d f
+root :: (Show f, Ord f, Algebra.Algebraic.C f, SingI p, SingI q) => NPN.Integer -> Multivector p q f -> Multivector p q f
 root n (BladeSum [Blade s []]) = scalar $ Algebra.Algebraic.root n s
 root n a@(BladeSum _) = converge $ rootIterationsStart n a one
 
-rootIterationsStart ::(Ord f, Show f, Algebra.Algebraic.C f)=>  NPN.Integer -> Multivector d f -> Multivector d f -> [Multivector d f]
+rootIterationsStart ::(Ord f, Show f, Algebra.Algebraic.C f)=>  NPN.Integer -> Multivector p q f -> Multivector p q f -> [Multivector p q f]
 rootIterationsStart n a@(BladeSum (Blade s [] :xs)) one = rootHalleysIterations n a g where
                      g = if s >= NPN.zero then one else Algebra.Ring.one `e` [1,2] --BladeSum[Blade Algebra.Ring.one [1,2]]
 rootIterationsStart n a@(BladeSum _) g = rootHalleysIterations n a g
 
 
-rootNewtonIterations :: (Algebra.Field.C f, Ord f, SingI d) => NPN.Integer -> Multivector d f -> Multivector d f -> [Multivector d f]
+rootNewtonIterations :: (Algebra.Field.C f, Ord f, SingI p, SingI q) => NPN.Integer -> Multivector p q f -> Multivector p q f -> [Multivector p q f]
 rootNewtonIterations n a = iterate xkplus1 where
                      xkplus1 xk = xk + deltaxk xk
                      deltaxk xk = oneOverN * ((Numeric.Clifford.Multivector.inverse (xk ^ (n - one))* a)  - xk)
                      oneOverN = scalar $ NPN.recip $ fromInteger n
 
-rootHalleysIterations :: (Show a, Ord a, Algebra.Algebraic.C a, SingI d) => NPN.Integer -> Multivector d a -> Multivector d a -> [Multivector d a]
+rootHalleysIterations :: (Show a, Ord a, Algebra.Algebraic.C a, SingI p, SingI q) => NPN.Integer -> Multivector p q a -> Multivector p q a -> [Multivector p q a]
 rootHalleysIterations n a = halleysMethod f f' f'' where
     f x = a - (x^n)
     f' x = fromInteger (-n) * (x^(n-1))
@@ -361,7 +366,7 @@
     down = denominator ratio
 
 
-halleysMethod :: (Show a, Ord a, Algebra.Algebraic.C a, SingI d) => (Multivector d a -> Multivector d a) -> (Multivector d a -> Multivector d a) -> (Multivector d a -> Multivector d a) -> Multivector d a -> [Multivector d a]
+halleysMethod :: (Show a, Ord a, Algebra.Algebraic.C a, SingI p, SingI q) => (Multivector p q a -> Multivector p q a) -> (Multivector p q a -> Multivector p q a) -> (Multivector p q a -> Multivector p q a) -> Multivector p q a -> [Multivector p q a]
 halleysMethod f f' f'' = iterate update where
     update x = x - (numerator x * Numeric.Clifford.Multivector.inverse (denominator x)) where
         numerator x = multiplyList [2, fx, dfx]
diff --git a/src/Numeric/Clifford/NumericIntegration.lhs b/src/Numeric/Clifford/NumericIntegration.lhs
--- a/src/Numeric/Clifford/NumericIntegration.lhs
+++ b/src/Numeric/Clifford/NumericIntegration.lhs
@@ -97,8 +97,8 @@
 
 showOutput name x = trace ("output of " ++ name ++" is " ++ show x) x
 
-convergeTolLists :: (Ord f, Algebra.Absolute.C f, Algebra.Algebraic.C f, Show f, SingI d) 
-                   => f ->  [[Multivector d f]] -> [Multivector d f]
+convergeTolLists :: (Ord f, Algebra.Absolute.C f, Algebra.Algebraic.C f, Show f, SingI p, SingI q) 
+                   => f ->  [[Multivector p q f]] -> [Multivector p q f]
 convergeTolLists t [] = error "converge: empty list"
 convergeTolLists t xs = fromMaybe empty (convergeBy check Just xs)
     where
@@ -110,17 +110,17 @@
               magnitude (sumList $ (zipWith (\x y -> NPN.abs (x-y)) b c))) <= t) = showOutput ("convergence with tolerance "++ show t )$ Just c
       check _ = Nothing
 
-type RKStepper (d::Nat) t stateType = 
-    (Ord t, Show t, Algebra.Module.C t (Multivector d t), Algebra.Field.C t) => 
+type RKStepper (p::Nat) (q::Nat) t stateType = 
+    (Ord t, Show t, Algebra.Module.C t (Multivector p q t), Algebra.Field.C t) => 
      t -> (t -> stateType -> stateType) -> 
-    ([Multivector d t] -> stateType) -> 
-    (stateType ->[Multivector d t]) ->
+    ([Multivector p q t] -> stateType) -> 
+    (stateType ->[Multivector p q t]) ->
     (t,stateType) -> (t,stateType)
 data ButcherTableau f = ButcherTableau {_tableauA :: [[f]], _tableauB :: [f], _tableauC :: [f]}
 makeLenses ''ButcherTableau
 
 
-type ConvergerFunction f = forall (d::Nat) f . [[Multivector d f]] -> [Multivector d f]
+type ConvergerFunction f = forall (p::Nat) (q::Nat) f . [[Multivector p q f]] -> [Multivector p q f]
 type AdaptiveStepSizeFunction f state = f -> state -> f 
 
 data RKAttribute f state = Explicit
@@ -135,9 +135,9 @@
 
 $( derive makeIs ''RKAttribute)
 
-genericRKMethod :: forall (d::Nat) t stateType . 
-                  ( Ord t, Show t, Algebra.Module.C t (Multivector d t),Algebra.Absolute.C t, Algebra.Algebraic.C t, SingI d)
-                  =>  ButcherTableau t -> [RKAttribute t stateType] -> RKStepper d t stateType
+genericRKMethod :: forall (p::Nat) (q::Nat) t stateType . 
+                  ( Ord t, Show t, Algebra.Module.C t (Multivector p q t),Algebra.Absolute.C t, Algebra.Algebraic.C t, SingI p, SingI q)
+                  =>  ButcherTableau t -> [RKAttribute t stateType] -> RKStepper p q t stateType
 genericRKMethod tableau attributes = rkMethodImplicitFixedPoint where
     s :: Int
     s =  length (_tableauC tableau)
@@ -151,10 +151,10 @@
     b i = l !! (i - 1) where
         l = _tableauB tableau
     
-    sumListOfLists :: [[Multivector d t]] -> [Multivector d t]
+    sumListOfLists :: [[Multivector p q t]] -> [Multivector p q t]
     sumListOfLists = map sumList . transpose 
 
-    converger :: [[Multivector d t]] -> [Multivector d t]
+    converger :: [[Multivector p q t]] -> [Multivector p q t]
     converger = case  find (\x -> isConvergenceTolerance x || isConvergenceFunction x) attributes of
                   Just (ConvergenceFunction conv) -> conv
                   Just (ConvergenceTolerance tol) -> convergeTolLists (trace ("Convergence tolerance set to " ++ show tol)tol)
@@ -165,13 +165,13 @@
                         Just (AdaptiveStepSize sigma) -> sigma
                         Nothing -> (\_ _ -> one)
 
-    rkMethodImplicitFixedPoint :: RKStepper d t stateType
+    rkMethodImplicitFixedPoint :: RKStepper p q t stateType
     rkMethodImplicitFixedPoint h f project unproject (time, state) =
         (time + (stepSizeAdapter time state)*h*(c s), newState) where
-        zi :: Int -> [Multivector d t]
+        zi :: Int -> [Multivector p q t]
         zi i = (\out -> trace ("initialGuess is " ++ show initialGuess++" whereas the final one is " ++ show out) out) $
                assert (i <= s && i>= 1) $ converger $ iterate (zkp1 i) initialGuess where
-            initialGuess :: [Multivector d t]
+            initialGuess :: [Multivector p q t]
             initialGuess = if i == 1 || null (zi (i-1)) then map (h'*>) $ unproject $ f guessTime state else zi (i-1)
             adaptiveStepSizeFraction :: t
             adaptiveStepSizeFraction = stepSizeAdapter time state
@@ -179,18 +179,18 @@
             h' = adaptiveStepSizeFraction *  h * (c i)
             guessTime :: t
             guessTime = time + h'
-            zkp1 :: NPN.Int -> [Multivector d t] -> [Multivector d t]
+            zkp1 :: NPN.Int -> [Multivector p q t] -> [Multivector p q t]
             zkp1 i zk = map (h*>) (sumOfJs i zk) where
-                sumOfJs :: Int -> [Multivector d t] -> [Multivector d t]
+                sumOfJs :: Int -> [Multivector p q t] -> [Multivector p q t]
                 sumOfJs i zk =  sumListOfLists $ map (scaledByAij zk) (a i) where 
-                    scaledByAij :: [Multivector d t] -> t -> [Multivector d t]
+                    scaledByAij :: [Multivector p q t] -> t -> [Multivector p q t]
                     scaledByAij guess a = map (a*>) $ evalDerivatives guessTime $ elementAdd state' guess
 
-        state' = unproject state :: [Multivector d t]
+        state' = unproject state :: [Multivector p q t]
         newState :: stateType
         newState = project $ elementAdd state' (assert (not $  null dy) dy)
-        dy = sumListOfLists  [map ((b i) *>) (zi i) | i <- [1..s]] :: [Multivector d t]
-        evalDerivatives :: t -> [Multivector d t] -> [Multivector d t]
+        dy = sumListOfLists  [map ((b i) *>) (zi i) | i <- [1..s]] :: [Multivector p q t]
+        evalDerivatives :: t -> [Multivector p q t] -> [Multivector p q t]
         evalDerivatives time stateAtTime= unproject $ (f time) $ project stateAtTime
 
 
diff --git a/src/Numeric/Clifford/NumericIntegration/DefaultIntegrators.hs b/src/Numeric/Clifford/NumericIntegration/DefaultIntegrators.hs
--- a/src/Numeric/Clifford/NumericIntegration/DefaultIntegrators.hs
+++ b/src/Numeric/Clifford/NumericIntegration/DefaultIntegrators.hs
@@ -36,6 +36,11 @@
 import           Test.QuickCheck
 --trace _ a = a
 
+gaussLegendreFourthOrderTableau = ButcherTableau [[0.25::NPN.Double, 0.25 - (1.0 NPN./6)* sqrt 3], [0.25 + (1.0 NPN./ 6) * sqrt 3, 0.25]] [0.5, 0.5] [0.5 - (1.0 NPN./6)* sqrt 3, 0.5 + (1.0 NPN./ 6)* sqrt 3]
+gaussLegendreFourthOrder h f (t, state) = impl h f id id (t,state) where
+    impl= genericRKMethod gaussLegendreFourthOrderTableau []
+
+
 rk4ClassicalFromTableau h f (t,state) = impl h f id id (t, state) where
     impl = genericRKMethod rk4ClassicalTableau []
 implicitEulerMethod h f (t, state) = impl h f id id (t, state) where
@@ -55,7 +60,7 @@
 lobattoIIIBFourthOrder h f (t, state) = impl h f id id (t, state) where
     impl = genericRKMethod lobattoIIIBFourthOrderTableau []
 
-rk4Classical :: (Ord a, Algebra.Algebraic.C a, SingI d) =>  stateType -> a -> (stateType->stateType) -> ([Multivector d a] -> stateType) -> (stateType -> [Multivector d a]) -> stateType
+rk4Classical :: (Ord a, Algebra.Algebraic.C a, SingI p, SingI q) =>  stateType -> a -> (stateType->stateType) -> ([Multivector p q a] -> stateType) -> (stateType -> [Multivector p q a]) -> stateType
 rk4Classical state h f project unproject = project newState where
     update = map (\(k1', k2', k3', k4') -> sumList [k1',2*k2',2*k3',k4'] MV./ Algebra.Ring.fromInteger 6) $ zip4 k1 k2 k3 k4
     newState = zipWith (+) state' update