Library matrix

(c) Copyright Microsoft Corporation and Inria. All rights reserved. 
Basic linear algebra : definition of the Matrix type. Matrix is defined   
as a double indexed list of coefficients. This is done by using the       
finfun structure. The file contains the definitions of:                   
  fun_of_matrix : Coercion from matrix to FunClass. Allow to use matrices 
               as functions.                                              
  matrix_of_fun : the construction operators of a matrix form a given     
                  function. This is the RECOMMENDED interface to build    
                  an element of matrix type.                              
It defines also:                                                          
- row and column operations :  cut, drop and swap                         
- block operation : left cut, right cut and paste of matrices             
- trace : \tr A                                                           
- determinant : \det A. The definition is done with Leibniz formula       
- adjugate matrix : \adj A                                                
- algebraic operation for group, ring, module and unital ring             
- The Canonicals Structures for this algebraic structures are defined     
- LUP matrix decomposition                                                
In addition to the lemmas relevant to these definitions, this file also   
contains proofs :                                                         
- Determinant multilinear property                                        
- Laplace formulas : expand_det_row & expand_det_col                      
- Cramer rule : mulmx_adjr & mulmx_adjl                                   

Import Prenex Implicits.

Import GroupScope.
Import GRing.Theory.

Reserved Notation "''M_' n" (at level 8, n at level 2, format "''M_' n").
Reserved Notation "''M_' ( n )" (at level 8, only parsing).
Reserved Notation "''M_' ( m , n )" (at level 8, format "''M_' ( m , n )").

Reserved Notation "\matrix_ ( i , j ) E"
  (at level 36, E at level 36, i, j at level 50,
   format "\matrix_ ( i , j ) E").
Reserved Notation "\matrix_ ( i < m , j < n ) E"
  (at level 36, E at level 36, i, m, j, n at level 50,
   format "\matrix_ ( i < m , j < n ) E").
Reserved Notation "\matrix_ ( i , j < n ) E"
  (at level 36, E at level 36, i, j, n at level 50,
   format "\matrix_ ( i , j < n ) E").

Reserved Notation "x %:M" (at level 8, format "x %:M").
Reserved Notation "x *m: A" (at level 40, format "x *m: A").
Reserved Notation "A *m B" (at level 40, format "A *m B").
Reserved Notation "A ^T" (at level 8, format "A ^T").
Reserved Notation "\tr A" (at level 10, A at level 8, format "\tr A").
Reserved Notation "\det A" (at level 10, A at level 8, format "\det A").
Reserved Notation "\adj A" (at level 10, A at level 8, format "\adj A").

Delimit Scope matrix_scope with M.

Notation Local simp := (Monoid.Theory.simpm, oppr0).

Open Local Scope ring_scope.
Open Local Scope matrix_scope.


Section MatrixDef.

Variable R : Type.
Variables m n : nat.

Basic linear algebra (matrices).                                       
We use dependent types (ordinals) for the indices so that ranges are   
mostly inferred automatically                                          

CoInductive matrix : predArgType := Matrix of {ffun 'I_m * 'I_n -> R}.

Definition mx_val A := let: Matrix g := A in g.

Lemma mx_valK : cancel mx_val Matrix. Qed.

Definition matrix_of_fun F := locked Matrix [ffun ij => F ij.1 ij.2].

Definition fun_of_matrix A (i : 'I_m) (j : 'I_n) := mx_val A (i, j).

Coercion fun_of_matrix : matrix >-> Funclass.

Lemma mxE : forall F, matrix_of_fun F =2 F.

Lemma matrixP : forall A B : matrix, A =2 B <-> A = B.

End MatrixDef.


Notation "''M_' n" := (matrix _ n n) : type_scope.
Notation "''M_' ( n )" := 'M_n (only parsing) : type_scope.
Notation "''M_' ( m , n )" := (matrix _ m n) : type_scope.

Notation "\matrix_ ( i < m , j < n ) E" :=
  (@matrix_of_fun _ m n (fun i j => E)) (only parsing) : ring_scope.

Notation "\matrix_ ( i , j < n ) E" :=
  (\matrix_(i < n, j < n) E) (only parsing) : ring_scope.

Notation "\matrix_ ( i , j ) E" := (\matrix_(i < _, j < _) E) : ring_scope.

Definition matrix_eqMixin (R : eqType) m n :=
  CanEqMixin (@mx_valK R m n).
Canonical Structure matrix_eqType R m n:=
  Eval hnf in EqType (matrix_eqMixin R m n).
