;; ;; 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] (let {[$d (M.det m)]} (generate-tensor 2#(match m matrix {[ (if (even? (+ %1 %2)) (/ (M.det (M.join A B C D)) d) (* -1 (/ (M.det (M.join A B C D)) d)))]}) (tensor-size m))))) (define $trace (lambda [%t] (with-symbols {i} (contract + t~i_i)))) (define $matrix (matcher {[ [math-expr matrix matrix matrix] {[$tgt (match (tensor-size tgt) (list integer) {[> {[tgt_1_1 tgt_1_[2 n] tgt_[2 m]_1 tgt_[2 m]_[2 n]]}] [_ {}]})]}] [ [math-expr matrix matrix matrix matrix] {[$tgt (let* {[$ns (tensor-size tgt)] [$m (nth 1 ns)] [$n (nth 2 ns)]} {[tgt_i_j tgt_[1 (- i 1)]_[1 (- j 1)] tgt_[1 (- i 1)]_[(+ j 1) n] tgt_[(+ i 1) m]_[1 (- j 1)] tgt_[(+ i 1) m]_[(+ j 1) n] ]})]}] [,$val [] {[$tgt (if (eq? val tgt) {[]} {})]}] [$ [something] {[$tgt {tgt}]}] })) (define $M.join (lambda [%A %B %C %D] (let* {[$as (tensor-size A)] [$a1 (nth 1 as)] [$a2 (nth 2 as)] [$bs (tensor-size B)] [$b1 (nth 1 bs)] [$b2 (nth 2 bs)] [$cs (tensor-size C)] [$c1 (nth 1 cs)] [$c2 (nth 2 cs)] [$ds (tensor-size D)] [$d1 (nth 1 ds)] [$d2 (nth 2 ds)] [$m1 (max a1 b1)] [$m2 (max a2 c2)] [$n1 (max c1 d1)] [$n2 (max b2 d2)] } (generate-tensor 2#(match [%1 %2] [integer integer] {[[?(lte? $ a1) ?(lte? $ a2)] A_%1_%2] [[?(lte? $ m1) _] B_%1_(- %2 a2)] [[_ ?(lte? $ m2)] C_(- %1 a1)_%2] [[_ _] D_(- %1 m1)_(- %2 m2)]}) {(+ m1 n1) (+ m2 n2)})))) ;; ;; Determinant ;; (define $even-and-odd-permutations (lambda [$n] (let {[[$es $os] (even-and-odd-permutations' n)]} [(map 1#(lambda [$i] (nth i %1)) es) (map 1#(lambda [$i] (nth i %1)) os)]))) (define $even-and-odd-permutations' (lambda [$n] (match n integer {[,1 [{{1}} {}]] [,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) {[>>> {@hs y @ms x @ts}] [>>> {@hs x @ms y @ts}]}))) (define $M.determinant (lambda [%m] (match (tensor-size m) (list integer) {[>> 1] [>> (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) {[>> (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) {[>> (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) {[>> (let* {[$L (generate-tensor 2#(match (compare %1 %2) ordering {[ 0] [ 1] [ b_%1_%2]}) {2 2})] [$U (generate-tensor 2#(match (compare %1 %2) ordering {[ 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)])] [>> (let* {[$L (generate-tensor 2#(match (compare %1 %2) ordering {[ 0] [ 1] [ b_%1_%2]}) {3 3})] [$U (generate-tensor 2#(match (compare %1 %2) ordering {[ 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]}))) ;; ;; Utility ;; (define $generate-matrix-from-quadratic-expr (lambda [$f $xs] (generate-tensor 2#(coefficient2 f (nth %1 xs) (nth %2 xs)) {(length xs) (length xs)})))