Library mathcomp.algebra.mxpoly

(* (c) Copyright 2006-2016 Microsoft Corporation and Inria.                  
 Distributed under the terms of CeCILL-B.                                  *)

From HB Require Import structures.
From mathcomp Require Import ssreflect ssrfun ssrbool eqtype ssrnat choice seq.
From mathcomp Require Import div fintype tuple finfun bigop fingroup perm.
From mathcomp Require Import ssralg zmodp matrix mxalgebra poly polydiv.

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 resultant 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 determinant is char_poly A. companionmx p == a matrix whose char_poly is p 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. kermxpoly g p == the kernel of p(g) geigenspace g a == the generalized eigenspace of g for eigenvalue a := kermxpoly g ('X ^ n - a%:P) where g : 'M_n eigenpoly g p <=> p is an eigen polynomial for g, i.e. kermxpoly g p != 0 integralOver RtoK u <-> u is in the integral closure of the image of R under RtoK : R -> K, i.e. u is a root of the image of a monic polynomial in R. algebraicOver FtoE u <-> u : E is algebraic over E; it is a root of the image of a nonzero polynomial under FtoE; as F must be a fieldType, this is equivalent to integralOver FtoE u. integralRange RtoK <-> the integral closure of the image of R contains all of K (:= forall u, integralOver RtoK u). 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). conjmx V f := V *m f *m pinvmx V == the conjugation of f by V, i.e. "the" matrix of f in the basis of row vectors of V. Although this makes sense only when f stabilizes V, the definition can be stated more generally. restrictmx V := conjmx (row_base V) A ~P {in S'} == where P is a base change matrix, A is a matrix, and S is a boolean predicate representing a set of matrices, this states that conjmx P A is in S, which means A is similar to a matrix in S. From the latter, we derive several related notions: A ~P B := A ~P {in pred1 B} A is similar to B, with base change matrix P A ~_{in S} B := exists P, P \in S /\ A ~P B == A is similar to B, with a base change matrix in S A ~_{in S} {in S'} := exists P, P \in S /\ A ~P {in S'} == A is similar to a matrix in the class S', with a base change matrix in S all_simmx_in S As S' == all the matrices in the sequence As are similar to some matrix in the predicate S', with a base change matrix in S.
We also specialize the class S' to diagonalizability: diagonalizable_for P A := A ~P {in is_diag_mx}. diagonalizable_in S A := A ~_{in S} {in is_diag_mx}. diagonalizable A := diagonalizable_in unitmx A. codiagonalizable_in S As := all_simmx_in S As is_diag_mx. codiagonalizable As := codiagonalizable_in unitmx As.
The main results in diagnonalization theory are:
  • diagonalizablePeigen: a matrix is diagonalizable iff there is a sequence of scalars r, such that the sum of the associated eigenspaces is full.
  • diagonalizableP: a matrix is diagonalizable iff its minimal polynomial divides a split polynomial with simple roots.
  • codiagonalizableP: a sequence of matrices are diagonalizable in the same basis iff they are all diagonalizable and commute pairwize.
Naming conventions:
  • p, q are polynomials
  • A, B, C are matrices
  • f, g are matrices that are viewed as linear maps
  • V, W are matrices that are viewed as subspaces

Set Implicit Arguments.

Import GRing.Theory.
Import Monoid.Theory.

Local Open Scope ring_scope.

Import Pdiv.Idomain.
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 v k : (rVpoly v)`_k = if insub k is Some i then v 0 i else 0.

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

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

Lemma rVpolyK : cancel rVpoly poly_rV.

Lemma poly_rV_K p : size p d rVpoly (poly_rV p) = p.

Lemma poly_rV_is_linear : linear poly_rV.

Lemma rVpoly_is_linear : linear rVpoly.

End RowPoly.

Arguments poly_rV {R d}.
Arguments poly_rV_K {R d} [p] le_p_d.

Section Resultant.

Variables (R : ringType) (p q : {poly R}).

Let dS := ((size q).-1 + (size p).-1)%N.

Definition Sylvester_mx : 'M[R]_dS := col_mx (band p) (band q).

Lemma Sylvester_mxE (i j : 'I_dS) :
  let S_ r k := r`_(j - k) *+ (k j) in
  Sylvester_mx i j = match split i with inl kS_ p k | inr kS_ q k end.

Definition resultant := \det Sylvester_mx.

End Resultant.


Lemma resultant_in_ideal (R : comRingType) (p q : {poly R}) :
    size p > 1 size q > 1
  {uv : {poly R} × {poly R} | size uv.1 < size q size uv.2 < size p
  & (resultant p q)%:P = uv.1 × p + uv.2 × q}.

Lemma resultant_eq0 (R : idomainType) (p q : {poly R}) :
  (resultant p q == 0) = (size (gcdp p q) > 1).

Section HornerMx.

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

Section OneMatrix.

Variable A : 'M[R]_n.

Definition horner_mx := horner_morph (comm_mx_scalar^~ A).

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

Lemma horner_mx_X : horner_mx 'X = A.

Lemma horner_mxZ : scalable horner_mx.


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

Lemma horner_rVpoly m (u : 'rV_m) :
  horner_mx (rVpoly u) = vec_mx (u ×m powers_mx m).

End OneMatrix.

Lemma horner_mx_diag (d : 'rV[R]_n) (p : {poly R}) :
  horner_mx (diag_mx d) p = diag_mx (map_mx (horner p) d).

Lemma comm_mx_horner A B p : comm_mx A B comm_mx A (horner_mx B p).

Lemma comm_horner_mx A B p : comm_mx A B comm_mx (horner_mx A p) B.

Lemma comm_horner_mx2 A p q : GRing.comm (horner_mx A p) (horner_mx A q).

End HornerMx.

Lemma horner_mx_stable (K : fieldType) m n p
    (V : 'M[K]_(n.+1, m.+1)) (f : 'M_m.+1) :
  stablemx V f stablemx V (horner_mx f p).


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 diagA := [seq A i i | i <- index_enum _ & true].
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 : char_poly \is monic.

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

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

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

Lemma char_poly_trig {R : comRingType} n (A : 'M[R]_n) : is_trig_mx A
  char_poly A = \prod_(i < n) ('X - (A i i)%:P).

Definition companionmx {R : ringType} (p : seq R) (d := (size p).-1) :=
  \matrix_(i < d, j < d)
    if (i == d.-1 :> nat) then - p`_j else (i.+1 == j :> nat)%:R.

Lemma companionmxK {R : comRingType} (p : {poly R}) :
   p \is monic char_poly (companionmx p) = p.

Lemma mulmx_delta_companion (R : ringType) (p : seq R)
  (i: 'I_(size p).-1) (i_small : i.+1 < (size p).-1):
  delta_mx 0 i ×m companionmx p = delta_mx 0 (Ordinal i_small) :> 'rV__.

Lemma row'_col'_char_poly_mx {R : ringType} m i (M : 'M[R]_m) :
  row' i (col' i (char_poly_mx M)) = char_poly_mx (row' i (col' i M)).

Lemma char_block_diag_mx {R : ringType} m n (A : 'M[R]_m) (B : 'M[R]_n) :
  char_poly_mx (block_mx A 0 0 B) =
  block_mx (char_poly_mx A) 0 0 (char_poly_mx B).

Section MinPoly.

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

Fact degree_mxminpoly_proof : d, \rank (powers_mx A d.+1) d.
Definition degree_mxminpoly := ex_minn degree_mxminpoly_proof.

Lemma mxminpoly_nonconstant : d > 0.

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

Lemma minpoly_mx_free : row_free Ad.

Lemma horner_mx_mem 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 B : (B \in Ad)%MS horner_mx A (mx_inv_horner B) = B.

Lemma minpoly_mxM 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).

Lemma size_mxminpoly : size p_A = d.+1.

Lemma mxminpoly_monic : p_A \is monic.

Lemma size_mod_mxminpoly p : size (p %% p_A) d.

Lemma mx_root_minpoly : horner_mx A p_A = 0.

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

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

Lemma mxminpoly_min p : horner_mx A p = 0 p_A %| p.

Lemma mxminpoly_minP p : reflect (horner_mx A p = 0) (p_A %| p).

Lemma dvd_mxminpoly p : (p_A %| p) = (horner_mx A p == 0).

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

Lemma mxminpoly_linear_is_scalar : (d 1) = is_scalar_mx A.

Lemma mxminpoly_dvd_char : p_A %| char_poly A.

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

Lemma root_mxminpoly a : root p_A a = root (char_poly A) a.

End MinPoly.

Lemma mxminpoly_diag {F : fieldType} {n} (d : 'rV[F]_n.+1)
    (u := undup [seq d 0 i | i <- enum 'I_n.+1]) :
  mxminpoly (diag_mx d) = \prod_(r <- u) ('X - r%:P).

Arguments mx_inv_hornerK {F n' A} [B] AnB.
Arguments horner_rVpoly_inj {F n' A} [u1 u2] eq_u12A : rename.

Parametricity.
Section MapRingMatrix.

Variables (aR rR : ringType) (f : {rmorphism aR rR}).
Variables (d n : nat) (A : 'M[aR]_n).

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

Lemma map_poly_rV 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 MapRingMatrix.

Section MapResultant.

Lemma map_resultant (aR rR : ringType) (f : {rmorphism {poly aR} rR}) p q :
    f (lead_coef p) != 0 f (lead_coef q) != 0
  f (resultant p q)= resultant (map_poly f p) (map_poly f q).

End MapResultant.

Section MapComRing.

Variables (aR rR : comRingType) (f : {rmorphism aR rR}).
Variables (n' : nat) (A : 'M[aR]_n'.+1).

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

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

End MapComRing.

Section MapField.

Variables (aF rF : fieldType) (f : {rmorphism aF rF}).
Variables (n' : nat) (A : 'M[aF]_n'.+1) (p : {poly aF}).

Lemma map_mx_companion (e := congr1 predn (size_map_poly _ _)) :
  (companionmx p)^f = castmx (e, e) (companionmx (fp p)).

Lemma companion_map_poly (e := esym (congr1 predn (size_map_poly _ _))) :
  companionmx (fp p) = castmx (e, e) (companionmx p)^f.

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 u : fp (mx_inv_horner A u) = mx_inv_horner A^f u^f.

End MapField.

Section KernelLemmas.

Variable K : fieldType.

convertible to kermx (horner_mx g p) when n = n.+1
Definition kermxpoly n (g : 'M_n) (p : {poly K}) : 'M_n :=
  kermx ((if n is n.+1 then horner_mx^~ p : 'M_n.+1 'M_n.+1 else \0) g).

Lemma kermxpolyC n (g : 'M_n) c : c != 0 kermxpoly g c%:P = 0.

Lemma kermxpoly1 n (g : 'M_n) : kermxpoly g 1 = 0.

Lemma kermxpolyX n (g : 'M_n) : kermxpoly g 'X = kermx g.

Lemma kermxpoly_min n (g : 'M_n.+1) p :
  mxminpoly g %| p (kermxpoly g p :=: 1)%MS.

Lemma comm_mx_stable_kermxpoly n (f g : 'M_n) (p : {poly K}) : comm_mx f g
  stablemx (kermxpoly f p) g.

Lemma mxdirect_kermxpoly n (g : 'M_n) (p q : {poly K}) :
  coprimep p q (kermxpoly g p :&: kermxpoly g q = 0)%MS.

Lemma kermxpolyM n (g : 'M_n) (p q : {poly K}) : coprimep p q
  (kermxpoly g (p × q) :=: kermxpoly g p + kermxpoly g q)%MS.

Lemma kermxpoly_prod n (g : 'M_n)
    (I : finType) (P : {pred I}) (p_ : I {poly K}) :
  {in P &, i j, j != i coprimep (p_ i) (p_ j)}
  (kermxpoly g (\prod_(i | P i) p_ i) :=: \sum_(i | P i) kermxpoly g (p_ i))%MS.

Lemma mxdirect_sum_kermx n (g : 'M_n)
    (I : finType) (P : {pred I}) (p_ : I {poly K}) :
  {in P &, i j, j != i coprimep (p_ i) (p_ j)}
  mxdirect (\sum_(i | P i) kermxpoly g (p_ i))%MS.

Lemma eigenspace_poly n a (f : 'M_n) :
  eigenspace f a = kermxpoly f ('X - a%:P).

Definition geigenspace n (g : 'M_n) a := kermxpoly g (('X - a%:P) ^+ n).

Lemma geigenspaceE n' (g : 'M_n'.+1) a :
  geigenspace g a = kermx ((g - a%:M) ^+ n'.+1).

Lemma eigenspace_sub_geigen n (g : 'M_n) a :
  (eigenspace g a geigenspace g a)%MS.

Lemma mxdirect_sum_geigenspace
    (I : finType) (n : nat) (g : 'M_n) (P : {pred I}) (a_ : I K) :
  {in P &, injective a_} mxdirect (\sum_(i | P i) geigenspace g (a_ i)).

Definition eigenpoly n (g : 'M_n) : pred {poly K} :=
  (fun pkermxpoly g p != 0).

Lemma eigenpolyP n (g : 'M_n) (p : {poly K}) :
  reflect (exists2 v : 'rV_n, (v kermxpoly g p)%MS & v != 0) (eigenpoly g p).

Lemma eigenvalue_poly n a (f : 'M_n) : eigenvalue f a = eigenpoly f ('X - a%:P).

Lemma comm_mx_stable_geigenspace n (f g : 'M_n) a : comm_mx f g
  stablemx (geigenspace f a) g.

End KernelLemmas.

Section MapKermxPoly.
Variables (aF rF : fieldType) (f : {rmorphism aF rF}).

Lemma map_kermxpoly (n : nat) (g : 'M_n) (p : {poly aF}) :
  map_mx f (kermxpoly g p) = kermxpoly (map_mx f g) (map_poly f p).

Lemma map_geigenspace (n : nat) (g : 'M_n) (a : aF) :
  map_mx f (geigenspace g a) = geigenspace (map_mx f g) (f a).

Lemma eigenpoly_map n (g : 'M_n) (p : {poly aF}) :
  eigenpoly (map_mx f g) (map_poly f p) = eigenpoly g p.

End MapKermxPoly.

Section IntegralOverRing.

Definition integralOver (R K : ringType) (RtoK : R K) (z : K) :=
  exists2 p, p \is monic & root (map_poly RtoK p) z.

Definition integralRange R K RtoK := z, @integralOver R K RtoK z.

Variables (B R K : ringType) (BtoR : B R) (RtoK : {rmorphism R K}).

Lemma integral_rmorph x :
  integralOver BtoR x integralOver (RtoK \o BtoR) (RtoK x).

Lemma integral_id x : integralOver RtoK (RtoK x).

Lemma integral_nat n : integralOver RtoK n%:R.

Lemma integral0 : integralOver RtoK 0.

Lemma integral1 : integralOver RtoK 1.

Lemma integral_poly (p : {poly K}) :
  ( i, integralOver RtoK p`_i) {in p : seq K, integralRange RtoK}.

End IntegralOverRing.

Section IntegralOverComRing.

Variables (R K : comRingType) (RtoK : {rmorphism R K}).

Lemma integral_horner_root w (p q : {poly K}) :
    p \is monic root p w
    {in p : seq K, integralRange RtoK} {in q : seq K, integralRange RtoK}
  integralOver RtoK q.[w].

Lemma integral_root_monic u p :
    p \is monic root p u {in p : seq K, integralRange RtoK}
  integralOver RtoK u.

Let integral0_RtoK := integral0 RtoK.
Let integral1_RtoK := integral1 RtoK.
Let monicXsubC_K := @monicXsubC K.
Hint Resolve integral0_RtoK integral1_RtoK monicXsubC_K : core.

Let XsubC0 (u : K) : root ('X - u%:P) u.
Let intR_XsubC u :
  integralOver RtoK (- u) {in 'X - u%:P : seq K, integralRange RtoK}.

Lemma integral_opp u : integralOver RtoK u integralOver RtoK (- u).

Lemma integral_horner (p : {poly K}) u :
    {in p : seq K, integralRange RtoK} integralOver RtoK u
  integralOver RtoK p.[u].

Lemma integral_sub u v :
  integralOver RtoK u integralOver RtoK v integralOver RtoK (u - v).

Lemma integral_add u v :
  integralOver RtoK u integralOver RtoK v integralOver RtoK (u + v).

Lemma integral_mul u v :
  integralOver RtoK u integralOver RtoK v integralOver RtoK (u × v).

End IntegralOverComRing.

Section IntegralOverField.

Variables (F E : fieldType) (FtoE : {rmorphism F E}).

Definition algebraicOver (fFtoE : F E) u :=
  exists2 p, p != 0 & root (map_poly fFtoE p) u.

Notation mk_mon p := ((lead_coef p)^-1 *: p).

Lemma integral_algebraic u : algebraicOver FtoE u integralOver FtoE u.

Lemma algebraic_id a : algebraicOver FtoE (FtoE a).

Lemma algebraic0 : algebraicOver FtoE 0.

Lemma algebraic1 : algebraicOver FtoE 1.

Lemma algebraic_opp x : algebraicOver FtoE x algebraicOver FtoE (- x).

Lemma algebraic_add x y :
  algebraicOver FtoE x algebraicOver FtoE y algebraicOver FtoE (x + y).

Lemma algebraic_sub x y :
  algebraicOver FtoE x algebraicOver FtoE y algebraicOver FtoE (x - y).

Lemma algebraic_mul x y :
  algebraicOver FtoE x algebraicOver FtoE y algebraicOver FtoE (x × y).

Lemma algebraic_inv u : algebraicOver FtoE u algebraicOver FtoE u^-1.

Lemma algebraic_div x y :
  algebraicOver FtoE x algebraicOver FtoE y algebraicOver FtoE (x / y).

Lemma integral_inv x : integralOver FtoE x integralOver FtoE x^-1.

Lemma integral_div x y :
  integralOver FtoE x integralOver FtoE y integralOver FtoE (x / y).

Lemma integral_root p u :
    p != 0 root p u {in p : seq E, integralRange FtoE}
  integralOver FtoE u.

End IntegralOverField.

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, subtraction, scaling, etc, because they have simple matrix expressions.
Module MatrixFormula.

Section MatrixFormula.

Variable F : fieldType.


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

Definition mx_term := @map_mx F term GRing.Const.

Lemma eval_mx_term 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 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.


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'.+1fun 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 r m n (A : 'M_(m, n)) : qf_form (mxrank_form r A).

Lemma eval_mxrank 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 e m n (u : 'rV_(m × n)) :
  eval_mx e (vec_mx u) = vec_mx (eval_mx e u).

Lemma eval_mxvec 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 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 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 v : size (seq_of_rV v) = d.

Lemma nth_seq_of_rV 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 e k (i : 'I_d) : (row_env e)`_(k × d + i) = e`_k 0 i.

Lemma eval_row_var 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 (codom (fun i : 'I_dk × d + i)%N).

Lemma Exists_rowP e k f :
  d > 0
   (( 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.

Section ConjMx.
Context {F : fieldType}.

Definition conjmx (m n : nat)
 (V : 'M_(m, n)) (f : 'M[F]_n) : 'M_m := V ×m f ×m pinvmx V.
Notation restrictmx V := (conjmx (row_base V)).

Lemma stablemx_comp (m n p : nat) (V : 'M[F]_(m, n)) (W : 'M_(n, p)) (f : 'M_p) :
  stablemx W f stablemx V (conjmx W f) stablemx (V ×m W) f.

Lemma stablemx_restrict m n (A : 'M[F]_n) (V : 'M_n) (W : 'M_(m, \rank V)):
  stablemx V A stablemx W (restrictmx V A) = stablemx (W ×m row_base V) A.

Lemma conjmxM (m n : nat) (V : 'M[F]_(m, n)) :
   {in [pred f | stablemx V f] &, {morph conjmx V : f g / f ×m g}}.

Lemma conjMmx (m n p : nat) (V : 'M[F]_(m, n)) (W : 'M_(n, p)) (f : 'M_p) :
  row_free (V ×m W) stablemx W f stablemx V (conjmx W f)
  conjmx (V ×m W) f = conjmx V (conjmx W f).

Lemma conjuMmx (m n : nat) (V : 'M[F]_m) (W : 'M_(m, n)) (f : 'M_n) :
  V \in unitmx row_free W stablemx W f
  conjmx (V ×m W) f = conjmx V (conjmx W f).

Lemma conjMumx (m n : nat) (V : 'M[F]_(m, n)) (W f : 'M_n) :
  W \in unitmx row_free V stablemx V (conjmx W f)
  conjmx (V ×m W) f = conjmx V (conjmx W f).

Lemma conjuMumx (n : nat) (V W f : 'M[F]_n) :
  V \in unitmx W \in unitmx
  conjmx (V ×m W) f = conjmx V (conjmx W f).

Lemma conjmx_scalar (m n : nat) (V : 'M[F]_(m, n)) (a : F) :
  row_free V conjmx V a%:M = a%:M.

Lemma conj0mx (m n : nat) f : conjmx (0 : 'M[F]_(m, n)) f = 0.

Lemma conjmx0 (m n : nat) (V : 'M[F]_(m, n)) : conjmx V 0 = 0.

Lemma conjumx (n : nat) (V : 'M_n) (f : 'M[F]_n) : V \in unitmx
  conjmx V f = V ×m f ×m invmx V.

Lemma conj1mx (n : nat) (f : 'M[F]_n) : conjmx 1%:M f = f.

Lemma conjVmx (n : nat) (V : 'M_n) (f : 'M[F]_n) : V \in unitmx
  conjmx (invmx V) f = invmx V ×m f ×m V.

Lemma conjmxK (n : nat) (V f : 'M[F]_n) :
  V \in unitmx conjmx (invmx V) (conjmx V f) = f.

Lemma conjmxVK (n : nat) (V f : 'M[F]_n) :
  V \in unitmx conjmx V (conjmx (invmx V) f) = f.

Lemma horner_mx_conj m n p (V : 'M[F]_(n.+1, m.+1)) (f : 'M_m.+1) :
   row_free V stablemx V f
   horner_mx (conjmx V f) p = conjmx V (horner_mx f p).

Lemma horner_mx_uconj n p (V : 'M[F]_(n.+1)) (f : 'M_n.+1) :
   V \is a GRing.unit
   horner_mx (V ×m f ×m invmx V) p = V ×m horner_mx f p ×m invmx V.

Lemma horner_mx_uconjC n p (V : 'M[F]_(n.+1)) (f : 'M_n.+1) :
   V \is a GRing.unit
   horner_mx (invmx V ×m f ×m V) p = invmx V ×m horner_mx f p ×m V.

Lemma mxminpoly_conj m n (V : 'M[F]_(m.+1, n.+1)) (f : 'M_n.+1) :
  row_free V stablemx V f mxminpoly (conjmx V f) %| mxminpoly f.

Lemma mxminpoly_uconj n (V : 'M[F]_(n.+1)) (f : 'M_n.+1) :
  V \in unitmx mxminpoly (conjmx V f) = mxminpoly f.

Section fixed_stablemx_space.

Variables (m n : nat).

Implicit Types (V : 'M[F]_(m, n)) (A : 'M[F]_n).
Implicit Types (a : F) (p : {poly F}).

Section Sub.
Variable (k : nat).
Implicit Types (W : 'M[F]_(k, m)).

Lemma sub_kermxpoly_conjmx V f p W : stablemx V f row_free V
  (W kermxpoly (conjmx V f) p)%MS = (W ×m V kermxpoly f p)%MS.

Lemma sub_eigenspace_conjmx V f a W : stablemx V f row_free V
  (W eigenspace (conjmx V f) a)%MS = (W ×m V eigenspace f a)%MS.
End Sub.

Lemma eigenpoly_conjmx V f : stablemx V f row_free V
  {subset eigenpoly (conjmx V f) eigenpoly f}.

Lemma eigenvalue_conjmx V f : stablemx V f row_free V
  {subset eigenvalue (conjmx V f) eigenvalue f}.

Lemma conjmx_eigenvalue a V f : (V eigenspace f a)%MS row_free V
  conjmx V f = a%:M.

End fixed_stablemx_space.
End ConjMx.
Notation restrictmx V := (conjmx (row_base V)).

Definition simmx_to_for {F : fieldType} {m n}
  (P : 'M_(m, n)) A (S : {pred 'M[F]_m}) := S (conjmx P A).
Notation "A ~_ P '{' 'in' S '}'" := (simmx_to_for P A S)
  (at level 70, P at level 0, format "A ~_ P '{' 'in' S '}'") :
  ring_scope.

Notation simmx_for P A B := (A ¬_P {in PredOfSimpl.coerce (pred1 B)}).
Notation "A ~_ P B" := (simmx_for P A B)
  (at level 70, P at level 0, format "A ~_ P B").

Notation simmx_in S A B := (exists2 P, P \in S & A ¬_P B).
Notation "A '~_{in' S '}' B" := (simmx_in S A B)
  (at level 70, format "A '~_{in' S '}' B").

Notation simmx_in_to S A S' := (exists2 P, P \in S & A ¬_P {in S'}).
Notation "A '~_{in' S '}' '{' 'in' S' '}'" := (simmx_in_to S A S')
  (at level 70, format "A '~_{in' S '}' '{' 'in' S' '}'").

Notation all_simmx_in S As S' :=
  (exists2 P, P \in S & all [pred A | A ¬_P {in S'}] As).

Notation diagonalizable_for P A := (A ¬_P {in is_diag_mx}).
Notation diagonalizable_in S A := (A ¬_{in S} {in is_diag_mx}).
Notation diagonalizable A := (diagonalizable_in unitmx A).
Notation codiagonalizable_in S As := (all_simmx_in S As is_diag_mx).
Notation codiagonalizable As := (codiagonalizable_in unitmx As).

Section Simmxity.
Context {F : fieldType}.

Lemma simmxPp m n {P : 'M[F]_(m, n)} {A B} :
  stablemx P A A ¬_P B P ×m A = B ×m P.

Lemma simmxW m n {P : 'M[F]_(m, n)} {A B} : row_free P
  P ×m A = B ×m P A ¬_P B.

Section Simmx.

Context {n : nat}.
Implicit Types (A B P : 'M[F]_n) (As : seq 'M[F]_n) (d : 'rV[F]_n).

Lemma simmxP {P A B} : P \in unitmx
  reflect (P ×m A = B ×m P) (A ¬_P B).

Lemma simmxRL {P A B} : P \in unitmx
  reflect (B = P ×m A ×m invmx P) (A ¬_P B).

Lemma simmxLR {P A B} : P \in unitmx
  reflect (A = conjmx (invmx P) B) (A ¬_P B).

End Simmx.

Lemma simmx_minpoly {n} {P A B : 'M[F]_n.+1} : P \in unitmx
  A ¬_P B mxminpoly A = mxminpoly B.

Lemma diagonalizable_for_row_base m n (P : 'M[F]_(m, n)) (A : 'M_n) :
  diagonalizable_for (row_base P) A = is_diag_mx (restrictmx P A).

Lemma diagonalizable_forPp m n (P : 'M[F]_(m, n)) A :
  reflect ( i j : 'I__, i != j :> nat conjmx P A i j = 0)
          (diagonalizable_for P A).

Lemma diagonalizable_forP n (P : 'M[F]_n) A : P \in unitmx
  reflect ( i j : 'I__, i != j :> nat (P ×m A ×m invmx P) i j = 0)
          (diagonalizable_for P A).

Lemma diagonalizable_forPex {m} {n} {P : 'M[F]_(m, n)} {A} :
  reflect ( D, A ¬_P (diag_mx D)) (diagonalizable_for P A).

Lemma diagonalizable_forLR n {P : 'M[F]_n} {A} : P \in unitmx
  reflect ( D, A = conjmx (invmx P) (diag_mx D)) (diagonalizable_for P A).

Lemma diagonalizable_for_mxminpoly {n} {P A : 'M[F]_n.+1}
  (rs := undup [seq conjmx P A i i | i <- enum 'I_n.+1]) :
  P \in unitmx diagonalizable_for P A
  mxminpoly A = \prod_(r <- rs) ('X - r%:P).

End Simmxity.

Lemma diagonalizable_for_sum (F : fieldType) (m n : nat) (p_ : 'I_n nat)
      (V_ : i, 'M[F]_(p_ i, m)) (A : 'M[F]_m) :
    mxdirect (\sum_i <<V_ i>>)
    ( i, stablemx (V_ i) A)
    ( i, row_free (V_ i))
  diagonalizable_for (\mxcol_i V_ i) A = [ i, diagonalizable_for (V_ i) A].

Section Diag.
Variable (F : fieldType).

Lemma codiagonalizable1 n (A : 'M[F]_n) :
  codiagonalizable [:: A] diagonalizable A.

Definition codiagonalizablePfull n (As : seq 'M[F]_n) :
  codiagonalizable As
   m, exists2 P : 'M_(m, n), row_full P &
        all [pred A | diagonalizable_for P A] As.

Lemma codiagonalizable_on m n (V_ : 'I_n 'M[F]_m) (As : seq 'M[F]_m) :
    (\sum_i V_ i :=: 1%:M)%MS mxdirect (\sum_i V_ i)
    ( i, all (fun Astablemx (V_ i) A) As)
    ( i, codiagonalizable (map (restrictmx (V_ i)) As))
  codiagonalizable As.

Lemma diagonalizable_diag {n} (d : 'rV[F]_n) : diagonalizable (diag_mx d).
Hint Resolve diagonalizable_diag : core.

Lemma diagonalizable_scalar {n} (a : F) : diagonalizable (a%:M : 'M_n).
Hint Resolve diagonalizable_scalar : core.

Lemma diagonalizable0 {n} : diagonalizable (0 : 'M[F]_n).
Hint Resolve diagonalizable0 : core.

Lemma diagonalizablePeigen {n} {A : 'M[F]_n} :
  diagonalizable A
  exists2 rs, uniq rs & (\sum_(r <- rs) eigenspace A r :=: 1%:M)%MS.

Lemma diagonalizableP n' (n := n'.+1) (A : 'M[F]_n) :
  diagonalizable A
  exists2 rs, uniq rs & mxminpoly A %| \prod_(x <- rs) ('X - x%:P).

Lemma diagonalizable_conj_diag m n (V : 'M[F]_(m, n)) (d : 'rV[F]_n) :
  stablemx V (diag_mx d) row_free V diagonalizable (conjmx V (diag_mx d)).

Lemma codiagonalizableP n (As : seq 'M[F]_n) :
  {in As &, A B, comm_mx A B} {in As, A, diagonalizable A}
   codiagonalizable As.

End Diag.