Definition matrix_choiceMixin (R : choiceType) m n :=
  CanChoiceMixin (@mx_valK R m n).
Canonical Structure matrix_choiceType R m n :=
  Eval hnf in ChoiceType (matrix_choiceMixin R m n).


Section Slicing.

Variable R : Type.

Row/Column sub matrices of a matrix 
Definition mx_row m n i0 (A : 'M_(m, n)) :=
  \matrix_(i < 1, j < n) (A i0 j : R).
Definition mx_col m n j0 (A : 'M_(m, n)) :=
  \matrix_(i < m, j < 1) (A i j0 : R).

Removing a row/column from a matrix 
Definition mx_row' m n i0 (A : 'M_(m, n)) :=
  \matrix_(i, j) (A (lift i0 i) j : R).
Definition mx_col' m n j0 (A : 'M_(m, n)) :=
  \matrix_(i, j) (A i (lift j0 j) : R).

Swaping two row/column of a matrix 
Definition rswap m n (A : 'M_(m, n)) i1 i2 :=
  \matrix_(i, j) (A (tperm i1 i2 i) j : R).
Definition cswap m n (A : 'M_(m, n)) i1 i2 :=
  \matrix_(i, j) (A i (tperm i1 i2 j) : R).

Transpose 
Definition trmx m n (A : 'M_(m, n)) := \matrix_(i, j) (A j i : R).

Lemma trmxK : forall m n, cancel (@trmx m n) (@trmx n m).

Lemma trmx_inj : forall m n, injective (@trmx m n).

Notation "A ^T" := (trmx A) : ring_scope.

Lemma trmx_row : forall m n i0 (A : 'M_(m, n)),
  (mx_row i0 A)^T = mx_col i0 A^T.

Lemma trmx_row' : forall m n i0 (A : 'M_(m, n)),
  (mx_row' i0 A)^T = mx_col' i0 A^T.

Lemma trmx_col : forall m n j0 (A : 'M_(m, n)),
  (mx_col j0 A)^T = mx_row j0 A^T.

Lemma trmx_col' : forall m n j0 (A : 'M_(m, n)),
  (mx_col' j0 A)^T = mx_row' j0 A^T.

Lemma trmx_cswap : forall m n (A : 'M_(m, n)) i1 i2,
  (cswap A i1 i2)^T = rswap A^T i1 i2.

Lemma trmx_rswap : forall m n (A : 'M_(m, n)) i1 i2,
  (rswap A i1 i2)^T = cswap A^T i1 i2.

Lemma mx_row_id : forall n (A : 'M_(1, n)), mx_row 0 A = A.

Lemma mx_row_eq : forall m1 m2 n i1 i2 (A1 : 'M_(m1, n)) (A2 : 'M_(m2, n)),
  mx_row i1 A1 = mx_row i2 A2 -> A1 i1 =1 A2 i2.

Lemma mx_row'_eq : forall m n i0 (A B : 'M_(m, n)),
  mx_row' i0 A = mx_row' i0 B -> {in predC1 i0, A =2 B}.

Section CutPaste.

Variables m n1 n2 : nat.

Left/Right cuting of colums of a matrix 
The shape of the (dependent) width parameter of the type of A 
determines where the cut is made 
Definition lcutmx (A : 'M_(m, n1 + n2)):=
  \matrix_(i < m, j < n1) (A i (lshift n2 j) : R).

Definition rcutmx (A : 'M_(m, n1 + n2)) :=
  \matrix_(i < m, j < n2) (A i (rshift n1 j) : R).

Concatenating two matrices 
Definition pastemx (A1 : 'M_(m, n1)) (A2 : 'M_(m, n2)) :=
   \matrix_(i < m, j < n1 + n2)
      (match split j with inl j1 => A1 i j1 | inr j2 => A2 i j2 end : R).

Lemma pastemxEl : forall A1 A2 i j, pastemx A1 A2 i (lshift n2 j) = A1 i j.

Lemma pastemxEr : forall A1 A2 i j, pastemx A1 A2 i (rshift n1 j) = A2 i j.

Lemma pastemxKl : forall A1 A2, lcutmx (pastemx A1 A2) = A1.

Lemma pastemxKr : forall A1 A2, rcutmx (pastemx A1 A2) = A2.

Lemma cutmxK : forall A, pastemx (lcutmx A) (rcutmx A) = A.

End CutPaste.

Lemma mx_row_paste : forall m n1 n2 i0 (A1 : 'M_(m, n1)) (A2 : 'M_(m, n2)),
  mx_row i0 (pastemx A1 A2) = pastemx (mx_row i0 A1) (mx_row i0 A2).

Lemma mx_row'_paste : forall m n1 n2 i0 (A1 : 'M_(m, n1)) (A2 : 'M_(m, n2)),
  mx_row' i0 (pastemx A1 A2) = pastemx (mx_row' i0 A1) (mx_row' i0 A2).

Lemma mx_col_lshift : forall m n1 n2 j1 (A1 : 'M_(m, n1)) A2,
  mx_col (lshift n2 j1) (pastemx A1 A2) = mx_col j1 A1.

Lemma mx_col_rshift : forall m n1 n2 j2 A1 (A2 : 'M_(m, n2)),
  mx_col (rshift n1 j2) (pastemx A1 A2) = mx_col j2 A2.

Lemma mx_col'_lshift : forall m n1 n2 j1 (A1 : 'M_(m, n1.+1)) A2,
  mx_col' (lshift n2 j1) (pastemx A1 A2) = pastemx (mx_col' j1 A1) A2.

Lemma mx_col'_rcast : forall n1 n2, 'I_n2 -> (n1 + n2.-1)%N = (n1 + n2).-1.

Lemma paste_mx_col' : forall m n1 n2 j2 A1 (A2 : 'M_(m, n2)),
  pastemx A1 (mx_col' j2 A2)
    = eq_rect _ (matrix R m) (mx_col' (rshift n1 j2) (pastemx A1 A2))
              _ (esym (mx_col'_rcast n1 j2)).

Lemma mx_col'_rshift : forall m n1 n2 j2 A1 (A2 : 'M_(m, n2)),
  mx_col' (rshift n1 j2) (pastemx A1 A2)
    = eq_rect _ (matrix R m) (pastemx A1 (mx_col' j2 A2))
              _ (mx_col'_rcast n1 j2).

Section Block.

Variables m1 m2 n1 n2 : nat.

Building a block matrix from 4 matrices :            
 up left, up right, low left and low right component 
Definition block_mx Aul Aur All Alr : 'M_(m1 + m2, n1 + n2) :=
  (pastemx (pastemx Aul Aur)^T (pastemx All Alr)^T)^T.

Section CutBlock.

Variable A : matrix R (m1 + m2) (n1 + n2).

Definition ulsubmx := lcutmx (lcutmx A^T)^T.
Definition ursubmx := rcutmx (lcutmx A^T)^T.
Definition llsubmx := lcutmx (rcutmx A^T)^T.
Definition lrsubmx := rcutmx (rcutmx A^T)^T.

Lemma submxK : block_mx ulsubmx ursubmx llsubmx lrsubmx = A.

End CutBlock.

Section PasteBlock.

Variables (Aul : matrix R m1 n1) (Aur : matrix R m1 n2).
Variables (All : matrix R m2 n1) (Alr : matrix R m2 n2).

Let A := block_mx Aul Aur All Alr.

Lemma block_mxEul : forall i j, A (lshift m2 i) (lshift n2 j) = Aul i j.
Lemma block_mxKul : ulsubmx A = Aul.

Lemma block_mxEur : forall i j, A (lshift m2 i) (rshift n1 j) = Aur i j.
Lemma block_mxKur : ursubmx A = Aur.

Lemma block_mxEll : forall i j, A (rshift m1 i) (lshift n2 j) = All i j.
Lemma block_mxKll : llsubmx A = All.

Lemma block_mxElr : forall i j, A (rshift m1 i) (rshift n1 j) = Alr i j.
Lemma block_mxKlr : lrsubmx A = Alr.

End PasteBlock.

End Block.

Section TrBlock.

Variables m1 m2 n1 n2 : nat.

Section TrCut.

Variable A : matrix R (m1 + m2) (n1 + n2).

Lemma trmx_ulsub : (ulsubmx A)^T = ulsubmx A^T.

Lemma trmx_ursub : (ursubmx A)^T = llsubmx A^T.

Lemma trmx_llsub : (llsubmx A)^T = ursubmx A^T.

Lemma trmx_lrsub : (lrsubmx A)^T = lrsubmx A^T.

End TrCut.

Lemma trmx_block : forall Aul Aur All Alr,
 (block_mx Aul Aur All Alr)^T =
    block_mx Aul^T All^T Aur^T Alr^T :> 'M_(n1 + n2, m1 + m2).

End TrBlock.

End Slicing.

Notation "A ^T" := (trmx A) : ring_scope.


Definition of operations for matrices over a ring 
Section MatrixAlgebraOps.

Variable R : ringType.

Section ZmodOps.
The Zmodule structure 

Variables m n : nat.
Implicit Types A B C : matrix R m n.

Definition null_mx := \matrix_(i < m, j < n) (0 : R).
Definition oppmx A := \matrix_(i < m, j < n) (- A i j).
Definition addmx A B := \matrix_(i < m, j < n) (A i j + B i j).
Definition scalemx x A := \matrix_(i < m, j < n) (x * A i j).

Lemma addmxA : associative addmx.

Lemma addmxC : commutative addmx.

Lemma add0mx : left_id null_mx addmx.

Lemma addNmx : left_inverse null_mx oppmx addmx.

Definition matrix_zmodMixin := ZmodMixin addmxA addmxC add0mx addNmx.
Canonical Structure matrix_zmodType := Eval hnf in ZmodType matrix_zmodMixin.

Lemma summxE : forall I r (P : pred I) (E : I -> 'M_(m, n)) i j,
  (\sum_(k <- r | P k) E k) i j = \sum_(k <- r | P k) E k i j.

Vector space structure 

Notation "x *m: A" := (scalemx x A) : ring_scope.

Lemma scale0mx : forall A, 0 *m: A = 0.

Lemma scalemx0 : forall x, x *m: 0 = 0.

Lemma scale1mx : forall A, 1 *m: A = A.

Lemma scaleNmx : forall x A, (- x) *m: A = - (x *m: A).

Lemma scalemxN : forall x A, x *m: (- A) = - (x *m: A).

Lemma scalemx_addl : forall x y A, (x + y) *m: A = x *m: A + y *m: A.

Lemma scalemx_addr : forall x A B, x *m: (A + B) = x *m: A + x *m: B.

Lemma scalemx_subl : forall x y A, (x - y) *m: A = x *m: A - y *m: A.

Lemma scalemx_subr : forall x A B, x *m: (A - B) = x *m: A - x *m: B.

Lemma scalemxA : forall x y A, x *m: (y *m: A) = (x * y) *m: A.

Basis 

Definition delta_mx i0 j0 :=
  \matrix_(i < m, j < n) (((i == i0) && (j == j0))%:R : R).

Lemma matrix_sum_delta : forall A,
  A = \sum_(i < m) \sum_(j < n) A i j *m: delta_mx i j.

End ZmodOps.

Notation "x *m: A" := (scalemx x A) : ring_scope.

Lemma trmx0 : forall m n, 0^T = 0 :> matrix R n m.

Lemma trmx_add : forall m n A B, (A + B)^T = A^T + B^T :> matrix R n m.

Lemma trmx_scale : forall m n a A, (a *m: A)^T = a *m: A^T :> matrix R n m.

Lemma mx_row0 : forall m n i0, @mx_row R m n i0 0 = 0.

Lemma mx_col0 : forall m n j0, @mx_col R m n j0 0 = 0.

Lemma mx_row'0 : forall m n i0, @mx_row' R m n i0 0 = 0.

Lemma mx_col'0 : forall m n j0, @mx_col' R m n j0 0 = 0.

Lemma pastemx0 : forall m n1 n2,
  pastemx 0 0 = 0 :> matrix R m (n1 + n2).

Lemma addmx_paste : forall m n1 n2 A1 A2 B1 B2,
  pastemx A1 A2 + pastemx B1 B2 = pastemx (A1 + B1) (A2 + B2)
                              :> matrix R m (n1 + n2).

Lemma scalemx_paste : forall m n1 n2 a A1 A2,
  a *m: pastemx A1 A2 = pastemx (a *m: A1) (a *m: A2) :> matrix R m (n1 + n2).

Lemma block_mx0 : forall m1 m2 n1 n2,
  block_mx 0 0 0 0 = 0 :> matrix R (m1 + m2) (n1 + n2).

Lemma addmx_block : forall m1 m2 n1 n2 Aul Aur All Alr Bul Bur Bll Blr,
  block_mx Aul Aur All Alr + block_mx Bul Bur Bll Blr
    = block_mx (Aul + Bul) (Aur + Bur) (All + Bll) (Alr + Blr)
    :> matrix R (m1 + m2) (n1 + n2).

Lemma scalemx_block : forall m1 m2 n1 n2 a Aul Aur All Alr,
  a *m: block_mx Aul Aur All Alr
     = block_mx (a *m: Aul) (a *m: Aur) (a *m: All) (a *m: Alr)
    :> matrix R (m1 + m2) (n1 + n2).

Scalar matrix : a diagonal matrix with a constant on the diagonal 
Definition scalar_mx n x := \matrix_(i , j < n) (if i == j then x else 0 : R).

Matrix multiplication with the bigops 
Definition mulmx m n p (A : 'M_(m, n)) (B : 'M_(n, p)) :=
  \matrix_(i < m, k < p) \sum_(j < n) (A i j * B j k : R).

Notation "x %:M" := (scalar_mx _ x) : ring_scope.
Notation "A *m B" := (mulmx A B) : ring_scope.

Lemma scalar_mx0 : forall n, 0%:M = 0 :> 'M_n.

Lemma scalar_mx_opp : forall n a, (- a)%:M = - a%:M :> 'M_n.

Lemma scalar_mx_add : forall n a b, (a + b)%:M = a%:M + b%:M :> 'M_n.

Lemma mulmx_scalar : forall m n a A, a%:M *m A = a *m: A :> 'M_(m, n).

Lemma scalar_mx_mul : forall n a b, (a * b)%:M = a%:M *m b%:M :> 'M_n.

Lemma trmx_scalar : forall n a, (a%:M)^T = a%:M :> 'M_n.

Lemma mul1mx : forall m n A, 1%:M *m A = A :> 'M_(m, n).

Lemma mulmx_addl : forall m n p (A1 A2 : 'M_(m, n)) (B : 'M_(n, p)),
  (A1 + A2) *m B = A1 *m B + A2 *m B.

Lemma scalemx_add : forall n a1 a2, (a1 + a2)%:M = a1%:M + a2%:M :> 'M_n.

Lemma scalemxAl : forall m n p a (A : 'M_(m, n)) (B : 'M_(n, p)),
  a *m: (A *m B) = (a *m: A) *m B.

Lemma mul0mx : forall m n p (A : 'M_(n, p)), 0 *m A = 0 :> 'M_(m, p).

Lemma mulmx0 : forall m n p (A : 'M_(m, n)), A *m 0 = 0 :> 'M_(m, p).

Lemma mulmx1 : forall m n (A : 'M_(m, n)), A *m 1%:M = A.

Lemma mulmx_addr : forall m n p (A : 'M_(m, n)) (B1 B2 : 'M_(n, p)),
  A *m (B1 + B2) = A *m B1 + A *m B2.

Lemma mulmxA : forall m n p q (A : 'M_(m, n)) (B : 'M_(n, p)) (C : 'M_(p, q)),
  A *m (B *m C) = A *m B *m C.

The trace, in 1/4 line. 
Definition mx_trace n (A : 'M_n) := \sum_i A i i : R.
Notation "'\tr' A" := (mx_trace A) : ring_scope.

Lemma mx_trace0 : forall n, \tr (0 : 'M_n) = 0.

Lemma mx_trace_scale : forall n a (A : 'M_n), \tr (a *m: A) = a * \tr A.

Lemma mx_trace_scalar : forall n a, \tr (a%:M : 'M_n) = a *+ n.

Lemma mx_trace_add : forall n A B, \tr (A + B : 'M_n) = \tr A + \tr B.

Lemma mx_trace_tr : forall n (A : 'M_n), \tr A^T = \tr A.

Lemma mx_trace_block : forall n1 n2 Aul Aur All Alr,
  \tr (block_mx Aul Aur All Alr : 'M_(n1 + n2)) = \tr Aul + \tr Alr.

Lemma mulmx_paste : forall m n p1 p2 (A : 'M_(m, n)) B1 B2,
  A *m (pastemx B1 B2) = pastemx (A *m B1) (A *m B2) :> 'M_(m, p1 + p2).

Lemma dotmx_paste : forall m n1 n2 p A1 A2 B1 B2,
  (pastemx A1 A2 : 'M_(m, n1 + n2)) *m (pastemx B1 B2 : 'M_(p, n1 + n2))^T
    = A1 *m B1^T + A2 *m B2^T.

Matrix ring Structure 
Section MatrixRing.
Variable n : pos_nat.
Require that n > 0 to avoid the trivial ring 

Lemma matrix_nonzero1 : 1%:M != 0 :> 'M_n.

Definition matrix_ringMixin :=
  RingMixin (@mulmxA n n n n) (@mul1mx n n) (@mulmx1 n n)
            (@mulmx_addl n n n) (@mulmx_addr n n n) matrix_nonzero1.
Canonical Structure matrix_ringType := Eval hnf in RingType matrix_ringMixin.

Lemma mulmxE : forall A B : 'M_n, A *m B = A * B. Qed.

Lemma idmxE : 1%:M = 1 :> 'M_n. Qed.

End MatrixRing.

End MatrixAlgebraOps.

Notation "a *m: A" := (scalemx a A) : ring_scope.
Notation "a %:M" := (scalar_mx _ a) : ring_scope.
Notation "A *m B" := (mulmx A B) : ring_scope.
Notation "\tr A" := (mx_trace A) : ring_scope.


Permutation matrix 
Section PermMatrix.
Variable R : ringType.

Definition perm_mx n (s : 'S_n) :=
  \matrix_(i, j) (if s i == j then 1 else 0 : R).

Definition tperm_mx n i1 i2 := @perm_mx n (tperm i1 i2).

Lemma trmx_perm : forall n (s : 'S_n), (perm_mx s)^T = perm_mx s^-1.

Lemma trmx_tperm : forall n i1 i2, (@tperm_mx n i1 i2)^T = tperm_mx i1 i2.

Lemma mulmx_perm : forall n (s t : 'S_n),
 perm_mx s *m perm_mx t = perm_mx (s * t).

Lemma mul_tperm_mx : forall m n (A : matrix R m n) i1 i2,
  (tperm_mx i1 i2) *m A = rswap A i1 i2.

Lemma perm_mx1 : forall n, perm_mx 1%g = 1%:M :> 'M_n.

Definition is_perm_mx n (A : 'M_n) := existsb s, A == perm_mx s.

Lemma is_perm_mxP : forall n (A : 'M_n),
  reflect (exists s, A = perm_mx s) (is_perm_mx A).

Lemma perm_mx_is_perm : forall n (s : 'S_n), is_perm_mx (perm_mx s).

Lemma is_perm_mxMl : forall n (A B : 'M_n),
  is_perm_mx A -> is_perm_mx (A *m B) = is_perm_mx B.

Lemma is_perm_mxMr : forall n (A B : 'M_n),
  is_perm_mx B -> is_perm_mx (A *m B) = is_perm_mx A.

Definitions and lemmas on permutations lifting : useful for Cramer proof 
and LUP decomposition 
Definition lift_perm_fun n i j (s : 'S_n) k :=
  if @unlift n.+1 i k is Some k' then @lift n.+1 j (s k') else j.

Lemma lift_permK : forall n i j s,
  cancel (@lift_perm_fun n i j s) (lift_perm_fun j i s^-1%g).

Definition lift_perm n i j s := perm (can_inj (@lift_permK n i j s)).

Lemma lift_perm_id : forall n i j s, lift_perm i j s i = j :> 'I_n.+1.

Lemma lift_perm_lift : forall n i j s k,
  lift_perm i j s (lift i k) = lift j (s k) :> 'I_n.+1.

Lemma lift_permM : forall n i j k s t,
  (@lift_perm n i j s * lift_perm j k t)%g = lift_perm i k (s * t)%g.

Lemma lift_perm1 : forall n i, @lift_perm n i i 1 = 1%g.

Lemma lift_permV : forall n i j s,
  (@lift_perm n i j s)^-1%g = lift_perm j i s^-1.

Lemma odd_lift_perm : forall n i j s,
  @lift_perm n i j s = odd i (+) odd j (+) s :> bool.

Variable n : nat.
Implicit Type s t : 'S_n.
Implicit Type A : matrix R n n.

Definition lift0_perm s := lift_perm 0 0 s.

Lemma lift0_perm0 : forall s, lift0_perm s 0 = 0.

Lemma rshift1 : @rshift 1 n =1 lift (0 : 'I_n.+1).

Lemma split1 : forall i : 'I_n.+1,
  @split 1 n i = oapp (@inr _ _) (inl _ 0) (unlift 0 i).

Lemma lift0_perm_lift : forall s (i : 'I_n),
  lift0_perm s (lift (0 : 'I_n.+1) i) = lift (0 : 'I_n.+1)(s i).

Lemma lift0_permK : forall s, cancel (lift0_perm s) (lift0_perm s^-1).

Lemma lift0_perm_eq0 : forall s i, (lift0_perm s i == 0) = (i == 0).

Definition lift0_mx A := block_mx (1 : 'M_(1, 1)) 0 0 A.

Lemma lift0_mx_perm : forall s,
  lift0_mx (perm_mx s) = perm_mx (lift0_perm s).

Lemma lift0_mx_is_perm : forall s, is_perm_mx (lift0_mx (perm_mx s)).

Lemma is_perm_mx1 : is_perm_mx (1%:M : 'M_n).

End PermMatrix.

Section TrMul.

Variable R : ringType.

Lemma trmx_mul_rev : forall m n p (A : matrix R m n) (B : matrix R n p),
  (A *m B)^T = (B^T : matrix (RevRingType R) p n) *m A^T.

Lemma mulmx_block : forall m1 m2 n1 n2 p1 p2 Aul Aur All Alr Bul Bur Bll Blr,
  (block_mx Aul Aur All Alr : matrix R (m1 + m2) (n1 + n2))
   *m (block_mx Bul Bur Bll Blr : 'M_(n1 + n2, p1 + p2))
    = block_mx (Aul *m Bul + Aur *m Bll) (Aul *m Bur + Aur *m Blr)
               (All *m Bul + Alr *m Bll) (All *m Bur + Alr *m Blr).

Lemma mul_mx_tperm : forall m n (A : matrix R m n) i1 i2,
  A *m (tperm_mx R i1 i2) = cswap A i1 i2.

End TrMul.


Lemmas requiring that the coefficients are in a commutative ring 
Section ComMatrix.

Variable R : comRingType.

Lemma trmx_mul : forall m n p (A : matrix R m n) (B : 'M_(n, p)),
  (A *m B)^T = B^T *m A^T.

Lemma scalemxAr : forall m n p a (A : matrix R m n) (B : 'M_(n, p)),
  a *m: (A *m B) = A *m (a *m: B).

Lemma scalar_mx_comm : forall (n : pos_nat) a (A : matrix R n n),
  GRing.comm A a%:M.

Lemma mx_trace_mulC : forall m n (A : matrix R m n) B,
  \tr (A *m B) = \tr (B *m A).

End ComMatrix.

Section Determinant.

Variable R : comRingType.

The determinant, in one line with the Leibniz Formula 
Definition determinant n (A : 'M_n) :=
  \sum_(s : 'S_n) (-1) ^+ s * \prod_i A i (s i) : R.

Notation "'\det' A" := (determinant A).

The cofactor of a matrix on the indexes i and j 
Definition cofactor n A (i j : 'I_n) : R :=
  (-1) ^+ (i + j) * \det (mx_row' i (mx_col' j A)).

The adjugate matrix : defined as the transpose of the matrix of cofactors 
Definition adjugate n A := \matrix_(i, j < n) (cofactor A j i : R).

Lemma determinant_multilinear : forall n (A B C : 'M_n) i0 b c,
    mx_row i0 A = b *m: mx_row i0 B + c *m: mx_row i0 C ->
    mx_row' i0 B = mx_row' i0 A ->
    mx_row' i0 C = mx_row' i0 A ->
  \det A = b * \det B + c * \det C :> R.

Lemma alternate_determinant : forall n (A : 'M_n) i1 i2,
  i1 != i2 -> A i1 =1 A i2 -> \det A = 0.

Lemma det_trmx : forall n (A : 'M_n), \det A^T = \det A.

Lemma det_perm_mx : forall n (s : 'S_n), \det (perm_mx R s) = (-1) ^+s.

Lemma det1 : forall n, \det (1%:M : matrix R n n) = 1.

Lemma det_scalemx : forall n x (A : 'M_n),
  \det (x *m: A) = x ^+ n * \det A.

Lemma det_mulmx : forall n (A B : 'M_n), \det (A *m B) = \det A * \det B.

Laplace expansion lemma 
Lemma expand_cofactor : forall n (A : 'M_n) i j,
  cofactor A i j =
    \sum_(s : 'S_n | s i == j) (-1) ^+ s * \prod_(k | i != k) A k (s k).

Lemma expand_det_row : forall n (A : 'M_n) i0,
  \det A = \sum_j A i0 j * cofactor A i0 j.

Lemma cofactor_tr : forall n (A : 'M_n) i j,
  cofactor A^T i j = cofactor A j i.

Lemma expand_det_col : forall n (A : 'M_n) j0,
  \det A = \sum_i (A i j0 * cofactor A i j0).

Cramer Rule : adjugate on the left 
Lemma mulmx_adjr : forall n (A : 'M_n), A *m adjugate A = (\det A)%:M.

Lemma trmx_adj : forall n (A : 'M_n), (adjugate A)^T = adjugate A^T.

Cramer rule : adjugate on the right 
Lemma mulmx_adjl : forall n (A : 'M_n), adjugate A *m A = (\det A)%:M.

Lemma detM : forall (n : pos_nat) (A B : 'M_n), \det (A * B) = \det A * \det B.

Lemma mx11_scalar : forall A : matrix R 1 1, A = (A 0 0)%:M.

Lemma det_scalar : forall n a, \det (a%:M : 'M_n) = a ^+ n.

Lemma det_scalar1 : forall a, \det (a%:M : 'M_1) = a.

Lemma det_perm_mx_neq0 : forall n s, \det (@perm_mx R n s) != 0.

Lemma det_ublock : forall n1 n2 Aul Aur Alr,
  \det (block_mx Aul Aur 0 Alr : 'M_(n1 + n2)) = \det Aul * \det Alr.

Lemma det_lblock : forall n1 n2 Aul All Alr,
  \det (block_mx Aul 0 All Alr : 'M_(n1 + n2)) = \det Aul * \det Alr.

End Determinant.

Notation "\det A" := (determinant A) : ring_scope.
Notation "\adj A" := (adjugate A) : ring_scope.


Section MatrixInv.

Variables (R : comUnitRingType) (n : pos_nat).

Definition unit_mx : pred 'M_n := fun A => GRing.unit (\det A : R).
Definition invmx A := if unit_mx A then (\det A)^-1 *m: \adj A else A.

Lemma mulVmx : {in unit_mx, left_inverse 1 invmx *%R}.

Lemma mulmxV : {in unit_mx, right_inverse 1 invmx *%R}.

Lemma intro_unit_mx : forall A B : 'M_n, B * A = 1 /\ A * B = 1 -> unit_mx A.

Lemma invmx_out : {in predC unit_mx, invmx =1 id}.

Definition matrix_unitRingMixin :=
  UnitRingMixin mulVmx mulmxV intro_unit_mx invmx_out.
Canonical Structure matrix_unitRing :=
  Eval hnf in UnitRingType matrix_unitRingMixin.

Lemmas requiring that the coefficients are in a unit ring 

Lemma det_invmx : forall A : 'M_n, \det A^-1 = (\det A)^-1.

Lemma trmx_unit : forall A : 'M_n, GRing.unit A^T = GRing.unit A.

Lemma trmx_inv : forall A : 'M_n, A^-1^T = (A^T)^-1.

Lemma perm_mxV : forall s : 'S_n, perm_mx R s^-1%g = (perm_mx R s)^-1.

Lemma is_perm_mxV : forall A, @is_perm_mx R n A -> is_perm_mx A^-1.

End MatrixInv.


Section CormenLUP.

Variable F : fieldType.

Decomposition of the matrix A to P A = L U with :
- P a permutation matrix
- L a lower triangular matrix
- U an upper triangular matrix 


Fixpoint cormen_lup n : let M := 'M_n.+1 in M -> M * M * M :=
  match n return let M := 'M_(1 + n) in M -> M * M * M with
  | 0 => fun A => (1%:M, 1%:M, A)
  | n'.+1 => fun A =>
    let k := odflt (0 : 'I__) (pick [pred k | A k 0 != 0]) in
    let A1 := rswap A 0 k in
    let P1 := tperm_mx F 0 k in
    let Schur := ((A k 0)^-1 *m: llsubmx A1) *m ursubmx A1 in
    let: (P2, L2, U2) := cormen_lup (lrsubmx A1 - Schur) in
    let P := block_mx 1 0 0 P2 * P1 in
    let L := block_mx 1 0 ((A k 0)^-1 *m: (P2 *m llsubmx A1)) L2 in
    let U := block_mx (ulsubmx A1) (ursubmx A1) 0 U2 in
    (P, L, U)
  end.

End CormenLUP.

Section CormenLUPCorrect.

Variable F : fieldType.

Lemma cormen_lup_perm : forall n A, is_perm_mx (@cormen_lup F n A).1.1.

Lemma cormen_lup_correct : forall n A,
  let: (P, L, U) := @cormen_lup F n A in P * A = L * U.

Lemma cormen_lup_detL : forall n A, \det (@cormen_lup F n A).1.2 = 1.

Lemma cormen_lup_lower : forall n A (i j : 'I_n.+1),
  i <= j -> (cormen_lup A).1.2 i j = (i == j)%:R :> F.

Lemma cormen_lup_upper : forall n A (i j : 'I_n.+1),
  j < i -> (cormen_lup A).2 i j = 0 :> F.

End CormenLUPCorrect.