```\section{RSAGL.Matrix}

\begin{code}
module RSAGL.Matrix
(Matrix,
matrix,
rowMajorForm,
colMajorForm,
rowAt,
matrixAt,
identityMatrix,
translationMatrix,
rotationMatrix,
scaleMatrix,
xyzMatrix,
matrixMultiply,
matrixTranspose,
matrixInverse,
determinant,
matrixInversePrim)
where

import Data.List
import RSAGL.Angle
import RSAGL.Vector
import RSAGL.Auxiliary
\end{code}

The Matrix data structure stores copies of the matrix in both row-major and column-major form,
as well a copy of the Matrix's inverse.  Caching this information has been shown to lead to performance
increases.

In row-major form, each list represents one row; columns run between lists.
In column-major form, each list represents one column; rows run between lists.

In row-major form, the Haskell list representation of a matrix has the same orientation
as if written on paper.

\begin{code}
data Matrix = Matrix { rows, cols :: Int,
row_major :: [[Double]],
matrix_inverse :: Matrix,
matrix_determinant :: Double }

instance Eq Matrix where
x == y = rows x == rows y &&
cols x == cols y &&
row_major x == row_major y

instance Show Matrix where
show m = show \$ row_major m

rowMajorForm :: Matrix -> [[Double]]
rowMajorForm mat = row_major mat

colMajorForm :: Matrix -> [[Double]]
colMajorForm = transpose . row_major
\end{code}

rowAt answers the nth row of a matrix.

\begin{code}
rowAt :: Matrix -> Int -> [Double]
rowAt m n = (rowMajorForm m) !! n
\end{code}

matrixAt answers the (i'th,j'th) element of a matrix.

\begin{code}
matrixAt :: Matrix -> (Int,Int) -> Double
matrixAt m (i,j) = ((rowMajorForm m) !! j) !! i
\end{code}

\subsection{Constructing matrices}

matrix constructs a matrix from row major list form.  (Such a list form can be formatted correctly in
monospaced font and haskell syntax, so that it looks like a matrix as it would be normally written.)

\begin{code}
matrix :: [[Double]] -> Matrix
matrix dats = if all (== row_length) row_lengths
then m
else error "row lengths do not match"
where m = Matrix { rows=length dats, --the number of rows is the length of each column
cols=length \$ head dats, --the number of columns is the length of each row
row_major=dats,
matrix_inverse = (matrixInversePrim m) { matrix_inverse = m },
matrix_determinant = determinantPrim m }
row_lengths = map length dats
\end{code}

identityMatrix constructs the n by n identity matrix

\begin{code}
identityMatrix :: (Integral i) => i -> Matrix
identityMatrix n = matrix \$ map (\x -> genericReplicate x 0 ++ [1] ++ genericReplicate (n-1-x) 0) [0..n-1]
\end{code}

\begin{code}
translationMatrix :: Vector3D -> Matrix
translationMatrix (Vector3D x y z) = matrix [[1,0,0,x],
[0,1,0,y],
[0,0,1,z],
[0,0,0,1]]
\end{code}

\begin{code}
rotationMatrix :: Vector3D -> Angle -> Matrix
rotationMatrix vector angle = let s = sine angle
c = cosine angle
c' = 1 - c
(Vector3D x y z) = vectorNormalize vector
in matrix [[c+c'*x*x,     c'*y*x - s*z,   c'*z*x + s*y, 0],
[c'*x*y+s*z,   c+c'*y*y,       c'*z*y - s*x, 0],
[c'*x*z-s*y,   c'*y*z+s*x,     c+c'*z*z,     0],
[0,            0,              0,            1]]
\end{code}

\begin{code}
scaleMatrix :: Vector3D -> Matrix
scaleMatrix (Vector3D x y z) = matrix [[x, 0, 0, 0],
[0, y, 0, 0],
[0, 0, z, 0],
[0, 0, 0, 1]]
\end{code}

\texttt{xyzMatrix} constructs the matrix in which the x y and z axis are transformed to point in the direction of the specified
vectors.

\begin{code}
xyzMatrix :: Vector3D -> Vector3D -> Vector3D -> Matrix
xyzMatrix (Vector3D x1 y1 z1) (Vector3D x2 y2 z2) (Vector3D x3 y3 z3) =
matrix [[x1,x2,x3,0],
[y1,y2,y3,0],
[z1,z2,z3,0],
[0,0,0,1.0]]
\end{code}

\subsection{Matrix arithmetic}

\begin{code}
matrixAdd :: Matrix -> Matrix -> Matrix
matrixAdd m n = let new_row_major = (if and [rows m == rows n,cols m == cols n]
then map ((map (uncurry (+))).(uncurry zip)) \$ zip (row_major m) (row_major n)
in matrix new_row_major
\end{code}

\label{howtoUseMatrixMultiplyImpl}

\begin{code}
matrixMultiply :: Matrix -> Matrix -> Matrix
matrixMultiply m n | cols m /= rows n = error "matrixMultiply: dimension mismatch"
matrixMultiply m n = matrix \$ matrixMultiplyImpl (+) 0 (*) m_data n_data
where m_data = row_major m
n_data = transpose \$ row_major n

matrixTranspose :: Matrix -> Matrix
matrixTranspose = matrix . colMajorForm -- works because matrix expects row major form
\end{code}

\subsection{Inverse and determinant of a matrix}

\begin{code}
matrixInverse :: Matrix -> Matrix
matrixInverse = matrix_inverse
\end{code}

\begin{code}
determinant :: Matrix -> Double
determinant = matrix_determinant
\end{code}

reduceMatrix eliminates the row and column corresponding to a specific element of a matrix.

\begin{code}
reduceMatrix :: Matrix -> (Int,Int) -> Matrix
reduceMatrix m (i,j) =
let (above,below) = splitAt j \$ rowMajorForm m
(left,right) = splitAt i \$ transpose \$ above ++ tail below
in matrix \$ transpose \$ left ++ tail right
\end{code}

The minor of an element of a matrix is the determinant of the matrix that is formed by removing
the row and column corresponding to that element (see reduceMatrix).

\begin{code}
matrixMinor :: Matrix -> (Int,Int) -> Double
matrixMinor m ij = determinant \$ reduceMatrix m ij
\end{code}

The cofactor of m at (i,j) is the minor of m at (i,j), multiplied by -1 at checkerboarded elements.

\begin{code}
matrixCofactor :: Matrix -> (Int,Int) -> Double
matrixCofactor m (0,0) | rows m == 1 && cols m == 1 = matrixAt m (0,0)
matrixCofactor m (i,j) = (-1)^(i+j) * matrixMinor m (i,j)
\end{code}

Implementation of determinant and matrix inverse for matrices of Rationals.

\begin{code}
matrixInversePrim :: Matrix -> Matrix
matrixInversePrim m | rows m /= cols m = error "matrixInversePrim: not a square matrix"
matrixInversePrim m | determinant m == 0 = error "matrixInversePrim: det m = 0"
matrixInversePrim m =
let scale_factor = 1 / determinant m
in matrixTranspose \$ matrix [[scale_factor * matrixCofactor m (i,j) | i <- [0..(cols m-1)]]
| j <- [0..(rows m-1)]]

determinantPrim :: Matrix -> Double
determinantPrim m | rows m /= cols m = error "determinantPrim: not a square matrix"
determinantPrim (Matrix { row_major=[[x]] }) = x
determinantPrim m = sum \$ zipWith (*) (rowAt m 0) \$ map (\x -> matrixCofactor m (x,0)) [0..(cols m - 1)]
\end{code}
```