```{-# LANGUAGE TypeSynonymInstances #-}
-- | A small simple matrix library.
module Algebra.Matrix
( Vector(Vec)
, unVec, lengthVec
, Matrix(M), matrix
, matrixToVector, vectorToMatrix, unMVec, unM
, identity, mulM, addM, transpose, isSquareMatrix, dimension
) where

import qualified Data.List as L
import Test.QuickCheck

import Algebra.Structures.IntegralDomain

-------------------------------------------------------------------------------
-- | Row vectors

newtype Vector r = Vec [r] deriving (Eq)

instance Show r => Show (Vector r) where
show (Vec vs) = show vs

instance (Ring r, Arbitrary r, Eq r) => Arbitrary (Vector r) where
arbitrary = do n <- choose (1,10) :: Gen Int
liftM Vec \$ gen n
where
gen 0 = return []
gen n = do x <- arbitrary
xs <- gen (n-1)
if x == zero then return (one:xs) else return (x:xs)

{-
instance Ring r => Ring (Vector r) where
(Vec xs) <+> (Vec ys) | length xs == length ys = Vec (zipWith (<+>) xs ys)
| otherwise = error "Bad dimensions in vector addition"
(Vec xs) <*> (Vec ys) | length xs == length ys = Vec (zipWith (<*>) xs ys)
| otherwise = error "Bad dimensions in vector multiplication"
-- In order to do these we need to know the length of the vector in advance...
-- Give me dependent types!
one  = ?
zero = ?
-}

unVec :: Vector r -> [r]
unVec (Vec vs) = vs

lengthVec :: Vector r -> Int
lengthVec = length . unVec

-------------------------------------------------------------------------------
-- | Matrices

newtype Matrix r = M [Vector r]
deriving (Eq)

instance Show r => Show (Matrix r) where
show xs = case unlines (map show (unMVec xs)) of
[] -> "[]"
xs -> init xs ++ "\n"

instance (Eq r, Arbitrary r, Ring r) => Arbitrary (Matrix r) where
arbitrary = do n <- choose (1,10) :: Gen Int
m <- choose (1,10) :: Gen Int
xs <- sequence [ liftM Vec (gen n) | _ <- [1..m]]
return (M xs)
where
gen 0 = return []
gen n = do x <- arbitrary
xs <- gen (n-1)
if x == zero then return (one:xs) else return (x:xs)

-- | Construct a mxn matrix.
matrix :: [[r]] -> Matrix r
matrix xs =
let m = fromIntegral \$ length xs
n = fromIntegral \$ length (head xs)
in if length (filter (\x -> fromIntegral (length x) == n) xs) == length xs
then M (map Vec xs)
else error "matrix: Bad dimensions"

unM :: Matrix r -> [Vector r]
unM (M xs) = xs

unMVec :: Matrix r -> [[r]]
unMVec = map unVec . unM

vectorToMatrix :: Vector r -> Matrix r
vectorToMatrix = matrix . (:[]) . unVec

matrixToVector :: Matrix r -> Vector r
matrixToVector m | fst (dimension m) == 1 = head (unM m)
| otherwise              = error "matrixToVector: Bad dimension"

-- | Compute the dimension of a matrix.
dimension :: Matrix r -> (Int, Int)
dimension (M xs) | null xs   = (0,0)
| otherwise = (length xs, length (unVec (head xs)))

isSquareMatrix :: Matrix r -> Bool
isSquareMatrix (M xs) = all (== length xs) (map lengthVec xs)

-- | Transpose a matrix.
transpose :: Matrix r -> Matrix r
transpose (M xs) = matrix (L.transpose (map unVec xs))

-- | Matrix addition.
addM :: Ring r => Matrix r -> Matrix r -> Matrix r
addM (M xs) (M ys)
| dimension (M xs) == dimension (M ys) = m
| otherwise = error "Bad dimensions in matrix addition"
where
m = matrix (zipWith (zipWith (<+>)) (map unVec xs) (map unVec ys))

-- | Matrix multiplication.
mulM :: Ring r => Matrix r -> Matrix r -> Matrix r
mulM (M xs) (M ys)
| snd (dimension (M xs)) == fst (dimension (M ys)) = m
| otherwise = error "Bad dimensions in matrix multiplication"
where
m = matrix [ [ foldr1 (<+>) (zipWith (<*>) x y)
| y <- L.transpose (map unVec ys) ]
| x <- map unVec xs ]

{-
-- In order to do this the size of the matrix need to be encoded in the type
-- There is also a problem with the fact that it is not possible to add or
-- multiply matrices with bad dimensions, so the generation of matrices has to be better...
instance Ring r => Ring (Matrix r) where
(<*>) = mul
neg (Vec xs d) = Vec [ map neg x | x <- xs ] d
zero  = undefined
-}

-- | Construct a nxn identity matrix.
identity :: IntegralDomain r => Int -> Matrix r
identity n = matrix (xs 0)
where
xs x | x == n    = []
| otherwise = (replicate x zero ++ [one] ++
replicate (n-x-1) zero) : xs (x+1)

-- Specification of identity.
propLeftIdentity :: (IntegralDomain r, Eq r) => Matrix r -> Bool
propLeftIdentity a = a == identity n `mulM` a
where n = fst (dimension a)

propRightIdentity :: (IntegralDomain r, Eq r) => Matrix r -> Bool
propRightIdentity a = a == a `mulM` identity m
where m = snd (dimension a)
```