    -- 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
    -}