-- Copyright (c) 2008-2015, David Amos. All rights reserved.


-- |A module providing elementary operations involving scalars, vectors, and matrices

-- over a ring or field. Vectors are represented as [a], matrices as [[a]].

-- (No distinction is made between row and column vectors.)

-- It is the caller's responsibility to ensure that the lists have the correct number of elements.

--

-- The mnemonic for many of the arithmetic operations is that the number of angle brackets

-- on each side indicates the dimension of the argument on that side. For example,

-- v \<*\>\> m is multiplication of a vector on the left by a matrix on the right.

module Math.Algebra.LinearAlgebra where

import Prelude hiding ( (*>), (<*>) )

import qualified Data.List as L
import Math.Core.Field -- not actually used in this module



infixr 8 *>, *>>
infixr 7 <<*>
infixl 7 <.>, <*>, <<*>>, <*>>
infixl 6 <+>, <->, <<+>>, <<->>

-- The mnemonic for these operations is that the number of angle brackets on each side indicates the dimension of the argument on that side



-- vector operations


-- |u \<+\> v returns the sum u+v of vectors

(<+>) :: (Num a) => [a] -> [a] -> [a]
u :: [a]
u <+> :: [a] -> [a] -> [a]
<+> v :: [a]
v = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(+) [a]
u [a]
v

-- |u \<-\> v returns the difference u-v of vectors

(<->) :: (Num a) => [a] -> [a] -> [a]
u :: [a]
u <-> :: [a] -> [a] -> [a]
<-> v :: [a]
v = (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [a]
u [a]
v

-- |k *\> v returns the product k*v of the scalar k and the vector v

(*>) :: (Num a) => a -> [a] -> [a]
k :: a
k *> :: a -> [a] -> [a]
*> v :: [a]
v = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a
ka -> a -> a
forall a. Num a => a -> a -> a
*) [a]
v

-- |u \<.\> v returns the dot product of vectors (also called inner or scalar product)

(<.>) :: (Num a) => [a] -> [a] -> a
u :: [a]
u <.> :: [a] -> [a] -> a
<.> v :: [a]
v = [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
u [a]
v)

-- |u \<*\> v returns the tensor product of vectors (also called outer or matrix product)

(<*>) :: (Num a) => [a] -> [a] -> [[a]]
u :: [a]
u <*> :: [a] -> [a] -> [[a]]
<*> v :: [a]
v = [ [a
aa -> a -> a
forall a. Num a => a -> a -> a
*a
b | a
b <- [a]
v] | a
a <- [a]
u]


-- matrix operations


-- |a \<\<+\>\> b returns the sum a+b of matrices

