Library mxpoly

   This file provides basic support for formal computation with matrices,   
 mainly results combining matrices and univariate polynomials, such as the  
 Cayley-Hamilton theorem; it also contains an extension of the first order  
 representation of algebra introduced in ssralg (GRing.term/formula).       
      rVpoly v == the little-endian decoding of the row vector v as a       
                  polynomial p = \sum_i (v 0 i)%:P * 'X^i.                  
     poly_rV p == the partial inverse to rVpoly, for polynomials of degree  
                  less than d to 'rV_d (d is inferred from the context).    
 Sylvester_mx p q == the Sylvester matrix of p and q.                       
 resultant p q == the rusultant of p and q, i.e., \det (Sylvester_mx p q).  
   horner_mx A == the morphism from {poly R} to 'M_n (n of the form n'.+1)  
                  mapping a (scalar) polynomial p to the value of its       
                  scalar matrix interpretation at A (this is an instance of 
                  the generic horner_morph construct defined in poly).      
 powers_mx A d == the d x (n ^ 2) matrix whose rows are the mxvec encodings 
                  of the first d powers of A (n of the form n'.+1). Thus,   
                  vec_mx (v *m powers_mx A d) = horner_mx A (rVpoly v).     
   char_poly A == the characteristic polynomial of A.                       
 char_poly_mx A == a matrix whose detereminant is char_poly A.              
   mxminpoly A == the minimal polynomial of A, i.e., the smallest monic     
                  polynomial that annihilates A (A must be nontrivial).     
 degree_mxminpoly A == the (positive) degree of mxminpoly A.                
 mx_inv_horner A == the inverse of horner_mx A for polynomials of degree    
                  smaller than degree_mxminpoly A.                          
 This toolkit for building formal matrix expressions is packaged in the     
 MatrixFormula submodule, and comprises the following:                      
     eval_mx e == GRing.eval lifted to matrices (:= map_mx (GRing.eval e)). 
     mx_term A == GRing.Const lifted to matrices.                           
 mulmx_term A B == the formal product of two matrices of terms.             
 mxrank_form m A == a GRing.formula asserting that the interpretation of    
                  the term matrix A has rank m.                             
 submx_form A B == a GRing.formula asserting that the row space of the      
                  interpretation of the term matrix A is included in the    
                  row space of the interpretation of B.                     
   seq_of_rV v == the seq corresponding to a row vector.                    
     row_env e == the flattening of a tensored environment e : seq 'rV_d.   
 row_var F d k == the term vector of width d such that for e : seq 'rV[F]_d 
                  we have eval e 'X_k = eval_mx (row_env e) (row_var d k).  


Import GRing.Theory.
Import Monoid.Theory.

Open Local Scope ring_scope.

 Row vector <-> bounded degree polynomial bijection 
Section RowPoly.

Variables (R : ringType) (d : nat).
Implicit Types u v : 'rV[R]_d.
Implicit Types p q : {poly R}.

Definition rVpoly v := \poly_(k < d) (if insub k is Some i then v 0 i else 0).
Definition poly_rV p := \row_(i < d) p`_i.

Lemma coef_rVpoly : forall v k,
  (rVpoly v)`_k = if insub k is Some i then v 0 i else 0.

Lemma coef_rVpoly_ord : forall v (i : 'I_d), (rVpoly v)`_i = v 0 i.

Lemma rVpoly_delta : forall i, rVpoly (delta_mx 0 i) = 'X^i.

Lemma rVpolyK : cancel rVpoly poly_rV.

Lemma poly_rV_K : forall p, size p <= d -> rVpoly (poly_rV p) = p.

Lemma poly_rV_is_linear : linear poly_rV.
Canonical Structure poly_rV_additive := Additive poly_rV_is_linear.
Canonical Structure poly_rV_linear := Linear poly_rV_is_linear.

Lemma rVpoly_is_linear : linear rVpoly.
Canonical Structure rVpoly_additive := Additive rVpoly_is_linear.
Canonical Structure rVpoly_linear := Linear rVpoly_is_linear.

End RowPoly.

Implicit Arguments poly_rV [R d].

Section Resultant.

Variables (F : fieldType) (p q : {poly F}).
Implicit Types m d r : {poly F}.


Local Notation band r := (lin1_mx (poly_rV \o r \o* rVpoly)).

Definition Sylvester_mx : 'M[F]_(dq + dp) := col_mx (band p) (band q).

Lemma Sylvester_mxE : forall i j : 'I_(dq + dp),
    let S_ r k := r`_(j - k) *+ (k <= j) in
  Sylvester_mx i j = match split i with inl k => S_ p k | inr k => S_ q k end.

Definition resultant := \det Sylvester_mx.

Lemma resultant_eq0 : (resultant == 0) = (size (gcdp p q) > 1).

End Resultant.

Section HornerMx.

Variables (R : comRingType) (n' : nat).
Local Notation n := n'.+1.
Variable A : 'M[R]_n.
Implicit Types p q : {poly R}.

Definition horner_mx := horner_morph (fun a => scalar_mx_comm a A).
Canonical Structure horner_mx_additive := [additive of horner_mx].
Canonical Structure horner_mx_rmorphism := [rmorphism of horner_mx].

Lemma horner_mx_C : forall a, horner_mx a%:P = a%:M.

Lemma horner_mx_X : horner_mx 'X = A.

Lemma horner_mxZ : scalable horner_mx.

Canonical Structure horner_mx_linear := AddLinear horner_mxZ.
Canonical Structure horner_mx_lrmorphism := [lrmorphism of horner_mx].

Definition powers_mx d := \matrix_(i < d) mxvec (A ^+ i).

Lemma horner_rVpoly : forall m (u : 'rV_m),
  horner_mx (rVpoly u) = vec_mx (u *m powers_mx m).

End HornerMx.

Section CharPoly.

Variables (R : ringType) (n : nat) (A : 'M[R]_n).
Implicit Types p q : {poly R}.

Definition char_poly_mx := 'X%:M - map_mx (@polyC R) A.
Definition char_poly := \det char_poly_mx.

Let size_diagA : size diagA = n.

Let split_diagA :
  exists2 q, \prod_(x <- diagA) ('X - x%:P) + q = char_poly & size q <= n.-1.

Lemma size_char_poly : size char_poly = n.+1.

Lemma char_poly_monic : monic char_poly.

Lemma char_poly_trace : n > 0 -> char_poly`_n.-1 = - \tr A.

Lemma char_poly_det : char_poly`_0 = (- 1) ^+ n * \det A.

End CharPoly.

Lemma mx_poly_ring_isom : forall (R : ringType) n' (n := n'.+1),
  exists phi : {rmorphism 'M[{poly R}]_n -> {poly 'M[R]_n}},
  [/\ bijective phi,
      forall p, phi p%:M = map_poly scalar_mx p,
      forall A, phi (map_mx polyC A) = A%:P
    & forall A i j k, (phi A)`_k i j = (A i j)`_k].

Theorem Cayley_Hamilton : forall (R : comRingType) n' (A : 'M[R]_n'.+1),
  horner_mx A (char_poly A) = 0.

Lemma eigenvalue_root_char : forall (F : fieldType) n (A : 'M[F]_n) a,
  eigenvalue A a = root (char_poly A) a.

Section MinPoly.

Variables (F : fieldType) (n' : nat).
Local Notation n := n'.+1.
Variable A : 'M[F]_n.
Implicit Types p q : {poly F}.

Fact degree_mxminpoly_proof : exists d, \rank (powers_mx A d.+1) <= d.
Definition degree_mxminpoly := ex_minn degree_mxminpoly_proof.
Local Notation d := degree_mxminpoly.
Local Notation Ad := (powers_mx A d).

Lemma mxminpoly_nonconstant : d > 0.

Lemma minpoly_mx1 : (1%:M \in Ad)%MS.

Lemma minpoly_mx_free : row_free Ad.

Lemma horner_mx_mem : forall p, (horner_mx A p \in Ad)%MS.

Definition mx_inv_horner B := rVpoly (mxvec B *m pinvmx Ad).

Lemma mx_inv_horner0 : mx_inv_horner 0 = 0.

Lemma mx_inv_hornerK : forall B,
  (B \in Ad)%MS -> horner_mx A (mx_inv_horner B) = B.

Lemma minpoly_mxM : forall B C, (B \in Ad -> C \in Ad -> B * C \in Ad)%MS.

Lemma minpoly_mx_ring : mxring Ad.

Definition mxminpoly := 'X^d - mx_inv_horner (A ^+ d).
Local Notation p_A := mxminpoly.

Lemma size_mxminpoly : size p_A = d.+1.

Lemma mxminpoly_monic : monic p_A.

Lemma size_mod_mxminpoly : forall p, size (p %% p_A) <= d.

Lemma mx_root_minpoly : horner_mx A p_A = 0.

Lemma horner_rVpolyK : forall u : 'rV_d,
  mx_inv_horner (horner_mx A (rVpoly u)) = rVpoly u.

Lemma horner_mxK : forall p, mx_inv_horner (horner_mx A p) = p %% p_A.

Lemma mxminpoly_min : forall p, horner_mx A p = 0 -> p_A %| p.

Lemma horner_rVpoly_inj : @injective 'M_n 'rV_d (horner_mx A \o rVpoly).

Lemma mxminpoly_linear_is_scalar : (d <= 1) = is_scalar_mx A.

Lemma mxminpoly_dvd_char : p_A %| char_poly A.

Lemma eigenvalue_root_min : forall a, eigenvalue A a = root p_A a.

End MinPoly.

 Parametricity. 
Section MapRing.

Variables (aR rR : ringType) (f : {rmorphism aR -> rR}).
Local Notation "A ^f" := (map_mx f A) : ring_scope.
Local Notation fp := (map_poly f).
Variables (d n : nat) (A : 'M[aR]_n).

Lemma map_rVpoly : forall u : 'rV_d, fp (rVpoly u) = rVpoly u^f.

Lemma map_poly_rV : forall p, (poly_rV p)^f = poly_rV (fp p) :> 'rV_d.

Lemma map_char_poly_mx : map_mx fp (char_poly_mx A) = char_poly_mx A^f.

Lemma map_char_poly : fp (char_poly A) = char_poly A^f.

End MapRing.

Section MapComRing.

Variables (aR rR : comRingType) (f : {rmorphism aR -> rR}).
Local Notation "A ^f" := (map_mx f A) : ring_scope.
Local Notation fp := (map_poly f).
Variables (n' : nat) (A : 'M[aR]_n'.+1).

Lemma map_powers_mx : forall e, (powers_mx A e)^f = powers_mx A^f e.

Lemma map_horner_mx : forall p, (horner_mx A p)^f = horner_mx A^f (fp p).

End MapComRing.

Section MapField.

Variables (aF rF : fieldType) (f : {rmorphism aF -> rF}).
Local Notation "A ^f" := (map_mx f A) : ring_scope.
Local Notation fp := (map_poly f).
Variables (n' : nat) (A : 'M[aF]_n'.+1).

Lemma degree_mxminpoly_map : degree_mxminpoly A^f = degree_mxminpoly A.

Lemma mxminpoly_map : mxminpoly A^f = fp (mxminpoly A).

Lemma map_mx_inv_horner : forall u,
  fp (mx_inv_horner A u) = mx_inv_horner A^f u^f.

End MapField.

 Lifting term, formula, envs and eval to matrices. Wlog, and for the sake  
 of simplicity, we only lift (tensor) envs to row vectors; we can always   
 use mxvec/vec_mx to store and retrieve matrices.                          
 We don't provide definitions for addition, substraction, scaling, etc,    
 because they have simple matrix expressions.                              
Module MatrixFormula.

Section MatrixFormula.

Variable F : fieldType.

Local Notation False := GRing.False.
Local Notation True := GRing.True.
Local Notation And := GRing.And (only parsing).
Local Notation Add := GRing.Add (only parsing).
Local Notation Bool b := (GRing.Bool b%bool).
Local Notation term := (GRing.term F).
Local Notation form := (GRing.formula F).
Local Notation eval := GRing.eval.
Local Notation holds := GRing.holds.
Local Notation qf_form := GRing.qf_form.
Local Notation qf_eval := GRing.qf_eval.

Definition eval_mx (e : seq F) := map_mx (eval e).

Definition mx_term := map_mx (@GRing.Const F).

Lemma eval_mx_term : forall e m n (A : 'M_(m, n)), eval_mx e (mx_term A) = A.

Definition mulmx_term m n p (A : 'M[term]_(m, n)) (B : 'M_(n, p)) :=
  \matrix_(i, k) (\big[Add/0]_j (A i j * B j k))%T.

Lemma eval_mulmx : forall e m n p (A : 'M[term]_(m, n)) (B : 'M_(n, p)),
  eval_mx e (mulmx_term A B) = eval_mx e A *m eval_mx e B.

Local Notation morphAnd := (@big_morph _ _ _ true _ andb).

Let Schur m n (A : 'M[term]_(1 + m, 1 + n)) (a := A 0 0) :=
  \matrix_(i, j) (drsubmx A i j - a^-1 * dlsubmx A i 0%R * ursubmx A 0%R j)%T.

Fixpoint mxrank_form (r m n : nat) : 'M_(m, n) -> form :=
  match m, n return 'M_(m, n) -> form with
  | m'.+1, n'.+1 => fun A : 'M_(1 + m', 1 + n') =>
    let nzA k := A k.1 k.2 != 0 in
    let xSchur k := Schur (xrow k.1 0%R (xcol k.2 0%R A)) in
    let recf k := Bool (r > 0) /\ mxrank_form r.-1 (xSchur k) in
    GRing.Pick nzA recf (Bool (r == 0%N))
  | _, _ => fun _ => Bool (r == 0%N)
  end%T.

Lemma mxrank_form_qf : forall r m n (A : 'M_(m, n)), qf_form (mxrank_form r A).

Lemma eval_mxrank : forall e r m n (A : 'M_(m, n)),
  qf_eval e (mxrank_form r A) = (\rank (eval_mx e A) == r).

Lemma eval_vec_mx : forall e m n (u : 'rV_(m * n)),
  eval_mx e (vec_mx u) = vec_mx (eval_mx e u).

Lemma eval_mxvec : forall e m n (A : 'M_(m, n)),
  eval_mx e (mxvec A) = mxvec (eval_mx e A).

Section Subsetmx.

Variables (m1 m2 n : nat) (A : 'M[term]_(m1, n)) (B : 'M[term]_(m2, n)).

Definition submx_form :=
  \big[And/True]_(r < n.+1) (mxrank_form r (col_mx A B) ==> mxrank_form r B)%T.

Lemma eval_col_mx : forall e,
  eval_mx e (col_mx A B) = col_mx (eval_mx e A) (eval_mx e B).

Lemma submx_form_qf : qf_form submx_form.

Lemma eval_submx : forall e,
  qf_eval e submx_form = (eval_mx e A <= eval_mx e B)%MS.

End Subsetmx.

Section Env.

Variable d : nat.

Definition seq_of_rV (v : 'rV_d) : seq F := fgraph [ffun i => v 0 i].

Lemma size_seq_of_rV : forall v, size (seq_of_rV v) = d.

Lemma nth_seq_of_rV : forall x0 v (i : 'I_d), nth x0 (seq_of_rV v) i = v 0 i.

Definition row_var k : 'rV[term]_d := \row_i ('X_(k * d + i))%T.

Definition row_env (e : seq 'rV_d) := flatten (map seq_of_rV e).

Lemma nth_row_env : forall e k (i : 'I_d), (row_env e)`_(k * d + i) = e`_k 0 i.

Lemma eval_row_var : forall e k,
  eval_mx (row_env e) (row_var k) = e`_k :> 'rV_d.

Definition Exists_row_form k (f : form) :=
  foldr GRing.Exists f (map (fun i : 'I_d => k * d + i)%N (enum 'I_d)).

Lemma Exists_rowP : forall e k f,
  d > 0 ->
   ((exists v : 'rV[F]_d, holds (row_env (set_nth 0 e k v)) f)
      <-> holds (row_env e) (Exists_row_form k f)).

End Env.

End MatrixFormula.

End MatrixFormula.