;;
;; Matrices
;;

(define $M.*
  (cambda $ms
    (foldl M.b.* (car ms) (cdr ms))))

(define $M.b.*
  (lambda [%m1 %m2]
    (with-symbols {j}
      (. m1~#~j m2_j))))

(define $M.*'
  (cambda $ms
    (foldl M.b.*' (car ms) (cdr ms))))

(define $M.b.*'
  (lambda [%m1 %m2]
    (with-symbols {j}
      (.' m1~#~j m2_j))))

(define $M.power
  (lambda [%m $n]
    (repeated-squaring M.* m n)))
                       
(define $M.comm
  (lambda [%m1 %m2]
    (with-symbols {i j k}
      (- (. m1~i~j m2_j_k) (. m2~i~j m1_j_k)))))

(define $M.inverse
  (lambda [%m]
    (match (tensor-size m) (list integer)
      {[<cons ,2 <cons ,2 <nil>>>
        (T.map (/ $ (M.det m)) (tensor {2 2} {m_2_2 (* -1 m_1_2) (* -1 m_2_1) m_1_1}))]
       [_ undefined]})))

(define $trace (lambda [%t] (with-symbols {i} (contract + t~i_i))))

(define $matrix
  (matcher
    {[<quad-cons $ $ $ $> [math-expr matrix matrix matrix]
      {[$tgt (match (tensor-size tgt) (list integer)
               {[<cons $m <cons $n _>>
                 {[tgt_1_1 tgt_1_[2 n] tgt_[2 m]_1 tgt_[2 m]_[2 n]]}]
                [_ {}]})]}]
     [,$val []
      {[$tgt (if (eq? val tgt) {[]} {})]}]
     [$ [something]
      {[$tgt {tgt}]}]
     }))

;;
;; Determinant
;;

(define $even-and-odd-permutations
  (lambda [$n]
    (match n integer
      {[,2 [{{1 2}} {{2 1}}]]
       [_ (let* {[[$es $os] (even-and-odd-permutations (- n 1))]
                 [$es' (map 1#{@%1 n} es)]
                 [$os' (map 1#{@%1 n} os)]}
            [{@es'
              @(concat (map (lambda [$i] (map (permutate i n $) os')) (between 1 (- n 1))))
              }
             {@os'
              @(concat (map (lambda [$i] (map (permutate i n $) es')) (between 1 (- n 1))))
              }
             ]
            )]})))

(define $permutate
  (lambda [$x $y $xs]
    (match xs (list eq)
      {[<join $hs <cons ,x <join $ms <cons ,y $ts>>>>
        {@hs y @ms x @ts}]
       [<join $hs <cons ,y <join $ms <cons ,x $ts>>>>
        {@hs x @ms y @ts}]})))

(define $M.determinant
  (lambda [%m]
    (match (tensor-size m) (list integer)
      {[<cons $n <cons ,n <nil>>>
        (let {[[$es $os] (even-and-odd-permutations n)]}
          (- (sum (map (lambda [$e]
                         (product (map2 (lambda [$i $j] m_i_j)
                                        (between 1 n)
                                        e)))
                       es))
             (sum (map (lambda [$o]
                         (product (map2 (lambda [$i $j] m_i_j)
                                        (between 1 n)
                                        o)))
                       os))))]
       [_ undefined]})))

(define $M.det M.determinant)

;;;
;;; Eigenvalues and eigenvectors
;;;

(define $M.eigenvalues
  (lambda [%m]
    (match (tensor-size m) (list integer)
      {[<cons ,2 <cons ,2 <nil>>>
        (let {[[$e1 $e2] (q-f (M.det (T.- m (scalar-to-tensor x {2 2}))) x)]}
          {e1 e2})]
       [_ undefined]})))

(define $M.eigenvectors
  (lambda [%m]
    (match (tensor-size m) (list integer)
      {[<cons ,2 <cons ,2 <nil>>>
        (let {[[$e1 $e2] (q-f (M.det (T.- m (scalar-to-tensor x {2 2}))) x)]}
          {[e1 (clear-index (T.- m (scalar-to-tensor e1 {2 2}))_i_1)]
           [e2 (clear-index (T.- m (scalar-to-tensor e2 {2 2}))_i_1)]})
        ]
       [_ undefined]})))

;;
;; LU decomposition
;;

(define $M.LU
  (lambda [%x]
    (match (tensor-size x) (list integer)
      {[<cons ,2 <cons ,2 <nil>>>
        (let* {[$L (generate-tensor 2#(match (compare %1 %2) ordering {[<less> 0] [<equal> 1] [<greater> b_%1_%2]}) {2 2})]
               [$U (generate-tensor 2#(match (compare %1 %2) ordering {[<greater> 0] [_ c_%1_%2]}) {2 2})]
               [$m (M.* L U)]
               [$ret (solve {[m_1_1 x_1_1 c_1_1] [m_1_2 x_1_2 c_1_2]
                             [m_2_1 x_2_1 b_2_1] [m_2_2 x_2_2 c_2_2]})]}
          [(substitute ret L) (substitute ret U)])]
       [<cons ,3 <cons ,3 <nil>>>
        (let* {[$L (generate-tensor 2#(match (compare %1 %2) ordering {[<less> 0] [<equal> 1] [<greater> b_%1_%2]}) {3 3})]
               [$U (generate-tensor 2#(match (compare %1 %2) ordering {[<greater> 0] [_ c_%1_%2]}) {3 3})]
               [$m (M.* L U)]
               [$ret (solve {[m_1_1 x_1_1 c_1_1] [m_1_2 x_1_2 c_1_2] [m_1_3 x_1_3 c_1_3]
                             [m_2_1 x_2_1 b_2_1] [m_2_2 x_2_2 c_2_2] [m_2_3 x_2_3 c_2_3]
                             [m_3_1 x_3_1 b_3_1] [m_3_2 x_3_2 b_3_2] [m_3_3 x_3_3 c_3_3]})]}
          [(substitute ret L) (substitute ret U)])]
       [_ undefined]})))