diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) 2014, Sophie Taylor
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * 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.
+
+    * Neither the name of Sophie Taylor nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"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 COPYRIGHT
+OWNER OR CONTRIBUTORS 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/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,6 @@
+haskell-clifford
+================
+
+Clifford algebra for Haskell! :D
+
+This is initially going to just be algebraic stuff, but I'll probably add things such as numerical differentiation/integration eventually.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/clifford.cabal b/clifford.cabal
new file mode 100644
--- /dev/null
+++ b/clifford.cabal
@@ -0,0 +1,71 @@
+-- Initial clifford.cabal generated by cabal init.  For further 
+-- documentation, see http://haskell.org/cabal/users-guide/
+
+-- The name of the package.
+name:                clifford
+
+-- The package version.  See the Haskell package versioning policy (PVP) 
+-- for standards guiding when and how versions should be incremented.
+-- http://www.haskell.org/haskellwiki/Package_versioning_policy
+-- PVP summary:      +-+------- breaking API changes
+--                   | | +----- non-breaking API additions
+--                   | | | +--- code changes with no API change
+version:             0.1.0.0
+
+-- A short (one-line) description of the package.
+synopsis:            A Clifford algebra library
+
+-- A longer description of the package.
+-- description:         
+
+-- URL for the project homepage or repository.
+homepage:            http://github.com/spacekitteh/haskell-clifford
+
+-- The license under which the package is released.
+license:             BSD3
+
+-- The file containing the license text.
+license-file:        LICENSE
+
+-- The package author(s).
+author:              Sophie Taylor
+
+-- An email address to which users can send suggestions, bug reports, and 
+-- patches.
+maintainer:          sophie@traumapony.org
+
+-- A copyright notice.
+-- copyright:           
+
+category:            Math, Numerical
+
+build-type:          Simple
+
+-- Extra files to be distributed with the package, such as examples or a 
+-- README.
+extra-source-files:  README.md
+
+-- 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
+  
+  -- Modules included in this library but not exported.
+  -- other-modules:       
+  
+  -- LANGUAGE extensions used by modules in this package.
+  -- other-extensions:
+  
+  -- 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, 
+                       deepseq >= 1.3.0.1 && < 1.4, vector >= 0.10.0.1 && < 0.11, stream-fusion >= 0.1 && < 0.2, criterion >= 0.8.0.0 && < 0.9, derive, QuickCheck, nats, tagged, cereal
+  
+  -- Directories containing source files.
+  hs-source-dirs:      src
+  
+  -- Base language which the package is written in.
+  default-language:    Haskell2010
+  
diff --git a/src/Numeric/Clifford/Blade.lhs b/src/Numeric/Clifford/Blade.lhs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Clifford/Blade.lhs
@@ -0,0 +1,224 @@
+\documentclass{article}
+%include polycode.fmt
+\usepackage{fontspec}
+\usepackage{amsmath}
+\usepackage{unicode-math}
+\usepackage{lualatex-math}
+\setmainfont{latinmodern-math.otf}
+\setmathfont{latinmodern-math.otf}
+\usepackage{verbatim}
+\author{Sophie Taylor}
+\title{haskell-clifford: A Haskell Clifford algebra dynamics library}
+\begin{document}
+
+So yeah. This is a Clifford number representation. I will fill out the documentation more fully and stuff once the design has stabilised.
+
+I am basing the design of this on Issac Trotts' geometric algebra library.\cite{hga}
+
+Let us  begin. We are going to use the Numeric Prelude because it is (shockingly) nicer for numeric stuff.
+
+\begin{code}
+{-# LANGUAGE NoImplicitPrelude, FlexibleContexts, RankNTypes, ScopedTypeVariables, DeriveDataTypeable #-}
+{-# LANGUAGE NoMonomorphismRestriction, UnicodeSyntax, GADTs #-}
+{-# LANGUAGE FlexibleInstances,  UnicodeSyntax, GADTs, KindSignatures, DataKinds #-}
+{-# LANGUAGE TemplateHaskell, StandaloneDeriving #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+\end{code}
+%if False
+\begin{code}
+{-# OPTIONS_GHC -fllvm -fexcess-precision -optlo-O3 -O3 -optlc-O=3 -Wall #-}
+-- OPTIONS_GHC -Odph -fvectorise -package dph-lifted-vseg 
+--  LANGUAGE ParallelArrays
+\end{code}
+%endif
+Clifford algebras are a module over a ring. They also support all the usual transcendental functions.
+\begin{code}
+module Numeric.Clifford.Blade where
+
+import NumericPrelude hiding (iterate, head, map, tail, reverse, scanl, zipWith, drop, (++), filter, null, length, foldr, foldl1, zip, foldl, concat, (!!), concatMap,any, repeat, replicate, elem, replicate)
+import Algebra.Laws
+import Algebra.Absolute
+import Algebra.Additive
+import Algebra.Ring
+import Data.Serialize
+import Data.Word
+import Data.List.Stream
+import Data.Permute (sort, isEven)
+import Numeric.Natural
+import qualified NumericPrelude.Numeric as NPN
+import Test.QuickCheck
+import Control.Lens hiding (indices)
+import Data.DeriveTH
+import GHC.TypeLits hiding (isEven)
+import GHC.Real (fromIntegral)
+import Algebra.Field
+import Debug.Trace
+--trace _ a = a
+
+\end{code}
+
+
+The first problem: How to represent basis blades. One way to do it is via generalised Pauli matrices. Another way is to use lists, which we will do because this is Haskell. >:0
+
+\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
+
+-- makeLenses ''Blade
+scale :: Lens' (Blade n f) f
+scale = lens _scale (\blade v -> blade {_scale = v})
+indices :: Lens' (Blade n 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)
+bScale b =  b^.scale
+bIndices b = b^.indices
+instance(Show f) =>  Show (Blade n f) where
+    --TODO: Do this with HaTeX
+    show  (Blade scale indices) = pref ++  if null indices then "" else basis where
+                        pref = show scale
+                        basis =  foldr (++) "" textIndices
+                        textIndices = map vecced indices
+                        vecced index = "\\vec{e_{" ++ show index ++ "}}"
+                                                
+                        
+instance (Algebra.Additive.C f, Eq f) => Eq (Blade n f) where
+   a == b = aScale == bScale && aIndices == bIndices where
+                 (Blade aScale aIndices) = bladeNormalForm a
+                 (Blade bScale bIndices) = bladeNormalForm b
+
+\end{code}
+
+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 d = Blade d []
+
+zeroBlade :: (Algebra.Field.C f, SingI n) => Blade n f
+zeroBlade = scalarBlade Algebra.Additive.zero
+
+bladeNonZero b = b^.scale /= Algebra.Additive.zero
+
+bladeNegate b = b&scale%~negate --Blade (Algebra.Additive.negate$ b^.scale) (b^.indices)
+
+bladeScaleLeft s (Blade f ind) = Blade (s * f) ind
+bladeScaleRight s (Blade f ind) = Blade (f * s) ind
+\end{code}
+
+However, the plain data constructor should never be used, for it doesn't order them by default. It also needs to represent the vectors in an ordered form for efficiency and niceness. Further, due to skew-symmetry, if the vectors are in an odd permutation compared to the normal form, then the scale is negative. Additionally, since $\vec{e}_k^2 = 1$, pairs of them should be removed.
+
+\begin{align}
+\vec{e}_1∧...∧\vec{e}_k∧...∧\vec{e}_k∧... = 0\\
+\vec{e}_2∧\vec{e}_1 = -\vec{e}_1∧\vec{e}_2\\
+\vec{e}_k^2 = 1
+\end{align}
+
+
+\begin{code}
+bladeNormalForm :: forall (n::Nat) f. Blade n f -> Blade n 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)
+             numOfIndices = length indices
+             (sorted, perm) = Data.Permute.sort numOfIndices indices
+             scale' = if isEven perm then scale else negate scale
+             uniqueSorted = removeDupPairs sorted
+                            where
+                              removeDupPairs [] = []
+                              removeDupPairs [x] = [x]
+                              removeDupPairs (x:y:rest) | x == y = removeDupPairs rest
+                                                        | otherwise = x : removeDupPairs (y:rest)
+\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 
+
+bladeIsOfGrade :: Blade n f -> Integer -> Bool
+blade `bladeIsOfGrade` k = grade blade == k
+
+bladeGetGrade :: Integer -> Blade n f -> Blade n f
+bladeGetGrade k blade@(Blade _ _) =
+    if blade `bladeIsOfGrade` k then blade else zeroBlade
+\end{code}
+
+
+
+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
+multiplyBladeList [] = error "Empty blade list!"
+multiplyBladeList (a:[]) = a
+multiplyBladeList a = bladeNormalForm $ Blade scale indices where
+    indices = concatMap bIndices a
+    scale = foldl1 (*) (map bScale a)
+
+
+\end{code}
+
+Next up: The outer (wedge) product, denoted by $∧$ :3
+
+\begin{code}
+bWedge :: Blade n f -> Blade n f -> Blade n f
+bWedge x y = bladeNormalForm $ bladeGetGrade k xy
+             where
+               k = grade x + grade y
+               xy = bladeMul x y
+
+\end{code}
+
+Now let's do the inner (dot) product, denoted by $⋅$ :D
+
+
+\begin{code}
+bDot ::Blade n f -> Blade n f -> Blade n f
+bDot x y = bladeNormalForm $ bladeGetGrade k xy
+          where
+            k = Algebra.Absolute.abs $ grade x - grade y
+            xy = bladeMul x y
+
+propBladeDotAssociative = Algebra.Laws.associative bDot
+
+\end{code}
+
+These are the three fundamental operations on basis blades.
+
+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
+    compare a b | bIndices a == bIndices b = compare (bScale a) (bScale b)
+                | otherwise =  compare (bIndices a) (bIndices b)
+
+
+instance Arbitrary Natural where
+    arbitrary = sized $ \n ->
+                let n' = NPN.abs n in
+                 fmap (toNatural . (\x -> (GHC.Real.fromIntegral x)::Word)) (choose (0, n'))
+    shrink = shrinkIntegral
+-- $(derive makeArbitrary ''Blade)
+\end{code}
+
+Now for some testing of algebraic laws.
+
+\begin{code}
+
+{- Note: Figure out what this is meant to be lol
+skewcommutative op x y = x `op` y == (bladeScaleLeft (fromInteger (-1))$ y `op` x)
+
+propAnticommutativeMultiplication :: (Eq f,Algebra.Ring.C f, Algebra.Additive.C f) => Blade f -> Blade f -> Bool
+propAnticommutativeMultiplication = anticommutative bladeMul
+-}
+propCommutativeAddition = commutative (+)
+\end{code}
+\bibliographystyle{IEEEtran}
+\bibliography{biblio.bib}
+\end{document}
diff --git a/src/Numeric/Clifford/ClassicalMechanics.lhs b/src/Numeric/Clifford/ClassicalMechanics.lhs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Clifford/ClassicalMechanics.lhs
@@ -0,0 +1,125 @@
+\documentclass{article}
+%include polycode.fmt
+\usepackage{fontspec}
+\usepackage{amsmath}
+\usepackage{unicode-math}
+\usepackage{lualatex-math}
+\setmainfont{latinmodern-math.otf}
+\setmathfont{latinmodern-math.otf}
+\usepackage{verbatim}
+\author{Sophie Taylor}
+\title{haskell-clifford: A Haskell Clifford algebra dynamics library}
+\begin{document}
+
+This is the classical mechanics portion of the library. 
+
+\begin{code}
+{-# LANGUAGE NoImplicitPrelude, FlexibleContexts, RankNTypes, ScopedTypeVariables, DeriveDataTypeable #-}
+{-# LANGUAGE NoMonomorphismRestriction, UnicodeSyntax, GADTs, KindSignatures, DataKinds #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE TemplateHaskell #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+\end{code}
+%if False
+\begin{code}
+{-# OPTIONS_GHC -fllvm -fexcess-precision -optlo-O3 -O3 -optlc-O=3 -Wall #-}
+-- OPTIONS_GHC -Odph -fvectorise -package dph-lifted-vseg 
+--  LANGUAGE ParallelArrays
+\end{code}
+%endif
+
+\begin{code}
+module Numeric.Clifford.ClassicalMechanics where
+import Numeric.Clifford.Multivector as MV
+import Numeric.Clifford.Blade
+import GHC.TypeLits
+import Data.Proxy
+import NumericPrelude hiding (iterate, head, map, tail, reverse, scanl, zipWith, drop, (++), filter, null, length, foldr, foldl1, zip, foldl, concat, (!!), concatMap,any, repeat, replicate, elem, replicate, all)
+import Algebra.Absolute
+import Algebra.Algebraic
+import Algebra.Additive
+import Algebra.Ring
+import Algebra.ToInteger
+import Algebra.Module
+import Algebra.Field
+import Data.List.Stream
+import Numeric.Natural
+import qualified Data.Vector as V
+import NumericPrelude.Numeric (sum)
+import qualified NumericPrelude.Numeric as NPN 
+import Test.QuickCheck
+import Math.Sequence.Converge (convergeBy)
+import Number.Ratio hiding (scale)
+import Algebra.ToRational
+import Control.Lens hiding (indices)
+import Control.Exception (assert)
+import Data.Maybe
+import Data.DeriveTH
+import Data.Word
+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 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}
+
+makeLenses ''EnergyMethod
+makeLenses ''DynamicSystem
+
+--evaluateDerivative s = dq++ dp where
+--    dq = (s&energyFunction.dqs) -- s
+--    dp = (s&energyFunction.dps) -- s
+--    dq = map ($ s) ((dqs $ energyFunction) s) --s&energyFunction.dqs.traverse--map ($ s) ((dqs . energyFunction) s)
+--    dp = map ($ s) ((dps $ energyFunction) s)
+
+
+
+
+
+\end{code}
+
+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
+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}
+{-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)
+squishToDimension' d (BladeSum terms) = r' where
+    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,
+                              _mass :: f,
+                              _attitude :: PhysicalVector d f,
+                              _angularMomentum :: PhysicalVector d f,
+                              _inertia :: PhysicalVector d f
+                             } -> RigidBody d f
+
+--makeLenses ''RigidBody doesn't actually work
+{- Things to do: 
+4. create a 1-form type 
+5. figure a way to take exterior product of 1 forms at a type level so i can just go like: omega = df1 ^ df2 ^ df ; omega a b c
+-}
+
+{-data NDVector (n :: Nat) f where
+ NDVector :: (Algebra.Field.C f, Algebra.Module.C f (Multivector f)) => {value :: Multivector f} -> NDVector n f-}
+
+{-ndVector :: forall n.(n ~ Nat) => Proxy n -> (forall f.
+                  (Algebra.Field.C f, Algebra.Module.C f (Multivector f)) =>
+                  Multivector f -> NDVector (n) f)
+ndVector _ value = NDVector $ squishToDimension' (toNatural nummed) value where
+    nummed :: Word32
+    nummed = fromIntegral $ fromSing (sing :: Sing n)-}
+\end{code}
+\bibliographystyle{IEEEtran}
+\bibliography{biblio.bib}
+\end{document}
diff --git a/src/Numeric/Clifford/Multivector.lhs b/src/Numeric/Clifford/Multivector.lhs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Clifford/Multivector.lhs
@@ -0,0 +1,417 @@
+\documentclass{article}
+%include polycode.fmt
+\usepackage{fontspec}
+\usepackage{amsmath}
+\usepackage{unicode-math}
+\usepackage{lualatex-math}
+\setmainfont{latinmodern-math.otf}
+\setmathfont{latinmodern-math.otf}
+\usepackage{verbatim}
+\author{Sophie Taylor}
+\title{haskell-clifford: A Haskell Clifford algebra dynamics library}
+\begin{document}
+
+So yeah. This is a Clifford number representation. I will fill out the documentation more fully and stuff once the design has stabilised. 
+
+I am basing the design of this on Issac Trotts' geometric algebra library.\cite{hga}
+
+Let us  begin. We are going to use the Numeric Prelude because it is (shockingly) nicer for numeric stuff.
+
+\begin{code}
+{-# LANGUAGE NoImplicitPrelude, FlexibleContexts, RankNTypes, ScopedTypeVariables, DeriveDataTypeable #-}
+{-# LANGUAGE NoMonomorphismRestriction, UnicodeSyntax, GADTs#-}
+{-# LANGUAGE FlexibleInstances, StandaloneDeriving, KindSignatures, DataKinds #-}
+{-# LANGUAGE TemplateHaskell #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+\end{code}
+%if False
+\begin{code}
+{-# OPTIONS_GHC -fllvm -fexcess-precision -optlo-O3 -O3 -optlc-O=3 -Wall #-}
+-- OPTIONS_GHC -Odph -fvectorise -package dph-lifted-vseg 
+--  LANGUAGE ParallelArrays
+\end{code}
+%endif
+Clifford algebras are a module over a ring. They also support all the usual transcendental functions.
+\begin{code}
+module Numeric.Clifford.Multivector where
+import Numeric.Clifford.Blade
+import NumericPrelude hiding (iterate, head, map, tail, reverse, scanl, zipWith, drop, (++), filter, null, length, foldr, foldl1, zip, foldl, concat, (!!), concatMap,any, repeat, replicate, elem, replicate, all)
+--import Algebra.Laws
+import Algebra.Absolute
+import Algebra.Algebraic
+import Algebra.Additive
+import Algebra.Ring
+import Algebra.OccasionallyScalar
+import Algebra.ToInteger
+import Algebra.Transcendental
+import Algebra.ZeroTestable
+import Algebra.Module
+import Algebra.Field
+import Data.Serialize
+import MathObj.Polynomial.Core (progression)
+import System.IO
+import Data.List.Stream
+import Data.Permute (sort, isEven)
+import Data.List.Ordered
+import Data.Ord
+import Data.Maybe
+--import Number.NonNegative
+import Numeric.Natural
+import qualified Data.Vector as V
+import NumericPrelude.Numeric (sum)
+import qualified NumericPrelude.Numeric as NPN
+import Test.QuickCheck
+import Math.Sequence.Converge (convergeBy)
+import Control.DeepSeq 
+import Number.Ratio hiding (scale)
+import Algebra.ToRational
+import qualified GHC.Num as PNum
+import Control.Lens hiding (indices)
+import Control.Exception (assert)
+import Data.Maybe
+import Data.Monoid
+import Data.Data
+import Data.DeriveTH
+import GHC.TypeLits
+import Data.Word
+import Debug.Trace
+--trace _ a = a
+
+\end{code}
+
+
+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
+
+deriving instance Eq (Multivector n f)
+deriving instance Ord (Multivector n f)
+deriving instance (Show f) => Show (Multivector n f)
+
+dimension :: forall (n::Nat) f. SingI n => Multivector n f ->  Natural
+dimension _ = toNatural  ((fromIntegral $ fromSing (sing :: Sing n))::Word)
+
+terms :: Lens' (Multivector n f) [Blade n f]
+terms = lens _terms (\bladeSum v -> bladeSum {_terms = v})
+
+mvNormalForm (BladeSum terms) = BladeSum $ if null resultant then [scalarBlade Algebra.Additive.zero] else resultant  where
+    resultant = filter bladeNonZero $ addLikeTerms' $ Data.List.Ordered.sortBy compare $  map bladeNormalForm $ terms
+mvTerms m = m^.terms
+
+addLikeTerms' = sumLikeTerms . groupLikeTerms
+
+groupLikeTerms ::Eq f =>  [Blade n f] -> [[Blade n f]]
+groupLikeTerms = groupBy (\a b -> a^.indices == b^.indices)
+
+compensatedSum' :: (Algebra.Additive.C f) => [f] -> f
+compensatedSum' xs = kahan zero zero xs where
+    kahan s _ [] = s
+    kahan s c (x:xs) = 
+        let y = x - c
+            t = s + y
+        in kahan t ((t-s)-y) xs
+
+--use this to sum taylor series et al with converge
+--compensatedRunningSum :: (Algebra.Additive.C f) => [f] -> [f]
+compensatedRunningSum xs=shanksTransformation . map fst $ scanl kahanSum (zero,zero) xs where
+    kahanSum (s,c) b = (t,newc) where
+        y = b - c
+        t = s + y
+        newc = (t - s) - y
+            
+--multiplyAdd a b c = a*b + c
+--twoProduct a b = (x,y) where
+--    x = a*b
+--z    y = multiplyAdd a b (negate x)
+
+--multiplyList [] = []
+--multiplyList a@(x:[])=a
+--multiplyList (a:b:xs) = loop a (b:xs) zero where
+--  loop pm [] ei = pm+ei
+--  loop pm1 (ai:remaining) eim1= loop pi remaining ei where
+--      (pi, pii) = twoProduct pm1 ai
+--      ei = multiplyAdd eim1 ai pii
+
+
+multiplyOutBlades :: (SingI n, Algebra.Ring.C a) => [Blade n a] -> [Blade n a] -> [Blade n a]
+multiplyOutBlades x y = [bladeMul l r | l <-x, r <- y]
+
+--multiplyList :: Algebra.Ring.C t => [Multivector t] -> Multivector t
+multiplyList [] = error "Empty list"
+--multiplyList a@(x:[]) = x
+multiplyList l = mvNormalForm $ BladeSum listOfBlades where
+    expandedBlades a = foldl1 multiplyOutBlades a
+    listOfBlades = expandedBlades $ map mvTerms l
+
+--things to test: is 1. adding blades into a map based on indices 2. adding errything together 3. sort results quicker than
+--                   1. sorting by indices 2. groupBy-ing on indices 3. adding the lists of identical indices
+
+sumList xs = mvNormalForm $ BladeSum $ concat $ map mvTerms xs
+
+sumLikeTerms :: (Algebra.Field.C f, SingI n) => [[Blade n f]] -> [Blade n 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
+    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
+    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
+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
+    a + b =  mvNormalForm $ BladeSum (mvTerms a ++ mvTerms b)
+    a - b =  mvNormalForm $ BladeSum (mvTerms a ++ map bladeNegate (mvTerms b))
+    zero = BladeSum [scalarBlade Algebra.Additive.zero]
+
+    
+\end{code}
+
+Now it is time for the Clifford product. :3
+
+\begin{code}
+
+instance (Algebra.Field.C f, Ord f, SingI n) => Algebra.Ring.C (Multivector n 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]
+    one = scalar Algebra.Ring.one
+    fromInteger i = scalar $ Algebra.Ring.fromInteger i    
+
+    a ^ 2 = a * a
+    a ^ 0 = one
+    a ^ 1 = a
+    --a ^ n  --n < 0 = Clifford.recip $ a ^ (negate n)
+    a ^ n  =  multiplyList (replicate (NPN.fromInteger n) a)
+
+two = fromInteger 2
+mul = (Algebra.Ring.*)
+\end{code}
+
+Clifford numbers have a magnitude and absolute value:
+
+\begin{code}
+
+--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
+    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
+--    (*>) 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
+(/) v d = BladeSum $ map (bladeScaleLeft (NPN.recip d)) $ mvTerms v --Algebra.Field.recip d *> v
+
+(</) n d = Numeric.Clifford.Multivector.inverse d * n
+(/>) n d = n * Numeric.Clifford.Multivector.inverse d
+(</>) n d = n /> d
+
+integratePoly c x = c : zipWith (Numeric.Clifford.Multivector./) x progression
+
+--converge :: (Eq f, Show f) => [f] -> f
+converge [] = error "converge: empty list"
+converge xs = fromMaybe empty (convergeBy checkPeriodic Just xs) 
+    where
+      empty = error "converge: error in implmentation"
+      checkPeriodic (a:b:c:_)
+          | (trace ("Converging at " ++ show a) a) == b = Just a
+          | a == c = Just a
+      checkPeriodic _ = Nothing
+
+
+
+aitkensAcceleration [] = []
+aitkensAcceleration a@(xn:[]) = a
+aitkensAcceleration a@(xn:xnp1:[]) = a
+aitkensAcceleration a@(xn:xnp1:xnp2:[]) = a
+aitkensAcceleration (xn:xnp1:xnp2:xs) | xn == xnp1 = [xnp1]
+                                      | xn == xnp2 = [xnp2]
+                                      | otherwise = xn - ((dxn ^ 2) /> ddxn) : aitkensAcceleration (xnp1:xnp2:xs) where
+    dxn = sumList [xnp1,negate xn]
+    ddxn = sumList [xn,  (-2) *  xnp1, xnp2]
+
+shanksTransformation [] = []
+shanksTransformation a@(xnm1:[]) = a
+shanksTransformation a@(xnm1:xn:[]) = a
+shanksTransformation (xnm1:xn:xnp1:xs) | xnm1 == xn = [xn]
+                                       | xnm1 == xnp1 = [xnm1]
+                                       | denominator == zero = [xnp1]
+                                       | otherwise = trace ("Shanks transformation input = " ++ show xn ++ "\nShanks transformation output = " ++ show out) out:shanksTransformation (xn:xnp1:xs) where
+                                       out = numerator />  denominator 
+                                       numerator = sumList [xnp1*xnm1, negate (xn^2)]
+                                       denominator = sumList [xnp1, (-2)*xn, xnm1] 
+
+
+--exp ::(Ord f, Show f, Algebra.Transcendental.C f)=> Multivector f -> Multivector f
+exp (BladeSum [ Blade s []]) = trace ("scalar exponential of " ++ show s) scalar $ Algebra.Transcendental.exp s
+exp x = trace ("Computing exponential of " ++ show x) convergeTerms x where --(expMag ^ expScaled) where
+    --todo: compute a ^ p via a^n where n = floor p then multiply remaining power
+    expMag = Algebra.Transcendental.exp mag
+    expScaled = converge $ shanksTransformation.shanksTransformation . compensatedRunningSum $ expTerms scaled 
+    convergeTerms terms = converge $ shanksTransformation.shanksTransformation.compensatedRunningSum $ expTerms terms
+    mag = trace ("In exponential, magnitude is " ++ show ( magnitude x)) magnitude x
+    scaled = let val = (Numeric.Clifford.Multivector./) x mag in trace ("In exponential, scaled is" ++ show val) val
+
+
+
+
+
+takeEvery nth xs = case drop (nth-1) xs of
+                     (y:ys) -> y : takeEvery nth ys
+                     [] -> []
+
+cosh x = converge $ shanksTransformation . compensatedRunningSum $ takeEvery 2 $ expTerms x
+
+sinh x = converge $ shanksTransformation . compensatedRunningSum $ takeEvery 2 $ tail $ expTerms x
+
+seriesPlusMinus (x:y:rest) = x:Algebra.Additive.negate y: seriesPlusMinus rest
+seriesMinusPlus (x:y:rest) = Algebra.Additive.negate x : y : seriesMinusPlus rest
+
+
+sin x = converge $ shanksTransformation $ compensatedRunningSum $ sinTerms x
+sinTerms x = seriesPlusMinus $ takeEvery 2 $ expTerms x
+cos x = converge $ shanksTransformation $ compensatedRunningSum (Algebra.Ring.one : cosTerms x)
+cosTerms x = seriesMinusPlus $ takeEvery 2 $ tail $ expTerms x
+
+expTerms x = map snd $ iterate (\(n,b) -> (n + 1, (x*b) Numeric.Clifford.Multivector./ fromInteger n )) (1::NPN.Integer,one)
+
+dot a b = mvNormalForm $ BladeSum [x `bDot` y | x <- mvTerms a, y <- mvTerms b]
+wedge a b = mvNormalForm $ BladeSum [x `bWedge` y | x <- mvTerms a, y <- mvTerms b]
+(∧) = wedge
+(⋅) = dot
+
+reverseBlade b = bladeNormalForm $ b & indices %~ reverse 
+reverseMultivector v = mvNormalForm $ v & terms.traverse%~ reverseBlade
+
+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
+    toScalar = bScale . bladeGetGrade 0 . head . mvTerms
+    toMaybeScalar (BladeSum [Blade s []]) = Just s
+    toMaybeScalar (BladeSum []) = Just Algebra.Additive.zero
+    toMaybeScalar _ = Nothing
+    fromScalar = scalar
+\end{code}
+
+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
+    (+) = (Algebra.Additive.+)
+    (-) = (Algebra.Additive.-)
+    (*) = (Algebra.Ring.*)
+    negate = NPN.negate
+    abs = scalar . magnitude 
+    fromInteger = Algebra.Ring.fromInteger
+    signum m = Numeric.Clifford.Multivector.inverse (scalar $ magnitude m) * m
+
+
+\end{code}
+ 
+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 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 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 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 n a = halleysMethod f f' f'' where
+    f x = a - (x^n)
+    f' x = fromInteger (-n) * (x^(n-1))
+    f'' x = fromInteger (-(n*(n-1))) * (x^(n-2))
+
+
+pow a p = (a ^ up) Numeric.Clifford.Multivector./> Numeric.Clifford.Multivector.root down a where
+    ratio = toRational p
+    up = numerator ratio
+    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 f f' f'' = iterate update where
+    update x = x - (numerator x * Numeric.Clifford.Multivector.inverse (denominator x)) where
+        numerator x = multiplyList [2, fx, dfx]
+        denominator x = multiplyList [2, dfx, dfx] - (fx * ddfx)
+        fx = f x
+        dfx = f' x
+        ddfx = f'' x
+
+
+secantMethod f x0 x1 = update x1 x0  where
+    update xm1 xm2 | xm1 == xm2 = [xm1]
+                   | otherwise = if x == xm1 then [x] else x : update x xm1 where
+      x = xm1 - f xm1 * (xm1-xm2) * Numeric.Clifford.Multivector.inverse (f xm1 - f xm2)
+
+
+\end{code}
+
+Now let's try logarithms by fixed point iteration. It's gonna be slow, but whatever!
+
+\begin{code}
+
+normalised a = a * (scalar $ NPN.recip $ magnitude a)
+
+log (BladeSum [Blade s []]) = scalar $ NPN.log s
+log a = scalar (NPN.log mag) + log' scaled where
+    scaled = normalised a
+    mag = magnitude a
+    log' a = converge $  halleysMethod f f' f'' (one `e` [1,2])  where
+         f x = a - Numeric.Clifford.Multivector.exp x
+         f' x = NPN.negate $ Numeric.Clifford.Multivector.exp x
+         f'' = f'
+\end{code}
+
+Now let's do (slow as fuck probably) numerical integration! :D~! Since this is gonna be used for physical applications, it's we're gonna start off with a Hamiltonian structure and then a symplectic integrator.
+
+\begin{code}
+
+
+
+
+{- $(derive makeSerialize ''Blade)
+$(derive makeSerialize ''Multivector)
+$(derive makeData ''Blade)
+$(derive makeTypeable ''Blade)
+$(derive makeData ''Multivector)
+$(derive makeTypeable ''Multivector)-}
+
+-- $(derive makeArbitrary ''Multivector)
+
+\end{code}
+\bibliographystyle{IEEEtran}
+\bibliography{biblio.bib}
+\end{document}
diff --git a/src/Numeric/Clifford/NumericIntegration.lhs b/src/Numeric/Clifford/NumericIntegration.lhs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Clifford/NumericIntegration.lhs
@@ -0,0 +1,212 @@
+\documentclass{article}
+%include polycode.fmt
+\usepackage{fontspec}
+\usepackage{amsmath}
+\usepackage{unicode-math}
+\usepackage{lualatex-math}
+\setmainfont{latinmodern-math.otf}
+\setmathfont{latinmodern-math.otf}
+\usepackage{verbatim}
+\author{Sophie Taylor}
+\title{haskell-clifford: A Haskell Clifford algebra dynamics library}
+\begin{document}
+
+This is the numeric integration portion of the library. 
+
+\begin{code}
+{-# LANGUAGE NoImplicitPrelude, FlexibleContexts, RankNTypes, ScopedTypeVariables, DeriveDataTypeable #-}
+{-# LANGUAGE NoMonomorphismRestriction, UnicodeSyntax, GADTs, DataKinds, KindSignatures #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE TemplateHaskell #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+\end{code}
+%if False
+\begin{code}
+{-# OPTIONS_GHC -fllvm -fexcess-precision -optlo-O3 -O3 -optlc-O=3 -Wall #-}
+-- OPTIONS_GHC -Odph -fvectorise -package dph-lifted-vseg 
+--  LANGUAGE ParallelArrays
+\end{code}
+%endif
+
+\begin{code}
+
+module Numeric.Clifford.NumericIntegration where
+import Numeric.Clifford.Multivector as MV
+import NumericPrelude hiding (iterate, head, map, tail, reverse, scanl, zipWith, drop, (++), filter, null, length, foldr, foldl1, zip, foldl, concat, (!!), concatMap,any, repeat, replicate, elem, replicate)
+import Algebra.Absolute
+import Algebra.Algebraic
+import Algebra.Additive hiding (elementAdd, elementSub)
+import Algebra.Ring
+import Data.Monoid
+import Algebra.ToInteger
+import Algebra.Module
+import Algebra.Field
+import Data.List.Stream
+import Numeric.Natural
+import qualified Data.Vector as V
+import NumericPrelude.Numeric (sum)
+import qualified NumericPrelude.Numeric as NPN 
+import Test.QuickCheck
+import Math.Sequence.Converge (convergeBy)
+import Number.Ratio hiding (scale)
+import Algebra.ToRational
+import Control.Lens hiding (indices)
+import Control.Exception (assert)
+import Data.Maybe
+import GHC.TypeLits
+import Data.DeriveTH
+import Debug.Trace
+--trace _ a = a
+
+
+elementAdd = zipWith (+)
+elementScale = zipWith (*>) 
+a `elementSub` b = zipWith (-) a b
+a `elementMul` b = zipWith (*) a b
+
+
+--rk4ClassicalTableau :: ButcherTableau NPN.Double
+rk4ClassicalTableau = ButcherTableau [[0,0,0,0],[0.5,0,0,0],[0,0.5,0,0],[0,0,1,0]] [1.0 NPN./6,1.0 NPN./3, 1.0 NPN./3, 1.0 NPN./6] [0, 0.5, 0.5, 1]
+implicitEulerTableau = ButcherTableau [[1.0::NPN.Double]] [1] [1]
+
+
+
+sumVector = sumList . V.toList 
+
+--systemRootSolver :: [Multivector f] -> [Multivector f] -> ratio -> [Multivector f] -> [Multivector f] -> [Multivector f]
+
+--This will stop as soon as one of the elements converges. This is bad. Need to make it skip convergent ones and focus on the remainig.
+systemBroydensMethod f x0 x1 = map fst $ update (x1,ident) x0  where
+    update (xm1,jm1) xm2 | zero `elem` dx =  [(xm1,undefined)]
+                   | otherwise = if x == xm1 then [(x,undefined)] else (x,j) : update (x,j) xm1 where
+      x = xm1 `elementSub` ( (fm1 `elementMul` dx) `elementMul` ody)
+      j = undefined
+      fm1 = f xm1
+      fm2 = f xm2
+      dx = elementSub xm1 xm2
+      dy = elementSub fm1 fm2
+      ody = map inverse dy
+    ident = undefined
+
+--TODO: implement Broyden-Fletcher-Goldfarb-Shanno method
+
+
+convergeList ::(Show f, Ord f) => [[f]] -> [f]
+convergeList = converge
+
+
+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 t [] = error "converge: empty list"
+convergeTolLists t xs = fromMaybe empty (convergeBy check Just xs)
+    where
+      empty = error "converge: error in impl"
+      check (a:b:c:_)
+          | (trace ("Converging at " ++ show a) a) == b = Just b
+          | a == c = Just c
+          | ((showOutput ("convergence check with tolerance " ++ show t) $ 
+              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) => 
+     t -> (t -> stateType -> stateType) -> 
+    ([Multivector d t] -> stateType) -> 
+    (stateType ->[Multivector d 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 AdaptiveStepSizeFunction f state = f -> state -> f 
+
+data RKAttribute f state = Explicit
+                 | HamiltonianFunction
+                 | AdaptiveStepSize {sigma :: AdaptiveStepSizeFunction f state}
+                 | ConvergenceTolerance {epsilon :: f}
+                 | ConvergenceFunction {converger :: ConvergerFunction f } 
+                 | RootSolver 
+                 | UseAutomaticDifferentiationForRootSolver
+                 | StartingGuessMethod 
+
+
+$( 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 tableau attributes = rkMethodImplicitFixedPoint where
+    s :: Int
+    s =  length (_tableauC tableau)
+    c :: Int -> t
+    c n = l !!  (n-1) where
+        l = _tableauC tableau
+    a :: Int -> [t]
+    a n = (l !! (n-1)) & filter (/= zero) where
+        l = _tableauA tableau
+    b :: Int -> t
+    b i = l !! (i - 1) where
+        l = _tableauB tableau
+    
+    sumListOfLists :: [[Multivector d t]] -> [Multivector d t]
+    sumListOfLists = map sumList . transpose 
+
+    converger :: [[Multivector d t]] -> [Multivector d 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)
+                  Nothing -> trace "No convergence tolerance specified, defaulting to equality" convergeList 
+
+    stepSizeAdapter :: AdaptiveStepSizeFunction t stateType
+    stepSizeAdapter = case find isAdaptiveStepSize attributes of
+                        Just (AdaptiveStepSize sigma) -> sigma
+                        Nothing -> (\_ _ -> one)
+
+    rkMethodImplicitFixedPoint :: RKStepper d t stateType
+    rkMethodImplicitFixedPoint h f project unproject (time, state) =
+        (time + (stepSizeAdapter time state)*h*(c s), newState) where
+        zi :: Int -> [Multivector d 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 = if i == 1 || null (zi (i-1)) then map (h'*>) $ unproject $ f guessTime state else zi (i-1)
+            adaptiveStepSizeFraction :: t
+            adaptiveStepSizeFraction = stepSizeAdapter time state
+            h' :: t
+            h' = adaptiveStepSizeFraction *  h * (c i)
+            guessTime :: t
+            guessTime = time + h'
+            zkp1 :: NPN.Int -> [Multivector d t] -> [Multivector d t]
+            zkp1 i zk = map (h*>) (sumOfJs i zk) where
+                sumOfJs :: Int -> [Multivector d t] -> [Multivector d t]
+                sumOfJs i zk =  sumListOfLists $ map (scaledByAij zk) (a i) where 
+                    scaledByAij :: [Multivector d t] -> t -> [Multivector d t]
+                    scaledByAij guess a = map (a*>) $ evalDerivatives guessTime $ elementAdd state' guess
+
+        state' = unproject state :: [Multivector d 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]
+        evalDerivatives time stateAtTime= unproject $ (f time) $ project stateAtTime
+
+
+\end{code}
+
+
+
+
+
+
+Look at creating an exponential integrator: https://en.wikipedia.org/wiki/Exponential_integrators
+
+
+
+
+
+\bibliographystyle{IEEEtran}
+\bibliography{biblio.bib}
+\end{document}
diff --git a/src/Numeric/Clifford/NumericIntegration/DefaultIntegrators.hs b/src/Numeric/Clifford/NumericIntegration/DefaultIntegrators.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Clifford/NumericIntegration/DefaultIntegrators.hs
@@ -0,0 +1,69 @@
+{-# LANGUAGE NoImplicitPrelude #-}
+module Numeric.Clifford.NumericIntegration.DefaultIntegrators where
+import           Algebra.Absolute
+import           Algebra.Additive                    hiding (elementAdd,
+                                                      elementSub)
+import           Algebra.Algebraic
+import           Algebra.Field
+import           Algebra.Module
+import           Algebra.Ring
+import           Algebra.ToInteger
+import           Algebra.ToRational
+import           Control.Exception                   (assert)
+import           Control.Lens                        hiding (indices)
+import           Data.DeriveTH
+import           Data.List.Stream
+import           Data.Maybe
+import           Data.Monoid
+import qualified Data.Vector                         as V
+import           Debug.Trace
+import           GHC.TypeLits
+import           Math.Sequence.Converge              (convergeBy)
+import           Number.Ratio                        hiding (scale)
+import           Numeric.Clifford.Multivector        as MV
+import           Numeric.Clifford.NumericIntegration
+import           Numeric.Natural
+import           NumericPrelude                      hiding (any, concat,
+                                                      concatMap, drop, elem,
+                                                      filter, foldl, foldl1,
+                                                      foldr, head, iterate,
+                                                      length, map, null, repeat,
+                                                      replicate, replicate,
+                                                      reverse, scanl, tail, zip,
+                                                      zipWith, (!!), (++))
+import           NumericPrelude.Numeric              (sum)
+import qualified NumericPrelude.Numeric              as NPN
+import           Test.QuickCheck
+--trace _ a = a
+
+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
+    impl = genericRKMethod implicitEulerTableau []
+
+lobattoIIIASecondOrderTableau = ButcherTableau [[0,0],[0.5::NPN.Double,0.5]] [0.5,0.5] [0,1]
+lobattoIIIASecondOrder h f (t, state) = impl h f id id  (t, state) where
+    impl = genericRKMethod lobattoIIIASecondOrderTableau []
+
+lobattoIIIAFourthOrderWithTol h f (t, state) = impl h f id id (t, state) where
+    impl = genericRKMethod lobattoIIIAFourthOrderTableau [ConvergenceTolerance 1.0e-8]
+lobattoIIIAFourthOrderTableau = ButcherTableau [[0,0,0],[((5 NPN./24)::NPN.Double),1 NPN./3,-1 NPN./24],[1 NPN./6,2 NPN./3,1 NPN./6]] [1 NPN./6,2 NPN./3,1 NPN./6] [0,0.5,1]
+lobattoIIIAFourthOrder h f (t, state) = impl h f id id (t, state) where
+    impl = genericRKMethod lobattoIIIAFourthOrderTableau []
+
+lobattoIIIBFourthOrderTableau = ButcherTableau [[1 NPN./6,(-1) NPN./6,0],[((1 NPN./6)::NPN.Double),1 NPN./3,0],[1 NPN./6,5 NPN./6, 0]] [1 NPN./6,2 NPN./3,1 NPN./6] [0,0.5,1]
+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 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
+    state' = unproject state
+    evalDerivatives x = unproject $ f $ project x
+    k1 = map (h*>) $ evalDerivatives state'
+    k2 = map (h*>) $ evalDerivatives . map (uncurry (+)) $ zip state' (map (MV./ two) k1)
+    k3 = map (h*>) $ evalDerivatives . map (uncurry (+)) $ zip state' (map (MV./ two) k2)
+    k4 = map (h*>) $ evalDerivatives . map (uncurry (+)) $ zip state' k3
+
+rk4ClassicalList state h f = rk4Classical state h f id id