(<<+>>) :: (Num a) => [[a]] -> [[a]] -> [[a]]
a :: [[a]]
a <<+>> :: [[a]] -> [[a]] -> [[a]]
<<+>> b :: [[a]]
b = (([a] -> [a] -> [a]) -> [[a]] -> [[a]] -> [[a]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (([a] -> [a] -> [a]) -> [[a]] -> [[a]] -> [[a]])
-> ((a -> a -> a) -> [a] -> [a] -> [a])
-> (a -> a -> a)
-> [[a]]
-> [[a]]
-> [[a]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith) a -> a -> a
forall a. Num a => a -> a -> a
(+) [[a]]
a [[a]]
b

-- |a \<\<-\>\> b returns the difference a-b of matrices

(<<->>) :: (Num a) => [[a]] -> [[a]] -> [[a]]
a :: [[a]]
a <<->> :: [[a]] -> [[a]] -> [[a]]
<<->> b :: [[a]]
b = (([a] -> [a] -> [a]) -> [[a]] -> [[a]] -> [[a]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (([a] -> [a] -> [a]) -> [[a]] -> [[a]] -> [[a]])
-> ((a -> a -> a) -> [a] -> [a] -> [a])
-> (a -> a -> a)
-> [[a]]
-> [[a]]
-> [[a]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith) (-) [[a]]
a [[a]]
b

-- |a \<\<*\>\> b returns the product a*b of matrices

(<<*>>) :: (Num a) => [[a]] -> [[a]] -> [[a]]
a :: [[a]]
a <<*>> :: [[a]] -> [[a]] -> [[a]]
<<*>> b :: [[a]]
b = [ [[a]
u [a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
<.> [a]
v | [a]
v <- [[a]] -> [[a]]
forall a. [[a]] -> [[a]]
L.transpose [[a]]
b] | [a]
u <- [[a]]
a]
 
-- |k *\>\> m returns the product k*m of the scalar k and the matrix m

(*>>) :: (Num a) => a -> [[a]] -> [[a]]
k :: a
k *>> :: a -> [[a]] -> [[a]]
*>> m :: [[a]]
m = (([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map (([a] -> [a]) -> [[a]] -> [[a]])
-> ((a -> a) -> [a] -> [a]) -> (a -> a) -> [[a]] -> [[a]]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map) (a
ka -> a -> a
forall a. Num a => a -> a -> a
*) [[a]]
m

-- |m \<\<*\> v is multiplication of a vector by a matrix on the left

(<<*>) :: (Num a) => [[a]] -> [a] -> [a]
m :: [[a]]
m <<*> :: [[a]] -> [a] -> [a]
<<*> v :: [a]
v = ([a] -> a) -> [[a]] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map ([a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
<.> [a]
v) [[a]]
m

-- |v \<*\>\> m is multiplication of a vector by a matrix on the right

(<*>>) :: (Num a) => [a] -> [[a]] -> [a]
v :: [a]
v <*>> :: [a] -> [[a]] -> [a]
<*>> m :: [[a]]
m = ([a] -> a) -> [[a]] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map ([a]
v [a] -> [a] -> a
forall a. Num a => [a] -> [a] -> a
<.>) ([[a]] -> [[a]]
forall a. [[a]] -> [[a]]
L.transpose [[a]]
m)


fMatrix :: t -> (t -> t -> a) -> [[a]]
fMatrix n :: t
n f :: t -> t -> a
f = [[t -> t -> a
f t
i t
j | t
j <- [1..t
n]] | t
i <- [1..t
n]] 

-- version with indices from zero

fMatrix' :: t -> (t -> t -> a) -> [[a]]
fMatrix' n :: t
n f :: t -> t -> a
f = [[t -> t -> a
f t
i t
j | t
j <- [0..t
nt -> t -> t
forall a. Num a => a -> a -> a
-1]] | t
i <- [0..t
nt -> t -> t
forall a. Num a => a -> a -> a
-1]] 


-- idMx n = fMatrix n (\i j -> if i == j then 1 else 0)


idMx :: Int -> [[a]]
idMx n :: Int
n = [[[a]]]
idMxs [[[a]]] -> Int -> [[a]]
forall a. [a] -> Int -> a
!! Int
n where
    idMxs :: [[[a]]]
idMxs = ((Int, [[a]]) -> [[a]]) -> [(Int, [[a]])] -> [[[a]]]
forall a b. (a -> b) -> [a] -> [b]
map (Int, [[a]]) -> [[a]]
forall a b. (a, b) -> b
snd ([(Int, [[a]])] -> [[[a]]]) -> [(Int, [[a]])] -> [[[a]]]
forall a b. (a -> b) -> a -> b
$ ((Int, [[a]]) -> (Int, [[a]])) -> (Int, [[a]]) -> [(Int, [[a]])]
forall a. (a -> a) -> a -> [a]
iterate (Int, [[a]]) -> (Int, [[a]])
forall a. Num a => (Int, [[a]]) -> (Int, [[a]])
next (0,[])
    next :: (Int, [[a]]) -> (Int, [[a]])
next (j :: Int
j,m :: [[a]]
m) = (Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
+1, (1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: Int -> a -> [a]
forall a. Int -> a -> [a]
replicate Int
j 0) [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: ([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map (0a -> [a] -> [a]
forall a. a -> [a] -> [a]
:) [[a]]
m)

-- |iMx n is the n*n identity matrix

iMx :: (Num t) => Int -> [[t]]
iMx :: Int -> [[t]]
iMx n :: Int
n = Int -> [[t]]
forall a. Num a => Int -> [[a]]
idMx Int
n

-- |jMx n is the n*n matrix of all 1s

jMx :: (Num t) => Int -> [[t]]
jMx :: Int -> [[t]]
jMx n :: Int
n = Int -> [t] -> [[t]]
forall a. Int -> a -> [a]
replicate Int
n (Int -> t -> [t]
forall a. Int -> a -> [a]
replicate Int
n 1)

-- |zMx n is the n*n matrix of all 0s

zMx :: (Num t) => Int -> [[t]]
zMx :: Int -> [[t]]
zMx n :: Int
n = Int -> [t] -> [[t]]
forall a. Int -> a -> [a]
replicate Int
n (Int -> t -> [t]
forall a. Int -> a -> [a]
replicate Int
n 0)

{-
-- VECTORS

data Vector d k = V [k] deriving (Eq,Ord,Show) 

instance (IntegerAsType d, Num k) => Num (Vector d k) where
    V a + V b = V $ a <+> b
    V a - V b = V $ a <-> b
    negate (V a) = V $ map negate a
    fromInteger 0 = V $ replicate d' 0 where d' = fromInteger $ value (undefined :: d)

V v <>> M m = V $ v <*>> m

M m <<> V v = V $ m <<*> v

k |> V v = V $ k *> v
-}

-- MATRICES


{-
-- Square matrices of dimension d over field k
data Matrix d k = M [[k]] deriving (Eq,Ord,Show)

instance (IntegerAsType d, Num k) => Num (Matrix d k) where
    M a + M b = M $ a <<+>> b
    M a - M b = M $ a <<->> b
    negate (M a) = M $ (map . map) negate a
    M a * M b = M $ a <<*>> b
    fromInteger 0 = M $ zMx d' where d' = fromInteger $ value (undefined :: d)
    fromInteger 1 = M $ idMx d' where d' = fromInteger $ value (undefined :: d)

instance (IntegerAsType d, Fractional a) => Fractional (Matrix d a) where
	recip (M a) = case inverse a of
		Nothing -> error "Matrix.recip: matrix is singular"
		Just a' -> M a'
-}


-- |The inverse of a matrix (over a field), if it exists

inverse :: (Eq a, Fractional a) => [[a]] -> Maybe [[a]]
inverse :: [[a]] -> Maybe [[a]]
inverse m :: [[a]]
m =
    let d :: Int
d = [[a]] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[a]]
m -- the dimension

        i :: [[a]]
i = Int -> [[a]]
forall a. Num a => Int -> [[a]]
idMx Int
d
        m' :: [[a]]
m' = ([a] -> [a] -> [a]) -> [[a]] -> [[a]] -> [[a]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
(++) [[a]]
m [[a]]
i
        i1 :: [[a]]
i1 = [[a]] -> [[a]]
forall a. (Eq a, Fractional a) => [[a]] -> [[a]]
inverse1 [[a]]
m'
        i2 :: [[a]]
i2 = [[a]] -> [[a]]
forall a. (Eq a, Num a) => [[a]] -> [[a]]
inverse2 [[a]]
i1
    in if [[a]] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[a]]
i1 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
d
       then [[a]] -> Maybe [[a]]
forall a. a -> Maybe a
Just [[a]]
i2
       else Maybe [[a]]
forall a. Maybe a
Nothing

-- given (M|I), use row operations to get to (U|A), where U is upper triangular with 1s on diagonal

inverse1 :: [[a]] -> [[a]]
inverse1 [] = []
inverse1 ((x :: a
x:xs :: [a]
xs):rs :: [[a]]
rs) =
    if a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= 0
    then let r' :: [a]
r' = (1a -> a -> a
forall a. Fractional a => a -> a -> a
/a
x) a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
*> [a]
xs
         in (1a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
r') [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: [[a]] -> [[a]]
inverse1 [[a]
ys [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
<-> a
y a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
*> [a]
r' | (y :: a
y:ys :: [a]
ys) <- [[a]]
rs]
    else case ([a] -> Bool) -> [[a]] -> [[a]]
forall a. (a -> Bool) -> [a] -> [a]
filter (\r' :: [a]
r' -> [a] -> a
forall a. [a] -> a
head [a]
r' a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= 0) [[a]]
rs of
         [] -> [] -- early termination, which will be detected in calling function

         r :: [a]
r:_ -> [[a]] -> [[a]]
inverse1 (((a
xa -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
xs) [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
<+> [a]
r) [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: [[a]]
rs)
-- This is basically row echelon form


-- given (U|A), use row operations to get to M^-1

inverse2 :: [[a]] -> [[a]]
inverse2 [] = []
inverse2 ((1:r :: [a]
r):rs :: [[a]]
rs) = [a] -> [[a]] -> [a]
forall a. (Eq a, Num a) => [a] -> [[a]] -> [a]
inverse2' [a]
r [[a]]
rs [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: [[a]] -> [[a]]
inverse2 [[a]]
rs where
    inverse2' :: [a] -> [[a]] -> [a]
inverse2' xs :: [a]
xs [] = [a]
xs
    inverse2' (x :: a
x:xs :: [a]
xs) ((1:r :: [a]
r):rs :: [[a]]
rs) = [a] -> [[a]] -> [a]
inverse2' ([a]
xs [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
<-> a
x a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
*> [a]
r) [[a]]
rs

xs :: [a]
xs ! :: [a] -> Int -> a
! i :: Int
i = [a]
xs [a] -> Int -> a
forall a. [a] -> Int -> a
!! (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-1) -- ie, a 1-based list lookup instead of 0-based


rowEchelonForm :: [[a]] -> [[a]]
rowEchelonForm [] = []
rowEchelonForm ((x :: a
x:xs :: [a]
xs):rs :: [[a]]
rs) =
    if a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= 0
    then let r' :: [a]
r' = (1a -> a -> a
forall a. Fractional a => a -> a -> a
/a
x) a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
*> [a]
xs
         in (1a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
r') [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: ([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map (0a -> [a] -> [a]
forall a. a -> [a] -> [a]
:) ([[a]] -> [[a]]
rowEchelonForm [[a]
ys [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
<-> a
y a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
*> [a]
r' | (y :: a
y:ys :: [a]
ys) <- [[a]]
rs])
    else case ([a] -> Bool) -> [[a]] -> [[a]]
forall a. (a -> Bool) -> [a] -> [a]
filter (\r' :: [a]
r' -> [a] -> a
forall a. [a] -> a
head [a]
r' a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= 0) [[a]]
rs of
         [] -> ([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map (0a -> [a] -> [a]
forall a. a -> [a] -> [a]
:) ([[a]] -> [[a]]
rowEchelonForm ([[a]] -> [[a]]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> a -> b
$ [a]
xs [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: ([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map [a] -> [a]
forall a. [a] -> [a]
tail [[a]]
rs)
         r :: [a]
r:_ -> [[a]] -> [[a]]
rowEchelonForm (((a
xa -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
xs) [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
<+> [a]
r) [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: [[a]]
rs)
rowEchelonForm zs :: [[a]]
zs@([]:_) = [[a]]
zs

reducedRowEchelonForm :: (Eq a, Fractional a) => [[a]] -> [[a]]
reducedRowEchelonForm :: [[a]] -> [[a]]
reducedRowEchelonForm m :: [[a]]
m = [[a]] -> [[a]]
forall a. [a] -> [a]
reverse ([[a]] -> [[a]]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> a -> b
$ [[a]] -> [[a]]
forall a. (Eq a, Num a) => [[a]] -> [[a]]
reduce ([[a]] -> [[a]]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> a -> b
$ [[a]] -> [[a]]
forall a. [a] -> [a]
reverse ([[a]] -> [[a]]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> a -> b
$ [[a]] -> [[a]]
forall a. (Eq a, Fractional a) => [[a]] -> [[a]]
rowEchelonForm [[a]]
m where
    reduce :: [[a]] -> [[a]]
reduce (r :: [a]
r:rs :: [[a]]
rs) = let r' :: [a]
r':rs' :: [[a]]
rs' = [[a]] -> [[a]]
forall a. (Eq a, Num a) => [[a]] -> [[a]]
reduceStep ([a]
r[a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
:[[a]]
rs) in [a]
r' [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: [[a]] -> [[a]]
reduce [[a]]
rs' -- is this scanl or similar?

    reduce [] = []
    reduceStep :: [[a]] -> [[a]]
reduceStep ((1:xs :: [a]
xs):rs :: [[a]]
rs) = (1a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
xs) [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: [ 0a -> [a] -> [a]
forall a. a -> [a] -> [a]
: ([a]
ys [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
<-> a
y a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
*> [a]
xs) | y :: a
y:ys :: [a]
ys <- [[a]]
rs]
    reduceStep rs :: [[a]]
rs@((0:_):_) = (a -> [a] -> [a]) -> [a] -> [[a]] -> [[a]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (:) (([a] -> a) -> [[a]] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map [a] -> a
forall a. [a] -> a
head [[a]]
rs) ([[a]] -> [[a]]
reduceStep ([[a]] -> [[a]]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> a -> b
$ ([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map [a] -> [a]
forall a. [a] -> [a]
tail [[a]]
rs)
    reduceStep rs :: [[a]]
rs = [[a]]
rs

-- Given a matrix m and (column) vector b, either find (column vector) x such that m x == b,

-- or indicate that there is none

solveLinearSystem :: [[a]] -> [a] -> Maybe [a]
solveLinearSystem m :: [[a]]
m b :: [a]
b =
    let augmented :: [[a]]
augmented = ([a] -> a -> [a]) -> [[a]] -> [a] -> [[a]]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\r :: [a]
r x :: a
x -> [a]
r [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
x]) [[a]]
m [a]
b -- augmented matrix

        trisystem :: [[a]]
trisystem = [[a]] -> [[a]]
forall a. (Eq a, Fractional a) => [[a]] -> [[a]]
inverse1 [[a]]
augmented -- upper triangular form

        solution :: [a]
solution = [a] -> [a]
forall a. [a] -> [a]
reverse ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ [[a]] -> [a]
forall a. Fractional a => [[a]] -> [a]
solveTriSystem ([[a]] -> [a]) -> [[a]] -> [a]
forall a b. (a -> b) -> a -> b
$ [[a]] -> [[a]]
forall a. [a] -> [a]
reverse ([[a]] -> [[a]]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> a -> b
$ ([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map [a] -> [a]
forall a. [a] -> [a]
reverse [[a]]
trisystem
    in if [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
solution Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
b then [a] -> Maybe [a]
forall a. a -> Maybe a
Just [a]
solution else Maybe [a]
forall a. Maybe a
Nothing
    where solveTriSystem :: [[a]] -> [a]
solveTriSystem ([v :: a
v,c :: a
c]:rs :: [[a]]
rs) =
              let x :: a
x = a
va -> a -> a
forall a. Fractional a => a -> a -> a
/a
c -- the first row tells us that cx == v

                  rs' :: [[a]]
rs' = ([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map (\(v' :: a
v':c' :: a
c':r :: [a]
r) -> (a
v'a -> a -> a
forall a. Num a => a -> a -> a
-a
c'a -> a -> a
forall a. Num a => a -> a -> a
*a
x)a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
r) [[a]]
rs
              in a
x a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [[a]] -> [a]
solveTriSystem [[a]]
rs'
          solveTriSystem [] = []
          solveTriSystem _ = [] -- abnormal termination - m wasn't invertible



isZero :: t a -> Bool
isZero v :: t a
v = (a -> Bool) -> t a -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (a -> a -> Bool
forall a. Eq a => a -> a -> Bool
==0) t a
v

-- inSpanRE m v returns whether the vector v is in the span of the matrix m, where m is required to be in row echelon form

inSpanRE :: [[a]] -> [a] -> Bool
inSpanRE ((1:xs :: [a]
xs):bs :: [[a]]
bs) (y :: a
y:ys :: [a]
ys) = [[a]] -> [a] -> Bool
inSpanRE (([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map [a] -> [a]
forall a. [a] -> [a]
tail [[a]]
bs) (if a
y a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== 0 then [a]
ys else [a]
ys [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
<-> a
y a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
*> [a]
xs)
inSpanRE ((0:xs :: [a]
xs):bs :: [[a]]
bs) (y :: a
y:ys :: [a]
ys) = if a
y a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== 0 then [[a]] -> [a] -> Bool
inSpanRE ([a]
xs [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: ([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map [a] -> [a]
forall a. [a] -> [a]
tail [[a]]
bs) [a]
ys else Bool
False
inSpanRE _ ys :: [a]
ys = [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a, Num a) => t a -> Bool
isZero [a]
ys

rank :: [[a]] -> Int
rank m :: [[a]]
m = [[a]] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([[a]] -> Int) -> [[a]] -> Int
forall a b. (a -> b) -> a -> b
$ ([a] -> Bool) -> [[a]] -> [[a]]
forall a. (a -> Bool) -> [a] -> [a]
filter (Bool -> Bool
not (Bool -> Bool) -> ([a] -> Bool) -> [a] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a, Num a) => t a -> Bool
isZero) ([[a]] -> [[a]]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> a -> b
$ [[a]] -> [[a]]
forall a. (Eq a, Fractional a) => [[a]] -> [[a]]
rowEchelonForm [[a]]
m

-- kernel of a matrix

-- returns basis for vectors v s.t m <<*> v == 0

kernel :: [[a]] -> [[a]]
kernel m :: [[a]]
m = [[a]] -> [[a]]
forall a. (Ord a, Num a) => [[a]] -> [[a]]
kernelRRE ([[a]] -> [[a]]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> a -> b
$ [[a]] -> [[a]]
forall a. (Eq a, Fractional a) => [[a]] -> [[a]]
reducedRowEchelonForm [[a]]
m

kernelRRE :: [[a]] -> [[a]]
kernelRRE m :: [[a]]
m =
    let nc :: Int
nc = [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ([a] -> Int) -> [a] -> Int
forall a b. (a -> b) -> a -> b
$ [[a]] -> [a]
forall a. [a] -> a
head [[a]]
m -- the number of columns

        is :: [Int]
is = Int -> [[a]] -> [Int]
forall a a. (Eq a, Num a, Num a) => a -> [[a]] -> [a]
findLeadingCols 1 ([[a]] -> [[a]]
forall a. [[a]] -> [[a]]
L.transpose [[a]]
m) -- these are the indices of the columns which have a leading 1

        js :: [Int]
js = [1..Int
nc] [Int] -> [Int] -> [Int]
forall a. Eq a => [a] -> [a] -> [a]
L.\\ [Int]
is
        freeCols :: [(Int, [a])]
freeCols = let m' :: [[a]]
m' = Int -> [[a]] -> [[a]]
forall a. Int -> [a] -> [a]
take ([Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
is) [[a]]
m -- discard zero rows

                   in [Int] -> [[a]] -> [(Int, [a])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int]
is ([[a]] -> [(Int, [a])]) -> [[a]] -> [(Int, [a])]
forall a b. (a -> b) -> a -> b
$ [[a]] -> [[a]]
forall a. [[a]] -> [[a]]
L.transpose [([a] -> a) -> [[a]] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a
forall a. Num a => a -> a
negate (a -> a) -> ([a] -> a) -> [a] -> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([a] -> Int -> a
forall a. [a] -> Int -> a
!Int
j)) [[a]]
m' | Int
j <- [Int]
js]
        boundCols :: [(Int, [a])]
boundCols = [Int] -> [[a]] -> [(Int, [a])]
forall a b. [a] -> [b] -> [(a, b)]
zip [Int]
js (Int -> [[a]]
forall a. Num a => Int -> [[a]]
idMx (Int -> [[a]]) -> Int -> [[a]]
forall a b. (a -> b) -> a -> b
$ [Int] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Int]
js)
    in [[a]] -> [[a]]
forall a. [[a]] -> [[a]]
L.transpose ([[a]] -> [[a]]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> a -> b
$ ((Int, [a]) -> [a]) -> [(Int, [a])] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map (Int, [a]) -> [a]
forall a b. (a, b) -> b
snd ([(Int, [a])] -> [[a]]) -> [(Int, [a])] -> [[a]]
forall a b. (a -> b) -> a -> b
$ [(Int, [a])] -> [(Int, [a])]
forall a. Ord a => [a] -> [a]
L.sort ([(Int, [a])] -> [(Int, [a])]) -> [(Int, [a])] -> [(Int, [a])]
forall a b. (a -> b) -> a -> b
$ [(Int, [a])]
freeCols [(Int, [a])] -> [(Int, [a])] -> [(Int, [a])]
forall a. [a] -> [a] -> [a]
++ [(Int, [a])]
boundCols
    where
    findLeadingCols :: a -> [[a]] -> [a]
findLeadingCols i :: a
i (c :: [a]
c@(1:_):cs :: [[a]]
cs) = a
i a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a -> [[a]] -> [a]
findLeadingCols (a
ia -> a -> a
forall a. Num a => a -> a -> a
+1) (([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map [a] -> [a]
forall a. [a] -> [a]
tail [[a]]
cs)
    findLeadingCols i :: a
i (c :: [a]
c@(0:_):cs :: [[a]]
cs) = a -> [[a]] -> [a]
findLeadingCols (a
ia -> a -> a
forall a. Num a => a -> a -> a
+1) [[a]]
cs
    findLeadingCols _ _ = []

-- m ^- n = recip m ^ n


-- t (M m) = M (L.transpose m)


-- |The determinant of a matrix (over a field)

det :: (Eq a, Fractional a) => [[a]] -> a
det :: [[a]] -> a
det [[x :: a
x]] = a
x
det ((x :: a
x:xs :: [a]
xs):rs :: [[a]]
rs) =
    if a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= 0
    then let r' :: [a]
r' = (1a -> a -> a
forall a. Fractional a => a -> a -> a
/a
x) a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
*> [a]
xs
         in a
x a -> a -> a
forall a. Num a => a -> a -> a
* [[a]] -> a
forall a. (Eq a, Fractional a) => [[a]] -> a
det [[a]
ys [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
<-> a
y a -> [a] -> [a]
forall a. Num a => a -> [a] -> [a]
*> [a]
r' | (y :: a
y:ys :: [a]
ys) <- [[a]]
rs]
    else case ([a] -> Bool) -> [[a]] -> [[a]]
forall a. (a -> Bool) -> [a] -> [a]
filter (\r' :: [a]
r' -> [a] -> a
forall a. [a] -> a
head [a]
r' a -> a -> Bool
forall a. Eq a => a -> a -> Bool
/= 0) [[a]]
rs of
         [] -> 0
         r :: [a]
r:_ -> [[a]] -> a
forall a. (Eq a, Fractional a) => [[a]] -> a
det (((a
xa -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
xs) [a] -> [a] -> [a]
forall a. Num a => [a] -> [a] -> [a]
<+> [a]
r) [a] -> [[a]] -> [[a]]
forall a. a -> [a] -> [a]
: [[a]]
rs)


{-
class IntegerAsType a where
    value :: a -> Integer

data Z
instance IntegerAsType Z where
    value _ = 0

data S a
instance IntegerAsType a => IntegerAsType (S a) where
    value _ = value (undefined :: a) + 1
-}