+-- | Returns True if the given matrix is upper-triangular, and False
+-- otherwise. The parameter @epsilon@ lets the caller choose a
+-- tolerance.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,1],[1e-12,1]] :: Mat2 Double
+-- >>> is_upper_triangular m
+-- False
+-- >>> is_upper_triangular' 1e-10 m
+-- True
+--
+-- TODO:
+--
+-- 1. Don't cheat with lists.
+--
+is_upper_triangular' :: (Ord a, Ring.C a, Absolute.C a, Arity m, Arity n)
+ => a -- ^ The tolerance @epsilon@.
+ -> Mat m n a
+ -> Bool
+is_upper_triangular' epsilon m =
+ and $ concat results
+ where
+ results = [[ test i j | i <- [0..(nrows m)-1]] | j <- [0..(ncols m)-1] ]
+
+ test :: Int -> Int -> Bool
+ test i j
+ | i <= j = True
+ -- use "less than or equal to" so zero is a valid epsilon
+ | otherwise = abs (m !!! (i,j)) <= epsilon
+
+
+-- | Returns True if the given matrix is upper-triangular, and False
+-- otherwise. A specialized version of 'is_upper_triangular\'' with
+-- @epsilon = 0@.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
+-- >>> is_upper_triangular m
+-- False
+--
+-- >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
+-- >>> is_upper_triangular m
+-- True
+--
+-- TODO:
+--
+-- 1. The Ord constraint is too strong here, Eq would suffice.
+--
+is_upper_triangular :: (Ord a, Ring.C a, Absolute.C a, Arity m, Arity n)
+ => Mat m n a -> Bool
+is_upper_triangular = is_upper_triangular' 0
+
+
+-- | Returns True if the given matrix is lower-triangular, and False
+-- otherwise. This is a specialized version of 'is_lower_triangular\''
+-- with @epsilon = 0@.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
+-- >>> is_lower_triangular m
+-- True
+--
+-- >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
+-- >>> is_lower_triangular m
+-- False
+--
+is_lower_triangular :: (Ord a,
+ Ring.C a,
+ Absolute.C a,
+ Arity m,
+ Arity n)
+ => Mat m n a
+ -> Bool
+is_lower_triangular = is_upper_triangular . transpose
+
+
+-- | Returns True if the given matrix is lower-triangular, and False
+-- otherwise. The parameter @epsilon@ lets the caller choose a
+-- tolerance.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,1e-12],[1,1]] :: Mat2 Double
+-- >>> is_lower_triangular m
+-- False
+-- >>> is_lower_triangular' 1e-12 m
+-- True
+--
+is_lower_triangular' :: (Ord a,
+ Ring.C a,
+ Absolute.C a,
+ Arity m,
+ Arity n)
+ => a -- ^ The tolerance @epsilon@.
+ -> Mat m n a
+ -> Bool
+is_lower_triangular' epsilon = (is_upper_triangular' epsilon) . transpose
+
+
+-- | Returns True if the given matrix is triangular, and False
+-- otherwise.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,0],[1,1]] :: Mat2 Int
+-- >>> is_triangular m
+-- True
+--
+-- >>> let m = fromList [[1,2],[0,3]] :: Mat2 Int
+-- >>> is_triangular m
+-- True
+--
+-- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
+-- >>> is_triangular m
+-- False
+--
+is_triangular :: (Ord a,
+ Ring.C a,
+ Absolute.C a,
+ Arity m,
+ Arity n)
+ => Mat m n a
+ -> Bool
+is_triangular m = is_upper_triangular m || is_lower_triangular m
+
+
+-- | Return the (i,j)th minor of m.
+--
+-- Examples:
+--
+-- >>> let m = fromList [[1,2,3],[4,5,6],[7,8,9]] :: Mat3 Int
+-- >>> minor m 0 0 :: Mat2 Int
+-- ((5,6),(8,9))
+-- >>> minor m 1 1 :: Mat2 Int
+-- ((1,3),(7,9))
+--
+minor :: (m ~ S r,
+ n ~ S t,
+ Arity r,
+ Arity t)
+ => Mat m n a
+ -> Int
+ -> Int
+ -> Mat r t a
+minor (Mat rows) i j = m
+ where
+ rows' = delete rows i
+ m = Mat $ V.map ((flip delete) j) rows'
+
+
+class (Eq a, Ring.C a) => Determined p a where
+ determinant :: (p a) -> a
+
+instance (Eq a, Ring.C a) => Determined (Mat (S Z) (S Z)) a where
+ determinant (Mat rows) = (V.head . V.head) rows
+
+instance (Ord a,
+ Ring.C a,
+ Absolute.C a,
+ Arity n,
+ Determined (Mat (S n) (S n)) a)
+ => Determined (Mat (S (S n)) (S (S n))) a where
+ -- | The recursive definition with a special-case for triangular matrices.
+ --
+ -- Examples:
+ --
+ -- >>> let m = fromList [[1,2],[3,4]] :: Mat2 Int
+ -- >>> determinant m
+ -- -1
+ --
+ determinant m
+ | is_triangular m = product [ m !!! (i,i) | i <- [0..(nrows m)-1] ]
+ | otherwise = determinant_recursive
+ where
+ m' i j = m !!! (i,j)
+
+ det_minor i j = determinant (minor m i j)
+
+ determinant_recursive =
+ sum [ (-1)^(toInteger j) NP.* (m' 0 j) NP.* (det_minor 0 j)
+ | j <- [0..(ncols m)-1] ]
+
+
+
+-- | Matrix multiplication.