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/np-linear.cabal b/np-linear.cabal
new file mode 100644
--- /dev/null
+++ b/np-linear.cabal
@@ -0,0 +1,30 @@
+name:                np-linear
+version:             0.1.1.1
+synopsis:            Linear algebra for the numeric-prelude framework
+-- description:         
+license:             BSD3
+author:              Arie Peterson
+maintainer:          ariep@xs4all.nl
+category:            Math
+stability:           experimental
+build-type:          Simple
+cabal-version:       >=1.10
+
+library
+  exposed-modules:
+    Algebra.Linear,
+    Algebra.Linear.Integral,
+    Algebra.Linear.Subspace,
+    Algebra.Module.Free
+  other-modules:
+    Auxiliary
+  build-depends:
+    base >= 4.5 && < 4.7,
+    binary >= 0.6.3 && < 0.8,
+    containers >= 0.5 && < 0.6,
+    numeric-prelude >= 0.3 && < 0.5,
+    -- bifunctors >= 4.1 && < 0.5,
+    reflection >= 1.3 && < 1.5,
+    tagged == 0.7.*
+  hs-source-dirs:      src
+  default-language:    Haskell2010
diff --git a/src/Algebra/Linear.hs b/src/Algebra/Linear.hs
new file mode 100644
--- /dev/null
+++ b/src/Algebra/Linear.hs
@@ -0,0 +1,228 @@
+{-# LANGUAGE NoImplicitPrelude,RebindableSyntax #-}
+{-# LANGUAGE ScopedTypeVariables,FlexibleContexts #-}
+module Algebra.Linear
+  (
+    Relation(..)
+  , satisfies
+  , solve
+  , dependencies
+  , equations
+  , inverseImage
+  
+  , Matrix
+  , Vector
+  , matrixProduct
+  , matrixVector
+  , innerProduct
+  
+  , identity
+  , matrixFromFunction
+  , affine
+  , permute
+  , rowSwap
+  
+  , invert
+  , determinant
+  , adjoint
+  , diagonal
+  ) where
+
+
+import           Auxiliary (δ,findAmong,headE,adorn)
+
+import qualified Algebra.Field
+import qualified Algebra.Lattice
+import qualified Algebra.Module
+import qualified Algebra.Ring
+import           NumericPrelude
+
+import           Control.Applicative ((<$>))
+import           Control.Arrow       (first,second)
+import           Data.List
+  (
+    partition
+  , find
+  , transpose
+  , genericLength
+  , genericTake
+  , genericDrop
+  , genericReplicate
+  )
+import qualified Data.Map as Map
+import           Data.Proxy          (Proxy(Proxy))
+import           Data.Reflection     (Reifies,reflect)
+
+
+type Vector k
+  = [k]
+
+type Matrix k
+  = [Vector k]
+
+-- A relation among vectors of degree 'd' over field 'k'
+data Relation k
+  = Relation { getRelation :: [k] }
+  deriving (Show)
+
+satisfies :: (Algebra.Ring.C k,Eq k) => Vector k -> Relation k -> Bool
+satisfies v (Relation r) = innerProduct v r == 0
+
+-- | Calculate all dependencies among the given vectors of degree d.
+dependencies :: (Algebra.Field.C k,Eq k) => Integer -> [Vector k] -> [Relation k]
+dependencies d = map (Relation . genericDrop d) . filter (all (== zero) . genericTake d) . (\ (r,_,_) -> r) . reduce . adorn
+
+-- | Calculate the equations satisfied by the subspace spanned by the given vectors of degree d.
+equations :: (Algebra.Field.C k,Eq k) => Integer -> [Vector k] -> [Relation k]
+equations d vs = dependencies (genericLength vs + 1) . transpose . (:) (genericReplicate d zero) $ vs
+
+-- | Solve the given equations, in the form of a basis of vectors of degree d.
+solve :: (Algebra.Field.C k,Eq k) => Integer -> [Relation k] -> [Vector k]
+solve d = map getRelation . equations d . map getRelation
+
+inverseImage :: (Algebra.Field.C k,Eq k) => Matrix k -> Vector k -> Vector k
+inverseImage a = solveUpperTriangular u . matrixVector b where
+  (u,b,_) = reduce a
+
+-- Gives the row-reduced matrix, together with the determinant of the applied row operations.
+-- rowReduce :: forall k. (Algebra.Field.C k,Eq k) => Matrix k -> (Matrix k,k)
+-- rowReduce []          = ([],1)
+-- rowReduce xs@([] : _) = (xs,1)
+-- rowReduce vs          = case findAmong ((/= 0) . headE "rr1" . snd) ivs of
+--   Nothing            -> recurse (map (tail . snd) ivs)
+--   Just ((i,v₀),rest) -> first (v₀' :) . second (* s) $ recurse vs' where
+--     λ   = headE "rr2" v₀
+--     v₀' = map (/ λ) v₀
+--     vs' = map (tail . f . snd) $ rest
+--     f v = zipWith (\ x y -> x - c * y) v v₀' where
+--       c = headE "rr3" v
+--     s = (if i == 0 then id else negate) λ
+--  where
+--   ivs = zip [0 ..] vs :: [(Int,Vector k)]
+--   recurse = first (map (0 :)) . rowReduce
+
+invert :: (Algebra.Field.C k,Eq k) => Matrix k -> Maybe (Matrix k)
+invert m = fmap strip . process . (\ (r,_,_) -> r) . reduce . adorn $ m where
+  n = length m
+  process x = go x where
+    go [v] = Just [v]
+    go (r@(pivot : r') : rest)
+      | pivot == 1 = do
+        m' <- go (map tail rest)
+        rowsToAdd <- sequence . zipWith3 f [0 ..] r' $ m'
+        return $ (1 : foldr (zipWith (+)) r' rowsToAdd) : map (0 :) m'
+      | otherwise = Nothing
+    f i c v
+      | v₀ == 0   = Nothing
+      | otherwise = Just $ map ((*) (negate c / v₀)) v
+      where
+        v₀ = v !! i
+  strip = map (drop n)
+
+determinant :: (Algebra.Field.C k,Eq k) => Matrix k -> k
+determinant [] = one
+determinant a
+  | m == n = product (diagonal rr) * σ
+ where
+  (rr,_,σ) = reduce a
+  m = length a
+  n = length (head a)
+
+adjoint :: (Algebra.Lattice.C k,Algebra.Field.C k,Eq k) => Matrix k -> Maybe (Matrix k)
+adjoint a = (abs (determinant a) *>) <$> invert a
+
+diagonal :: Matrix k -> [k]
+diagonal []               = []
+diagonal ((d : _) : rest) = d : diagonal (map tail rest)
+
+matrixProduct :: (Algebra.Ring.C k) => Matrix k -> Matrix k -> Matrix k
+matrixProduct a b = map (($ transpose b) . map . innerProduct) a 
+
+matrixVector :: (Algebra.Ring.C k) => Matrix k -> Vector k -> Vector k
+matrixVector a = map (\ [x] -> x) . matrixProduct a . map (: [])
+
+innerProduct :: (Algebra.Ring.C k) => Vector k -> Vector k -> k
+innerProduct a b = sum (zipWith (*) a b)
+
+matrixFromFunction :: Int -> Int -> (Int -> Int -> k) -> Matrix k
+matrixFromFunction m n f = [[f i j | j <- [1 .. n]] | i <- [1 .. m]]
+
+identity :: (Algebra.Ring.C k) => Int -> Matrix k
+identity n = matrixFromFunction n n δ
+
+
+-- Triangular decomposition
+
+solveUpperTriangular :: (Algebra.Field.C k,Eq k) => Matrix k -> Vector k -> Vector k
+solveUpperTriangular u x = reverse $ solveLowerTriangular (reverse . map reverse $ u) (reverse x)
+
+solveLowerTriangular :: (Algebra.Field.C k,Eq k) => Matrix k -> Vector k -> Vector k
+solveLowerTriangular l x = go l x [] where
+  go []        []        q = q
+  go (lᵢ : l') (xᵢ : x') q = go l' x' (q ++ [qᵢ]) where
+    (lFirst,lᵢᵢ : _) = splitAt (length q) lᵢ
+    y = xᵢ - innerProduct q lFirst
+    qᵢ
+      | lᵢᵢ == 0  = error "Linear.solveLowerTriangular: zero on diagonal"
+      | otherwise = y / lᵢᵢ
+
+-- Compute the row echelon form of the matrix,
+-- together with the basis transformation matrix,
+-- and its determinant.
+reduce :: (Algebra.Field.C k,Eq k) => Matrix k -> (Matrix k,Matrix k,k)
+reduce []          = ([],[],1)
+reduce xs@([] : _) = (xs,identity (length xs),1)
+reduce vs          = case nonZero of
+  []                    -> (\ (x,u,σ) -> (map (0 :) x,u,σ)) $ reduce (map tail vs)
+  (v@(v₀ : _),i) : []   -> let
+    (h,u,σ) = reduce (map (tail . fst) startZero)
+   in
+    ( (map (/ v₀) v :) . map (0 :) $ h
+    , normalisation v₀ `matrixProduct` rowSwap n (1,i) `matrixProduct` shift u
+    ,v₀ * σ
+    )
+  (v@(v₀ : _),i) : rest -> let
+    (reduced,translates) = unzip . flip map rest $ \ (x@(x₀ : _),j) -> let
+      c = x₀ / v₀ in ((zipWith (\ vᵢ xᵢ -> xᵢ - c * vᵢ) v x,j),(j,c))
+    (h,u,σ) = reduce $ v : map fst reduced ++ map fst startZero
+   in
+    ( h
+    , u
+        `matrixProduct` permute (i : map snd reduced ++ map snd startZero)
+        `matrixProduct` affine n i (map (second negate) translates)
+    , σ
+    )
+ where
+  (startZero,nonZero) = partition ((==) 0 . head . fst) . flip zip [1 ..] $ vs
+  normalisation v₀ = matrixFromFunction n n f where
+    f 1 1 = 1 / v₀
+    f i j = δ i j
+  shift u = (1 : replicate (n - 1) 0) : map (0 :) u
+  n = length vs
+
+rowSwap :: (Algebra.Ring.C k) => Int -> (Int,Int) -> Matrix k
+rowSwap n (k,l) = matrixFromFunction n n f where
+  f i j
+    | i == k && j == l           = one
+    | i == l && j == k           = one
+    | i == j && i /= k && i /= l = one
+    | otherwise                  = zero
+
+affine :: (Algebra.Ring.C k) => Int -> Int -> [(Int,k)] -> Matrix k
+affine n k ps = matrixFromFunction n n $ \ i j -> δ i j + c i j where
+  m = Map.fromList ps
+  c i j = case (Map.lookup i m,j == k) of
+    (Just cᵢ,True) -> cᵢ
+    _              -> zero
+
+permute :: (Algebra.Ring.C k) => [Int] -> Matrix k
+permute is = matrixFromFunction n n (\ i j -> δ (is !! (i - 1)) j) where
+  n = length is
+
+-- Example
+x :: Matrix Rational
+x = 
+  [
+    [0, 3,-6, 6,4,5 ]
+  , [3,-7, 8,-5,8,9 ]
+  , [3,-9,12,-9,6,15]
+  ]
diff --git a/src/Algebra/Linear/Integral.hs b/src/Algebra/Linear/Integral.hs
new file mode 100644
--- /dev/null
+++ b/src/Algebra/Linear/Integral.hs
@@ -0,0 +1,108 @@
+{-# LANGUAGE NoImplicitPrelude,RebindableSyntax #-}
+module Algebra.Linear.Integral
+  (
+    divModUpperTriangular
+  , divModLowerTriangular
+  , hermite
+  ) where
+
+
+import           Algebra.Linear
+import           Auxiliary     (minimumAmong,adorn)
+
+import qualified Algebra.Ring
+import           NumericPrelude
+
+import           Control.Arrow (first,second,(***))
+import           Data.Function (on)
+import           Data.List     (partition)
+
+
+divModUpperTriangular :: Vector Integer -> Matrix Integer -> (Vector Integer,Vector Integer)
+divModUpperTriangular x u = reverse *** reverse $ divModLowerTriangular (reverse x) (reverse . map reverse $ u)
+
+divModLowerTriangular :: Vector Integer -> Matrix Integer -> (Vector Integer,Vector Integer)
+divModLowerTriangular x l = go l x [] [] where
+  go []        []        q r = (q,r)
+  go (lᵢ : l') (xᵢ : x') q r = go l' x' (q ++ [qᵢ]) (r ++ [rᵢ]) where
+    (lFirst,lᵢᵢ : _) = splitAt (length q) lᵢ
+    y = xᵢ - innerProduct q lFirst
+    (qᵢ,rᵢ)
+      | lᵢᵢ == 0  = error "Algebra.Linear.divModLowerTriangular: zero on diagonal"
+      | otherwise = y `divMod` lᵢᵢ
+
+-- Compute the Hermite normal form of the matrix,
+-- together with the unimodular basis transformation matrix.
+hermite :: Matrix Integer -> (Matrix Integer,Matrix Integer)
+hermite []          = ([],[])
+hermite xs@([] : _) = (xs,identity (length xs))
+hermite vs          = case minimumAmong (compare `on` (abs . head . fst)) nonZero of
+  Nothing                    -> first (map (0 :)) $ hermite (map tail vs)
+  Just ((v,i),[])            -> case hermite (map (tail . fst) startZero) of
+    (h,u) -> ((positive v :) . map (0 :) $ h,changeSign v `matrixProduct` rowSwap n (1,i) `matrixProduct` shift u)
+  Just ((v@(v₀ : _),i),rest) -> let
+    (reduced,translates) = unzip . flip map rest $ \ (x@(x₀ : _),j) -> let
+      (c,_) = x₀ `divMod` v₀ in ((zipWith (\ vᵢ xᵢ -> xᵢ - c * vᵢ) v x,j),(j,c))
+    (h,u) = hermite $ v : map fst reduced ++ map fst startZero
+      in (h,
+          u
+            `matrixProduct` permute (i : map snd reduced ++ map snd startZero)
+            `matrixProduct` affine n i (map (second negate) translates)
+          )
+ where
+  (startZero,nonZero) = partition ((==) 0 . head . fst) . flip zip [1 ..] $ vs
+  positive v = map (* signum (head v)) v
+  changeSign v = matrixFromFunction n n f where
+    f i j
+      | i == j && i == 1 = signum (head v)
+      | i == j           = 1
+      | otherwise        = 0
+  shift u = (1 : replicate (n - 1) 0) : map (0 :) u
+  n = length vs
+
+
+-- adjoint :: Matrix Integer -> Matrix Integer
+-- adjoint m = strip . fst . hermite . adorn $ m where
+--   n = length m
+--   strip = map (drop n)
+
+-- Alternative implementation of Hermite normal form
+{-
+hermite' :: Matrix Integer -> Matrix Integer
+hermite' []          = []
+hermite' xs@([] : _) = xs
+hermite' vs          = case minimumAmong (compare `on` (abs . head)) nonZero of
+  Nothing                -> map (0 :) $ hermite' (map tail vs)
+  Just (v,[])            -> (positive v :) . map (0 :) $ hermite' (map tail startZero)
+  Just (v@(v₀ : _),rest) -> hermite' $ v : reduced ++ startZero where
+    reduced = map (\ x@(x₀ : _) -> let (c,_) = x₀ `divMod` v₀ in zipWith (\ vᵢ xᵢ -> xᵢ - c * vᵢ) v x) rest
+ where
+  (startZero,nonZero) = partition ((==) 0 . head) vs
+  positive v = map (* signum (head v)) v
+-}
+
+-- Examples
+
+la :: Matrix Integer
+la = 
+  [
+    [ 6, 2,-4,-4, 2,-2,-2]
+  , [ 2, 6,-4,-4, 2,-2,-2]                                                                                                                                                                
+  , [-4,-4, 8, 4,-4, 0, 4]                                                                                                                                                                 
+  , [-4,-4, 4, 8,-4, 4, 0]                                                                                                                                                                 
+  , [ 2, 2,-4,-4, 6,-2,-2]                                                                                                                                                                
+  , [-2,-2, 0, 4,-2, 6,-2]                                                                                                                                                                
+  , [-2,-2, 4, 0,-2,-2, 6]
+  ]
+
+l :: Matrix Integer
+l = 
+  [
+    [2,1,1,1,1,1,1]                                                                                                                                                                    
+  , [1,2,1,1,1,1,1]                                                                                                                                                                    
+  , [1,1,2,0,1,1,0]                                                                                                                                                                    
+  , [1,1,0,2,1,0,1]                                                                                                                                                                    
+  , [1,1,1,1,2,1,1]                                                                                                                                                                   
+  , [1,1,1,0,1,2,1]                                                                                                                                                                    
+  , [1,1,0,1,1,1,2]
+  ]
diff --git a/src/Algebra/Linear/Subspace.hs b/src/Algebra/Linear/Subspace.hs
new file mode 100644
--- /dev/null
+++ b/src/Algebra/Linear/Subspace.hs
@@ -0,0 +1,122 @@
+{-# LANGUAGE NoImplicitPrelude,RebindableSyntax #-}
+{-# LANGUAGE ScopedTypeVariables,FlexibleContexts #-}
+module Algebra.Linear.Subspace
+  (
+    Subspace
+  , fromGenerators
+  , span
+  , line
+  , empty
+  , fromRelations
+  
+  , inside
+  , basis
+  , union
+  , intersection
+  , pullback
+  , image
+  , kernel
+  , dimension
+  ) where
+
+
+import           Algebra.Linear
+  (
+    Vector
+  , Matrix
+  
+  , Relation(Relation)
+  , getRelation
+  
+  , satisfies
+  , solve
+  , equations
+  , matrixProduct
+  )
+
+import qualified Algebra.Field
+import           NumericPrelude hiding (span)
+
+import           Data.List       (transpose,genericLength)
+import           Data.Proxy      (Proxy(Proxy))
+import           Data.Reflection (Reifies,reflect)
+
+
+data Subspace k
+  = Generators Integer [Vector k]  -- These are linearly independent.
+  | Relations Integer [Relation k] -- These might be redundant.
+
+degree :: Subspace k -> Integer
+degree (Generators d _) = d
+degree (Relations  d _) = d
+
+instance (Show k) => Show (Subspace k) where
+  show (Generators d gs) = "Subspace (of space of dimension " ++ show d ++ ") generated by " ++ show gs
+  show (Relations d rs)  = "Subspace (of space of dimension " ++ show d ++ ") defined by relations " ++ show (map getRelation rs)
+
+fromGenerators :: (Algebra.Field.C k,Eq k) => Integer -> [Vector k] -> Subspace k
+fromGenerators d gs = Relations d (equations d gs)
+
+span :: (Algebra.Field.C k,Eq k) => [Vector k] -> Subspace k
+span gs = fromGenerators d gs where
+  d = case gs of
+    g : _ -> genericLength g
+    []    -> error "Algebra.Linear.Subspace.fromGenerators: no generators"
+
+line :: (Algebra.Field.C k,Eq k) => Vector k -> Subspace k
+line = span . (: [])
+
+empty :: Integer -> Subspace k
+empty d = Generators d []
+
+fromRelations :: Integer -> [Relation k] -> Subspace k
+fromRelations = Relations
+
+
+inside :: (Algebra.Field.C k,Eq k) => Subspace k -> Vector k -> Bool
+inside s v = all (v `satisfies`) (toRelations s)
+
+basis :: (Algebra.Field.C k,Eq k) => Subspace k -> [Vector k]
+basis (Generators _ gs) = gs
+basis (Relations d rs) = solve d rs
+
+toRelations :: (Algebra.Field.C k,Eq k) => Subspace k -> [Relation k]
+toRelations (Generators d gs) = equations d gs
+toRelations (Relations _ rs) = rs
+
+union :: (Algebra.Field.C k,Eq k) => Subspace k -> Subspace k -> Subspace k
+union a b = fromGenerators (sameDegree a b) $ basis a ++ basis b
+
+intersection :: (Algebra.Field.C k,Eq k) => Subspace k -> Subspace k -> Subspace k
+intersection a b = Relations (sameDegree a b) $ toRelations a ++ toRelations b where
+
+sameDegree :: Subspace k -> Subspace k -> Integer
+sameDegree a b
+  | degree a == degree b = degree a
+  | otherwise            = error "Algebra.Linear.Subspace.sameDegree: subspaces of different degree"
+
+pullback :: (Algebra.Field.C k,Eq k) => Matrix k -> Subspace k -> Subspace k
+pullback [] = error "Algebra.Linear.Subspace.pullback: empty matrix"
+pullback m  =
+    Relations d
+  . map Relation
+  . transpose
+  . matrixProduct (transpose m)
+  . transpose
+  . map getRelation
+  . toRelations
+  . intersection (image m)
+ where
+  d = genericLength (head m)
+
+image :: (Algebra.Field.C k,Eq k) => Matrix k -> Subspace k
+image rows = fromGenerators d . transpose $ rows where
+  d = genericLength rows
+
+kernel :: Matrix k -> Subspace k
+kernel []   = error "Algebra.Linear.Subspace.kernel: empty matrix"
+kernel rows = Relations d . map Relation $ rows where
+  d = genericLength (head rows)
+
+dimension :: (Algebra.Field.C k,Eq k) => Subspace k -> Integer
+dimension = genericLength . basis
diff --git a/src/Algebra/Module/Free.hs b/src/Algebra/Module/Free.hs
new file mode 100644
--- /dev/null
+++ b/src/Algebra/Module/Free.hs
@@ -0,0 +1,131 @@
+{-# LANGUAGE NoImplicitPrelude #-}
+{-# LANGUAGE FlexibleInstances,MultiParamTypeClasses #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE DeriveGeneric #-}
+module Algebra.Module.Free
+  (
+    Free
+  , first
+  , second
+    
+  , Map
+  , Domain
+  , Codomain
+  , apply
+  
+  , matrix
+  , vector
+  , fromVector
+  
+  , Enumerable
+  , enumerate
+  , coefficient
+  , basisVector
+  , fromList
+  , FreeBasis(FreeBasis)
+  ) where
+
+
+import           Algebra.Linear (Matrix,Vector)
+
+import qualified Algebra.Additive
+import qualified Algebra.Module
+import qualified Algebra.ModuleBasis
+import qualified Algebra.Ring
+import NumericPrelude
+
+-- import           Data.Bifunctor (Bifunctor,first,second)
+import           Data.Binary  (Binary)
+import           Data.List    (intercalate)
+import qualified Data.Map.Strict as Map
+import           Data.Proxy   (Proxy(Proxy))
+import           GHC.Generics (Generic)
+
+
+newtype Free k t
+  = Free (Map.Map t k)
+  deriving (Generic)
+
+-- The natural Bifunctor Free instance is not possible, because of the constraints.
+first :: (k₁ -> k₂) -> Free k₁ t -> Free k₂ t
+first f (Free m) = Free (fmap f m)
+second :: (Algebra.Additive.C k,Ord t₂) => (t₁ -> t₂) -> Free k t₁ -> Free k t₂
+second g (Free m) = Free (Map.mapKeysWith (+) g m)
+
+instance (Algebra.Additive.C k,Ord t) => Algebra.Additive.C (Free k t) where
+  zero = Free Map.empty
+  negate (Free m) = Free (fmap negate m)
+  Free m₁ + Free m₂ = Free $ Map.unionWith (+) m₁ m₂
+
+instance (Algebra.Ring.C k,Ord t) => Algebra.Module.C k (Free k t) where
+  c *> (Free m) = Free (fmap (c *) m)
+
+instance (Show t,Ord t,Show k) => Show (Free k t) where
+  show (Free m) = intercalate " + " . map (\ (x,c) -> show c ++ " · " ++ show x) $ Map.toList m
+
+instance (Binary t,Binary k) => Binary (Free k t)
+
+-- 'b' is a basis for the module 'm'
+class ModuleBasis b m where
+  type Scalar b m
+  type BasisElement b m
+  coefficient :: b -> BasisElement b m -> m -> Scalar b m
+  basisVector :: b -> BasisElement b m -> m
+
+fromList :: (ModuleBasis b m,Algebra.Module.C (Scalar b m) m) => b -> [(BasisElement b m,Scalar b m)] -> m
+fromList b = sum . map (\ (x,c) -> c *> basisVector b x)
+
+data FreeBasis
+  = FreeBasis
+
+instance (Algebra.Ring.C k,Ord t) => ModuleBasis FreeBasis (Free k t) where
+  type Scalar FreeBasis (Free k t) = k
+  type BasisElement FreeBasis (Free k t) = t
+  coefficient FreeBasis x (Free m) = Map.findWithDefault zero x m
+  basisVector FreeBasis x = Free (Map.singleton x one)
+
+class (Ord t) => Enumerable p t where
+  enumerate :: p -> [t]
+
+-- instance forall k t. (Enumerable t,Algebra.Ring.C k) => Algebra.ModuleBasis.C k (Free k t) where
+--   dimension _ _ = length $ enumerate (Proxy :: Proxy t)
+--   flatten (Free m) = map (flip coefficient m) $ enumerate (Proxy :: Proxy t)
+--   basis _ = map basisVector $ enumerate (Proxy :: Proxy t)
+
+
+class Map k m where
+  type Domain m
+  type Codomain m
+  apply :: m -> Domain m -> Codomain m
+
+instance (Algebra.Ring.C k,Ord t₁,Ord t₂) => Map k (Free k (t₁,t₂)) where
+  type Domain (Free k (t₁,t₂)) = Free k t₁
+  type Codomain (Free k (t₁,t₂)) = Free k t₂
+  apply (Free m) (Free v₁) = Free $ Map.foldrWithKey
+    (\ (x₁,x₂) c -> case Map.lookup x₁ v₁ of
+      Nothing -> id
+      Just d  -> Map.insertWith (+) x₂ (c * d)
+    )
+    Map.empty
+    m
+
+matrix :: forall k t₁ t₂ p₁ p₂. (Enumerable p₁ t₁,Enumerable p₂ t₂,Algebra.Ring.C k) => p₁ -> p₂ -> Free k (t₁,t₂) -> Matrix k
+matrix p₁ p₂ m = [[coefficient FreeBasis (x₁,x₂) m | x₁ <- enumerate p₁] | x₂ <- enumerate p₂]
+
+vector :: forall k t p. (Enumerable p t,Algebra.Ring.C k) => p -> Free k t -> Vector k
+vector p v = [coefficient FreeBasis x v | x <- enumerate p]
+
+fromVector :: forall k t p. (Enumerable p t,Algebra.Ring.C k) => p -> Vector k -> Free k t
+fromVector p = Free . Map.fromList . zip (enumerate p)
+
+-- data X2 = S1 | S2 deriving (Eq,Ord)
+-- instance Enumerable X2 where
+--   enumerate _ = [S1,S2]
+-- data X3 = T1 | T2 | T3 deriving (Eq,Ord)
+-- instance Enumerable X3 where
+--   enumerate _ = [T1,T2,T3]
+-- 
+-- m :: Free Integer (X2,X3)
+-- m = (2 :: Integer) *> basisVector FreeBasis (S2,T2) - (3 :: Integer) *> basisVector FreeBasis (S1,T3) :: Free Integer (X2,X3)
diff --git a/src/Auxiliary.hs b/src/Auxiliary.hs
new file mode 100644
--- /dev/null
+++ b/src/Auxiliary.hs
@@ -0,0 +1,39 @@
+{-# LANGUAGE NoImplicitPrelude #-}
+module Auxiliary where
+
+
+import           Algebra.Additive    (zero)
+import qualified Algebra.Ring
+import           Algebra.Ring        (one)
+import           NumericPrelude
+
+import           Control.Applicative ((<$>))
+import           Control.Arrow       (second)
+
+
+adorn :: (Algebra.Ring.C k,Eq k) => [[k]] -> [[k]]
+adorn vs = map f . zip [1 .. n] $ vs where
+  f (i,v) = v ++ map (δ i) [1 .. n]
+  n = length vs
+
+δ :: (Eq a,Algebra.Ring.C k) => a -> a -> k
+δ i j
+  | i == j    = one
+  | otherwise = zero
+
+findAmong :: (a -> Bool) -> [a] -> Maybe (a,[a])
+findAmong p [] = Nothing
+findAmong p (x : xs)
+  | p x       = Just (x,xs)
+  | otherwise = second (x :) <$> findAmong p xs
+
+minimumAmong :: (a -> a -> Ordering) -> [a] -> Maybe (a,[a])
+minimumAmong _ []  = Nothing
+minimumAmong _ [x] = Just (x,[])
+minimumAmong (~~) (x : y : xs) = fmap (second (l :)) $ minimumAmong (~~) (s : xs) where
+  (s,l) 
+    | x ~~ y == LT = (x,y)
+    | otherwise    = (y,x)
+
+headE _ (x : _) = x
+headE s _       = error s