packages feed

np-linear (empty) → 0.1.1.1

raw patch · 7 files changed

+660/−0 lines, 7 filesdep +basedep +binarydep +containerssetup-changed

Dependencies added: base, binary, containers, numeric-prelude, reflection, tagged

Files

+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ np-linear.cabal view
@@ -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
+ src/Algebra/Linear.hs view
@@ -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]+  ]
+ src/Algebra/Linear/Integral.hs view
@@ -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]+  ]
+ src/Algebra/Linear/Subspace.hs view
@@ -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
+ src/Algebra/Module/Free.hs view
@@ -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)
+ src/Auxiliary.hs view
@@ -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