Then an n x m matrix M just needs to also carry a list of units-of-measure T of length n, and another such list U of length m, such that the unit-of-measure of M[i][j] is T[i]U[j]. Matrix multiplication goes like

Mat T U ⇒ Mat U

^{-1}S ⇒ Mat T S

and inverses go

Mat T U ⇒ Mat U

^{-1}T

^{-1}

(where -

^{-1}is pointwise inversion of every unit-of-measure in the list) and the determinant of a square matrix Mat T U has type (ΠT)(ΠU) --- that is, the product of all the units of measure from both lists.

I started writing out Gauss-Jordan in sml to see how you'd actually type algorithms, and I think it'd be pretty cute:

(* Vec [] = unit Vec (u :: U) = u * Vec U Bivec U1 U2 = Vec U1 * Vec U2 Bimat [] U1 U2 = unit Bimat (t :: T) U1 U2 = Bivec (t . U1) (t . U2) * Bimat T U1 U2 t . [] = [] t . (u :: U) = (t . u) :: (t . U) Mat [] U = unit Mat (t :: T) U = Vec (t . U) * Mat T U *) (* aug_row : {n m : Nat} -> Vec T -> Vec 1s *) fun aug_row n m [] = [] | aug_row n m (h::tl) = (if n = m then 1.0 else 0.0) :: aug_row n (m+1) tl; (* aug : {n : Nat} {Z : Vec n *} -> Mat T U -> Bimat T U (Z @ (1 / T)) *) fun aug' n [] = [] | aug' n (h::tl) = (h, aug_row n 0 h) :: aug' (n+1) tl; (* aug : Mat T U -> Bimat T U (1 / T) *) fun aug m = aug' 0 m (* smul : t -> Vec T -> Vec (t * T) *) fun smul s [] = [] | smul (s:real) (h::tl) = s * h :: smul s tl (* vplus : Vec T -> Vec T -> Vec T *) fun vplus [] [] = [] | vplus ((h1:real)::tl1) (h2::tl2) = h1 + h2 :: (vplus tl1 tl2) | vplus _ _ = raise Match (* normalize : {n:Nat} -> Bivec T1 T2 -> Bivec (T1 / T1[n]) (T2 / T1[n]) *) fun normalize n (v1, v2) = let val s = 1.0 / List.nth(v1, n) in (smul s v1, smul s v2) end (* tri_pass : {n:Nat} -> Bivec _1 _2 -> Bimat _3 _4 _5 *) fun tri_pass _ _ [] = [] | tri_pass m (v1, v2) ((w1, w2)::tl) = let val s = ~(List.nth(w1, m)) in (vplus w1 (smul s v1), vplus w2 (smul s v2)) :: tri_pass m (v1, v2) tl end fun tri m (h::tl) = let val nv = normalize m h in nv :: tri (m+1) (tri_pass m nv tl) end | tri m [] = [] val x = aug_row 0 0 [1.0, 2.0]; val x = aug [[1.0, 2.0, 3.0], [5.0, 6.0, 7.0], [8.0, 9.0, 5.0]]; val z = tri 0 x;