module Math.Algebra.LinearAlgebra where
import qualified Data.List as L
import Math.Algebra.Field.Base
infixr 8 *>, *>>
infixr 7 <<*>
infixl 7 <.>, <*>, <<*>>, <*>>
infixl 6 <+>, <->, <<+>>, <<->>
u <+> v = zipWith (+) u v
u <-> v = zipWith () u v
k *> v = map (k*) v
k *>> m = (map . map) (k*) m
u <.> v = sum (zipWith (*) u v)
u <*> v = [ [a*b | b <- v] | a <- u]
a <<+>> b = (zipWith . zipWith) (+) a b
a <<->> b = (zipWith . zipWith) () a b
a <<*>> b = [ [u <.> v | v <- L.transpose b] | u <- a]
m <<*> v = map (<.> v) m
v <*>> m = map (v <.>) (L.transpose m)
fMatrix n f = [[f i j | j <- [1..n]] | i <- [1..n]]
fMatrix' n f = [[f i j | j <- [0..n1]] | i <- [0..n1]]
idMx n = idMxs !! n where
idMxs = map snd $ iterate next (0,[])
next (j,m) = (j+1, (1 : replicate j 0) : map (0:) m)
jMx n = replicate n (replicate n 1)
zMx n = replicate n (replicate n 0)
inverse m =
let d = length m
i = idMx d
m' = zipWith (++) m i
i1 = inverse1 m'
i2 = inverse2 i1
in if length i1 == d
then Just i2
else Nothing
inverse1 [] = []
inverse1 ((x:xs):rs) =
if x /= 0
then let r' = (1/x) *> xs
in (1:r') : inverse1 [ys <-> y *> r' | (y:ys) <- rs]
else case filter (\r' -> head r' /= 0) rs of
[] -> []
r:_ -> inverse1 (((x:xs) <+> r) : rs)
inverse2 [] = []
inverse2 ((1:r):rs) = inverse2' r rs : inverse2 rs where
inverse2' xs [] = xs
inverse2' (x:xs) ((1:r):rs) = inverse2' (xs <-> x *> r) rs
xs ! i = xs !! (i1)
rowEchelonForm [] = []
rowEchelonForm ((x:xs):rs) =
if x /= 0
then let r' = (1/x) *> xs
in (1:r') : map (0:) (rowEchelonForm [ys <-> y *> r' | (y:ys) <- rs])
else case filter (\r' -> head r' /= 0) rs of
[] -> map (0:) (rowEchelonForm $ xs : map tail rs)
r:_ -> rowEchelonForm (((x:xs) <+> r) : rs)
rowEchelonForm zs@([]:_) = zs
reducedRowEchelonForm m = reverse $ reduce $ reverse $ rowEchelonForm m where
reduce (r:rs) = let r':rs' = reduceStep (r:rs) in r' : reduce rs'
reduce [] = []
reduceStep ((1:xs):rs) = (1:xs) : [ 0: (ys <-> y *> xs) | y:ys <- rs]
reduceStep rs@((0:_):_) = zipWith (:) (map head rs) (reduceStep $ map tail rs)
reduceStep rs = rs
kernel m = kernelRRE $ reducedRowEchelonForm m
kernelRRE m =
let nc = length $ head m
is = findLeadingCols 1 (L.transpose m)
js = [1..nc] L.\\ is
freeCols = let m' = take (length is) m
in zip is $ L.transpose [map (negate . (!j)) m' | j <- js]
boundCols = zip js (idMx $ length js)
in L.transpose $ map snd $ L.sort $ freeCols ++ boundCols
where
findLeadingCols i (c@(1:_):cs) = i : findLeadingCols (i+1) (map tail cs)
findLeadingCols i (c@(0:_):cs) = findLeadingCols (i+1) cs
findLeadingCols _ _ = []
m ^- n = recip m ^ n
det [[x]] = x
det ((x:xs):rs) =
if x /= 0
then let r' = (1/x) *> xs
in x * det [ys <-> y *> r' | (y:ys) <- rs]
else case filter (\r' -> head r' /= 0) rs of
[] -> 0
r:_ -> det (((x:xs) <+> r) : rs)