Initial commit

1 files changed, 1109 insertions, 0 deletions

diff --git a/mathcomp/algebra/mxpoly.v b/mathcomp/algebra/mxpoly.v
new file mode 100644
index 0000000..7e889a0
--- /dev/null
+++ b/mathcomp/algebra/mxpoly.v
@@ -0,0 +1,1109 @@
+(* (c) Copyright Microsoft Corporation and Inria. All rights reserved. *)
+Require Import ssreflect ssrfun ssrbool eqtype ssrnat seq div fintype tuple.
+Require Import finfun bigop fingroup perm ssralg zmodp matrix mxalgebra.
+Require Import 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 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.                          *)
+(*  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).  *)
+(******************************************************************************)
+
+Set Implicit Arguments.
+Unset Strict Implicit.
+Unset Printing Implicit Defensive.
+
+Import GRing.Theory.
+Import Monoid.Theory.
+
+Open Local 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.
+Proof. by rewrite coef_poly; case: insubP => [i ->|]; rewrite ?if_same. Qed.
+
+Lemma coef_rVpoly_ord v (i : 'I_d) : (rVpoly v)`_i = v 0 i.
+Proof. by rewrite coef_rVpoly valK. Qed.
+
+Lemma rVpoly_delta i : rVpoly (delta_mx 0 i) = 'X^i.
+Proof.
+apply/polyP=> j; rewrite coef_rVpoly coefXn.
+case: insubP => [k _ <- | j_ge_d]; first by rewrite mxE.
+by case: eqP j_ge_d => // ->; rewrite ltn_ord.
+Qed.
+
+Lemma rVpolyK : cancel rVpoly poly_rV.
+Proof. by move=> u; apply/rowP=> i; rewrite mxE coef_rVpoly_ord. Qed.
+
+Lemma poly_rV_K p : size p <= d -> rVpoly (poly_rV p) = p.
+Proof.
+move=> le_p_d; apply/polyP=> k; rewrite coef_rVpoly.
+case: insubP => [i _ <- | ]; first by rewrite mxE.
+by rewrite -ltnNge => le_d_l; rewrite nth_default ?(leq_trans le_p_d).
+Qed.
+
+Lemma poly_rV_is_linear : linear poly_rV.
+Proof. by move=> a p q; apply/rowP=> i; rewrite !mxE coefD coefZ. Qed.
+Canonical poly_rV_additive := Additive poly_rV_is_linear.
+Canonical poly_rV_linear := Linear poly_rV_is_linear.
+
+Lemma rVpoly_is_linear : linear rVpoly.
+Proof.
+move=> a u v; apply/polyP=> k; rewrite coefD coefZ !coef_rVpoly.
+by case: insubP => [i _ _ | _]; rewrite ?mxE // mulr0 addr0.
+Qed.
+Canonical rVpoly_additive := Additive rVpoly_is_linear.
+Canonical rVpoly_linear := Linear rVpoly_is_linear.
+
+End RowPoly.
+
+Implicit Arguments poly_rV [R d].
+Prenex Implicits rVpoly poly_rV.
+
+Section Resultant.
+
+Variables (R : ringType) (p q : {poly R}).
+
+Let dS := ((size q).-1 + (size p).-1)%N.
+Local Notation band r := (lin1_mx (poly_rV \o r \o* rVpoly)).
+
+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 k => S_ p k | inr k => S_ q k end.
+Proof.
+move=> S_; rewrite mxE; case: {i}(split i) => i; rewrite !mxE /=;
+  by rewrite rVpoly_delta coefXnM ltnNge if_neg -mulrb.
+Qed.
+
+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}.
+Proof.
+move=> p_nc q_nc; pose dp := (size p).-1; pose dq := (size q).-1.
+pose S := Sylvester_mx p q; pose dS := (dq + dp)%N.
+have dS_gt0: dS > 0 by rewrite /dS /dq -(subnKC q_nc).
+pose j0 := Ordinal dS_gt0. 
+pose Ss0 := col_mx (p *: \col_(i < dq) 'X^i) (q *: \col_(i < dp) 'X^i).
+pose Ss := \matrix_(i, j) (if j == j0 then Ss0 i 0 else (S i j)%:P).
+pose u ds s := \sum_(i < ds) cofactor Ss (s i) j0 * 'X^i.
+exists (u _ (lshift dp), u _ ((rshift dq) _)).
+  suffices sz_u ds s: ds > 1 -> size (u ds.-1 s) < ds by rewrite !sz_u.
+  move/ltn_predK=> {2}<-; apply: leq_trans (size_sum _ _ _) _.
+  apply/bigmax_leqP=> i _.
+  have ->: cofactor Ss (s i) j0 = (cofactor S (s i) j0)%:P.
+    rewrite rmorphM rmorph_sign -det_map_mx; congr (_ * \det _).
+    by apply/matrixP=> i' j'; rewrite !mxE.
+  apply: leq_trans (size_mul_leq _ _) (leq_trans _ (valP i)).
+  by rewrite size_polyC size_polyXn addnS /= -add1n leq_add2r leq_b1.
+transitivity (\det Ss); last first.
+  rewrite (expand_det_col Ss j0) big_split_ord !big_distrl /=.
+  by congr (_ + _); apply: eq_bigr => i _;
+    rewrite mxE eqxx (col_mxEu, col_mxEd) !mxE mulrC mulrA mulrAC.
+pose S_ j1 := map_mx polyC (\matrix_(i, j) S i (if j == j0 then j1 else j)).
+pose Ss0_ i dj := \poly_(j < dj) S i (insubd j0 j).
+pose Ss_ dj := \matrix_(i, j) (if j == j0 then Ss0_ i dj else (S i j)%:P).
+have{Ss u} ->: Ss = Ss_ dS.
+  apply/matrixP=> i j; rewrite mxE [in X in _ = X]mxE; case: (j == j0) => {j}//.
+  apply/polyP=> k; rewrite coef_poly Sylvester_mxE mxE.
+  have [k_ge_dS | k_lt_dS] := leqP dS k.
+    case: (split i) => {i}i; rewrite !mxE coefMXn;
+    case: ifP => // /negbT; rewrite -ltnNge ltnS => hi.
+      apply: (leq_sizeP _ _ (leqnn (size p))); rewrite -(ltn_predK p_nc).
+      by rewrite ltn_subRL (leq_trans _ k_ge_dS) // ltn_add2r.
+    - apply: (leq_sizeP _ _ (leqnn (size q))); rewrite -(ltn_predK q_nc).
+      by rewrite ltn_subRL (leq_trans _ k_ge_dS) // addnC ltn_add2l.
+  by rewrite insubdK //; case: (split i) => {i}i;
+     rewrite !mxE coefMXn; case: leqP.
+elim: {-2}dS (leqnn dS) (dS_gt0) => // dj IHj dj_lt_dS _.
+pose j1 := Ordinal dj_lt_dS; pose rj0T (A : 'M[{poly R}]_dS) := row j0 A^T.
+have: rj0T (Ss_ dj.+1) = 'X^dj *: rj0T (S_ j1) + 1 *: rj0T (Ss_ dj).
+  apply/rowP=> i; apply/polyP=> k; rewrite scale1r !(Sylvester_mxE, mxE) eqxx.
+  rewrite coefD coefXnM coefC !coef_poly ltnS subn_eq0 ltn_neqAle andbC.
+  case: (leqP k dj) => [k_le_dj | k_gt_dj] /=; last by rewrite addr0.
+  rewrite Sylvester_mxE insubdK; last exact: leq_ltn_trans (dj_lt_dS).
+  by case: eqP => [-> | _]; rewrite (addr0, add0r).
+rewrite -det_tr => /determinant_multilinear->;
+  try by apply/matrixP=> i j; rewrite !mxE eq_sym (negPf (neq_lift _ _)).
+have [dj0 | dj_gt0] := posnP dj; rewrite ?dj0 !mul1r.
+  rewrite !det_tr det_map_mx addrC (expand_det_col _ j0) big1 => [|i _].
+    rewrite add0r; congr (\det _)%:P.
+    apply/matrixP=> i j; rewrite [in X in _ = X]mxE; case: eqP => // ->.
+    by congr (S i _); apply: val_inj.
+  by rewrite mxE /= [Ss0_ _ _]poly_def big_ord0 mul0r.
+have /determinant_alternate->: j1 != j0 by rewrite -val_eqE -lt0n.
+  by rewrite mulr0 add0r det_tr IHj // ltnW.
+by move=> i; rewrite !mxE if_same.
+Qed.
+
+Lemma resultant_eq0 (R : idomainType) (p q : {poly R}) :
+  (resultant p q == 0) = (size (gcdp p q) > 1).
+Proof.
+have dvdpp := dvdpp; set r := gcdp p q.
+pose dp := (size p).-1; pose dq := (size q).-1.
+have /andP[r_p r_q]: (r %| p) && (r %| q) by rewrite -dvdp_gcd.
+apply/det0P/idP=> [[uv nz_uv] | r_nonC].
+  have [p0 _ | p_nz] := eqVneq p 0.
+    have: dq + dp > 0.
+      rewrite lt0n; apply: contraNneq nz_uv => dqp0.
+      by rewrite dqp0 in uv *; rewrite [uv]thinmx0.
+    by rewrite /dp /dq /r p0 size_poly0 addn0 gcd0p -subn1 subn_gt0.
+  do [rewrite -[uv]hsubmxK -{1}row_mx0 mul_row_col !mul_rV_lin1 /=] in nz_uv *.
+  set u := rVpoly _; set v := rVpoly _; pose m := gcdp (v * p) (v * q).
+  have lt_vp: size v < size p by rewrite (polySpred p_nz) ltnS size_poly.
+  move/(congr1 rVpoly)/eqP; rewrite -linearD linear0 poly_rV_K; last first.
+    rewrite (leq_trans (size_add _ _)) // geq_max.
+    rewrite !(leq_trans (size_mul_leq _ _)) // -subn1 leq_subLR.
+      by rewrite addnC addnA leq_add ?leqSpred ?size_poly.
+    by rewrite addnCA leq_add ?leqSpred ?size_poly.
+  rewrite addrC addr_eq0 => /eqP vq_up.
+  have nz_v: v != 0.
+    apply: contraNneq nz_uv => v0; apply/eqP.
+    congr row_mx; apply: (can_inj (@rVpolyK _ _)); rewrite linear0 // -/u.
+    by apply: contra_eq vq_up; rewrite v0 mul0r -addr_eq0 add0r => /mulf_neq0->.
+  have r_nz: r != 0 := dvdpN0 r_p p_nz.
+  have /dvdpP [[c w] /= nz_c wv]: v %| m by rewrite dvdp_gcd !dvdp_mulr.
+  have m_wd d: m %| v * d -> w %| d.
+    case/dvdpP=> [[k f]] /= nz_k /(congr1 ( *:%R c)).
+    rewrite mulrC scalerA scalerAl scalerAr wv mulrA => /(mulIf nz_v)def_fw.
+    by apply/dvdpP; exists (c * k, f); rewrite //= mulf_neq0.
+  have w_r: w %| r by rewrite dvdp_gcd !m_wd ?dvdp_gcdl ?dvdp_gcdr.
+  have w_nz: w != 0 := dvdpN0 w_r r_nz.
+  have p_m: p %| m  by rewrite dvdp_gcd vq_up -mulNr !dvdp_mull.
+  rewrite (leq_trans _ (dvdp_leq r_nz w_r)) // -(ltn_add2l (size v)).
+  rewrite addnC -ltn_subRL subn1 -size_mul // mulrC -wv size_scale //.
+  rewrite (leq_trans lt_vp) // dvdp_leq // -size_poly_eq0.
+  by rewrite -(size_scale _ nz_c) size_poly_eq0 wv mulf_neq0.
+have [[c p'] /= nz_c p'r] := dvdpP _ _ r_p.
+have [[k q'] /= nz_k q'r] := dvdpP _ _ r_q.
+have def_r := subnKC r_nonC; have r_nz: r != 0 by rewrite -size_poly_eq0 -def_r.
+have le_p'_dp: size p' <= dp.
+  have [-> | nz_p'] := eqVneq p' 0; first by rewrite size_poly0.
+  by rewrite /dp -(size_scale p nz_c) p'r size_mul // addnC -def_r leq_addl.
+have le_q'_dq: size q' <= dq.
+  have [-> | nz_q'] := eqVneq q' 0; first by rewrite size_poly0.
+  by rewrite /dq -(size_scale q nz_k) q'r size_mul // addnC -def_r leq_addl.
+exists (row_mx (- c *: poly_rV q') (k *: poly_rV p')).
+  apply: contraNneq r_nz; rewrite -row_mx0; case/eq_row_mx=> q0 p0.
+  have{p0} p0: p = 0.
+    apply/eqP; rewrite -size_poly_eq0 -(size_scale p nz_c) p'r.
+    rewrite -(size_scale _ nz_k) scalerAl -(poly_rV_K le_p'_dp) -linearZ p0.
+    by rewrite linear0 mul0r size_poly0.
+  rewrite /r p0 gcd0p -size_poly_eq0 -(size_scale q nz_k) q'r.
+  rewrite -(size_scale _ nz_c) scalerAl -(poly_rV_K le_q'_dq) -linearZ.
+  by rewrite -[c]opprK scaleNr q0 !linear0 mul0r size_poly0.
+rewrite mul_row_col scaleNr mulNmx !mul_rV_lin1 /= !linearZ /= !poly_rV_K //.
+by rewrite !scalerCA p'r q'r mulrCA addNr.
+Qed.
+
+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 horner_mx_additive := [additive of horner_mx].
+Canonical horner_mx_rmorphism := [rmorphism of horner_mx].
+
+Lemma horner_mx_C a : horner_mx a%:P = a%:M.
+Proof. exact: horner_morphC. Qed.
+
+Lemma horner_mx_X : horner_mx 'X = A. Proof. exact: horner_morphX. Qed.
+
+Lemma horner_mxZ : scalable horner_mx.
+Proof.
+move=> a p /=; rewrite -mul_polyC rmorphM /=.
+by rewrite horner_mx_C [_ * _]mul_scalar_mx.
+Qed.
+
+Canonical horner_mx_linear := AddLinear horner_mxZ.
+Canonical horner_mx_lrmorphism := [lrmorphism of 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).
+Proof.
+rewrite mulmx_sum_row linear_sum [rVpoly u]poly_def rmorph_sum.
+apply: eq_bigr => i _.
+by rewrite valK !linearZ rmorphX /= horner_mx_X rowK /= mxvecK.
+Qed.
+
+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 diagA := [seq A i i | i : 'I_n].
+Let size_diagA : size diagA = n.
+Proof. by rewrite size_image card_ord. Qed.
+
+Let split_diagA :
+  exists2 q, \prod_(x <- diagA) ('X - x%:P) + q = char_poly & size q <= n.-1.
+Proof.
+rewrite [char_poly](bigD1 1%g) //=; set q := \sum_(s | _) _; exists q.
+  congr (_ + _); rewrite odd_perm1 mul1r big_map enumT; apply: eq_bigr => i _.
+  by rewrite !mxE perm1 eqxx.
+apply: leq_trans {q}(size_sum _ _ _) _; apply/bigmax_leqP=> s nt_s.
+have{nt_s} [i nfix_i]: exists i, s i != i.
+  apply/existsP; rewrite -negb_forall; apply: contra nt_s => s_1.
+  by apply/eqP; apply/permP=> i; apply/eqP; rewrite perm1 (forallP s_1).
+apply: leq_trans (_ : #|[pred j | s j == j]|.+1 <= n.-1).
+  rewrite -sum1_card (@big_mkcond nat) /= size_Msign.
+  apply: (big_ind2 (fun p m => size p <= m.+1)) => [| p mp q mq IHp IHq | j _].
+  - by rewrite size_poly1.
+  - apply: leq_trans (size_mul_leq _ _) _.
+    by rewrite -subn1 -addnS leq_subLR addnA leq_add.
+  rewrite !mxE eq_sym !inE; case: (s j == j); first by rewrite polyseqXsubC. 
+  by rewrite sub0r size_opp size_polyC leq_b1.
+rewrite -{8}[n]card_ord -(cardC (pred2 (s i) i)) card2 nfix_i !ltnS.
+apply: subset_leq_card; apply/subsetP=> j; move/(_ =P j)=> fix_j.
+rewrite !inE -{1}fix_j (inj_eq (@perm_inj _ s)) orbb.
+by apply: contraNneq nfix_i => <-; rewrite fix_j.
+Qed.   
+
+Lemma size_char_poly : size char_poly = n.+1.
+Proof.
+have [q <- lt_q_n] := split_diagA; have le_q_n := leq_trans lt_q_n (leq_pred n).
+by rewrite size_addl size_prod_XsubC size_diagA.
+Qed.
+
+Lemma char_poly_monic : char_poly \is monic.
+Proof.
+rewrite monicE -(monicP (monic_prod_XsubC diagA xpredT id)).
+rewrite !lead_coefE size_char_poly.
+have [q <- lt_q_n] := split_diagA; have le_q_n := leq_trans lt_q_n (leq_pred n).
+by rewrite size_prod_XsubC size_diagA coefD (nth_default 0 le_q_n) addr0.
+Qed.
+
+Lemma char_poly_trace : n > 0 -> char_poly`_n.-1 = - \tr A.
+Proof.
+move=> n_gt0; have [q <- lt_q_n] := split_diagA; set p := \prod_(x <- _) _.
+rewrite coefD {q lt_q_n}(nth_default 0 lt_q_n) addr0.
+have{n_gt0} ->: p`_n.-1 = ('X * p)`_n by rewrite coefXM eqn0Ngt n_gt0.
+have ->: \tr A = \sum_(x <- diagA) x by rewrite big_map enumT.
+rewrite -size_diagA {}/p; elim: diagA => [|x d IHd].
+  by rewrite !big_nil mulr1 coefX oppr0.
+rewrite !big_cons coefXM mulrBl coefB IHd opprD addrC; congr (- _ + _).
+rewrite mul_polyC coefZ [size _]/= -(size_prod_XsubC _ id) -lead_coefE. 
+by rewrite (monicP _) ?monic_prod_XsubC ?mulr1.
+Qed.
+
+Lemma char_poly_det : char_poly`_0 = (- 1) ^+ n * \det A.
+Proof.
+rewrite big_distrr coef_sum [0%N]lock /=; apply: eq_bigr => s _.
+rewrite -{1}rmorphN -rmorphX mul_polyC coefZ /=.
+rewrite mulrA -exprD addnC exprD -mulrA -lock; congr (_ * _).
+transitivity (\prod_(i < n) - A i (s i)); last by rewrite prodrN card_ord.
+elim: (index_enum _) => [|i e IHe]; rewrite !(big_nil, big_cons) ?coef1 //.
+by rewrite coefM big_ord1 IHe !mxE coefB coefC coefMn coefX mul0rn sub0r.
+Qed.
+
+End CharPoly.
+
+Lemma mx_poly_ring_isom (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].
+Proof.
+set M_RX := 'M[{poly R}]_n; set MR_X := ({poly 'M[R]_n}).
+pose Msize (A : M_RX) := \max_i \max_j size (A i j).
+pose phi (A : M_RX) := \poly_(k < Msize A) \matrix_(i, j) (A i j)`_k.
+have coef_phi A i j k: (phi A)`_k i j = (A i j)`_k.
+  rewrite coef_poly; case: (ltnP k _) => le_m_k; rewrite mxE // nth_default //.
+  apply: leq_trans (leq_trans (leq_bigmax i) le_m_k); exact: (leq_bigmax j).
+have phi_is_rmorphism : rmorphism phi.
+  do 2?[split=> [A B|]]; apply/polyP=> k; apply/matrixP=> i j; last 1 first.
+  - rewrite coef_phi mxE coefMn !coefC.
+    by case: (k == _); rewrite ?mxE ?mul0rn.
+  - by rewrite !(coef_phi, mxE, coefD, coefN).
+  rewrite !coef_phi !mxE !coefM summxE coef_sum.
+  pose F k1 k2 := (A i k1)`_k2 * (B k1 j)`_(k - k2).
+  transitivity (\sum_k1 \sum_(k2 < k.+1) F k1 k2); rewrite {}/F.
+    by apply: eq_bigr=> k1 _; rewrite coefM.
+  rewrite exchange_big /=; apply: eq_bigr => k2 _.
+  by rewrite mxE; apply: eq_bigr => k1 _; rewrite !coef_phi.
+have bij_phi: bijective phi.
+  exists (fun P : MR_X => \matrix_(i, j) \poly_(k < size P) P`_k i j) => [A|P].
+    apply/matrixP=> i j; rewrite mxE; apply/polyP=> k.
+    rewrite coef_poly -coef_phi.
+    by case: leqP => // P_le_k; rewrite nth_default ?mxE.
+  apply/polyP=> k; apply/matrixP=> i j; rewrite coef_phi mxE coef_poly.
+  by case: leqP => // P_le_k; rewrite nth_default ?mxE.
+exists (RMorphism phi_is_rmorphism).
+split=> // [p | A]; apply/polyP=> k; apply/matrixP=> i j.
+  by rewrite coef_phi coef_map !mxE coefMn.
+by rewrite coef_phi !mxE !coefC; case k; last rewrite /= mxE.
+Qed.
+
+Theorem Cayley_Hamilton (R : comRingType) n' (A : 'M[R]_n'.+1) :
+  horner_mx A (char_poly A) = 0.
+Proof.
+have [phi [_ phiZ phiC _]] := mx_poly_ring_isom R n'.
+apply/rootP/factor_theorem; rewrite -phiZ -mul_adj_mx rmorphM.
+by move: (phi _) => q; exists q; rewrite rmorphB phiC phiZ map_polyX.
+Qed.
+
+Lemma eigenvalue_root_char (F : fieldType) n (A : 'M[F]_n) a :
+  eigenvalue A a = root (char_poly A) a.
+Proof.
+transitivity (\det (a%:M - A) == 0).
+  apply/eigenvalueP/det0P=> [[v Av_av v_nz] | [v v_nz Av_av]]; exists v => //.
+    by rewrite mulmxBr Av_av mul_mx_scalar subrr.
+  by apply/eqP; rewrite -mul_mx_scalar eq_sym -subr_eq0 -mulmxBr Av_av.
+congr (_ == 0); rewrite horner_sum; apply: eq_bigr => s _.
+rewrite hornerM horner_exp !hornerE; congr (_ * _).
+rewrite (big_morph _ (fun p q => hornerM p q a) (hornerC 1 a)).
+by apply: eq_bigr => i _; rewrite !mxE !(hornerE, hornerMn).
+Qed.
+
+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.
+Proof. by exists (n ^ 2)%N; rewrite rank_leq_col. Qed.
+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.
+Proof.
+rewrite /d; case: ex_minnP; case=> //; rewrite leqn0 mxrank_eq0; move/eqP.
+move/row_matrixP; move/(_ 0); move/eqP; rewrite rowK row0 mxvec_eq0.
+by rewrite -mxrank_eq0 mxrank1.
+Qed.
+
+Lemma minpoly_mx1 : (1%:M \in Ad)%MS.
+Proof.
+by apply: (eq_row_sub (Ordinal mxminpoly_nonconstant)); rewrite rowK.
+Qed.
+
+Lemma minpoly_mx_free : row_free Ad.
+Proof.
+have:= mxminpoly_nonconstant; rewrite /d; case: ex_minnP; case=> // d' _.
+move/(_ d'); move/implyP; rewrite ltnn implybF -ltnS ltn_neqAle.
+by rewrite rank_leq_row andbT negbK.
+Qed.
+
+Lemma horner_mx_mem p : (horner_mx A p \in Ad)%MS.
+Proof.
+elim/poly_ind: p => [|p a IHp]; first by rewrite rmorph0 // linear0 sub0mx.
+rewrite rmorphD rmorphM /= horner_mx_C horner_mx_X.
+rewrite addrC -scalemx1 linearP /= -(mul_vec_lin (mulmxr_linear _ A)).
+case/submxP: IHp => u ->{p}.
+have: (powers_mx A (1 + d) <= Ad)%MS.
+  rewrite -(geq_leqif (mxrank_leqif_sup _)).
+    by rewrite (eqnP minpoly_mx_free) /d; case: ex_minnP.
+  rewrite addnC; apply/row_subP=> i.
+  by apply: eq_row_sub (lshift 1 i) _; rewrite !rowK.
+apply: submx_trans; rewrite addmx_sub ?scalemx_sub //.
+  by apply: (eq_row_sub 0); rewrite rowK.
+rewrite -mulmxA mulmx_sub {u}//; apply/row_subP=> i.
+rewrite row_mul rowK mul_vec_lin /= mulmxE -exprSr.
+by apply: (eq_row_sub (rshift 1 i)); rewrite rowK.
+Qed.
+
+Definition mx_inv_horner B := rVpoly (mxvec B *m pinvmx Ad).
+
+Lemma mx_inv_horner0 :  mx_inv_horner 0 = 0.
+Proof. by rewrite /mx_inv_horner !(linear0, mul0mx). Qed.
+
+Lemma mx_inv_hornerK B : (B \in Ad)%MS -> horner_mx A (mx_inv_horner B) = B.
+Proof. by move=> sBAd; rewrite horner_rVpoly mulmxKpV ?mxvecK. Qed.
+
+Lemma minpoly_mxM B C : (B \in Ad -> C \in Ad -> B * C \in Ad)%MS.
+Proof.
+move=> AdB AdC; rewrite -(mx_inv_hornerK AdB) -(mx_inv_hornerK AdC).
+by rewrite -rmorphM ?horner_mx_mem.
+Qed.
+
+Lemma minpoly_mx_ring : mxring Ad.
+Proof.
+apply/andP; split; first by apply/mulsmx_subP; exact: minpoly_mxM.
+apply/mxring_idP; exists 1%:M; split=> *; rewrite ?mulmx1 ?mul1mx //.
+  by rewrite -mxrank_eq0 mxrank1.
+exact: minpoly_mx1.
+Qed.
+
+Definition mxminpoly := 'X^d - mx_inv_horner (A ^+ d).
+Local Notation p_A := mxminpoly.
+
+Lemma size_mxminpoly : size p_A = d.+1.
+Proof. by rewrite size_addl ?size_polyXn // size_opp ltnS size_poly. Qed.
+
+Lemma mxminpoly_monic : p_A \is monic.
+Proof.
+rewrite monicE /lead_coef size_mxminpoly coefB coefXn eqxx /=.
+by rewrite nth_default ?size_poly // subr0.
+Qed.
+
+Lemma size_mod_mxminpoly p : size (p %% p_A) <= d.
+Proof.
+by rewrite -ltnS -size_mxminpoly ltn_modp // -size_poly_eq0 size_mxminpoly.
+Qed.
+
+Lemma mx_root_minpoly : horner_mx A p_A = 0.
+Proof.
+rewrite rmorphB -{3}(horner_mx_X A) -rmorphX /=.
+by rewrite mx_inv_hornerK ?subrr ?horner_mx_mem.
+Qed.
+
+Lemma horner_rVpolyK (u : 'rV_d) :
+  mx_inv_horner (horner_mx A (rVpoly u)) = rVpoly u.
+Proof.
+congr rVpoly; rewrite horner_rVpoly vec_mxK.
+by apply: (row_free_inj minpoly_mx_free); rewrite mulmxKpV ?submxMl.
+Qed.
+
+Lemma horner_mxK p : mx_inv_horner (horner_mx A p) = p %% p_A.
+Proof.
+rewrite {1}(Pdiv.IdomainMonic.divp_eq mxminpoly_monic p) rmorphD rmorphM /=.
+rewrite mx_root_minpoly mulr0 add0r.
+by rewrite -(poly_rV_K (size_mod_mxminpoly _)) horner_rVpolyK.
+Qed.
+
+Lemma mxminpoly_min p : horner_mx A p = 0 -> p_A %| p.
+Proof. by move=> pA0; rewrite /dvdp -horner_mxK pA0 mx_inv_horner0. Qed.
+
+Lemma horner_rVpoly_inj : @injective 'M_n 'rV_d (horner_mx A \o rVpoly).
+Proof.
+apply: can_inj (poly_rV \o mx_inv_horner) _ => u.
+by rewrite /= horner_rVpolyK rVpolyK.
+Qed.
+
+Lemma mxminpoly_linear_is_scalar : (d <= 1) = is_scalar_mx A.
+Proof.
+have scalP := has_non_scalar_mxP minpoly_mx1.
+rewrite leqNgt -(eqnP minpoly_mx_free); apply/scalP/idP=> [|[[B]]].
+  case scalA: (is_scalar_mx A); [by right | left].
+  by exists A; rewrite ?scalA // -{1}(horner_mx_X A) horner_mx_mem.
+move/mx_inv_hornerK=> <- nsB; case/is_scalar_mxP=> a defA; case/negP: nsB.
+move: {B}(_ B); apply: poly_ind => [|p c].
+  by rewrite rmorph0 ?mx0_is_scalar.
+rewrite rmorphD ?rmorphM /= horner_mx_X defA; case/is_scalar_mxP=> b ->.
+by rewrite -rmorphM horner_mx_C -rmorphD /= scalar_mx_is_scalar.
+Qed.
+
+Lemma mxminpoly_dvd_char : p_A %| char_poly A.
+Proof. by apply: mxminpoly_min; exact: Cayley_Hamilton. Qed.
+
+Lemma eigenvalue_root_min a : eigenvalue A a = root p_A a.
+Proof.
+apply/idP/idP=> Aa; last first.
+  rewrite eigenvalue_root_char !root_factor_theorem in Aa *.
+  exact: dvdp_trans Aa mxminpoly_dvd_char.
+have{Aa} [v Av_av v_nz] := eigenvalueP Aa.
+apply: contraR v_nz => pa_nz; rewrite -{pa_nz}(eqmx_eq0 (eqmx_scale _ pa_nz)).
+apply/eqP; rewrite -(mulmx0 _ v) -mx_root_minpoly.
+elim/poly_ind: p_A => [|p c IHp].
+  by rewrite rmorph0 horner0 scale0r mulmx0.
+rewrite !hornerE rmorphD rmorphM /= horner_mx_X horner_mx_C scalerDl.
+by rewrite -scalerA mulmxDr mul_mx_scalar mulmxA -IHp -scalemxAl Av_av.
+Qed.
+
+End MinPoly.
+
+(* Parametricity. *)
+Section MapRingMatrix.
+
+Variables (aR rR : ringType) (f : {rmorphism aR -> rR}).
+Local Notation "A ^f" := (map_mx (GRing.RMorphism.apply f) A) : ring_scope.
+Local Notation fp := (map_poly (GRing.RMorphism.apply f)).
+Variables (d n : nat) (A : 'M[aR]_n).
+
+Lemma map_rVpoly (u : 'rV_d) : fp (rVpoly u) = rVpoly u^f.
+Proof.
+apply/polyP=> k; rewrite coef_map !coef_rVpoly.
+by case: (insub k) => [i|]; rewrite  /=  ?rmorph0 // mxE.
+Qed.
+
+Lemma map_poly_rV p : (poly_rV p)^f = poly_rV (fp p) :> 'rV_d.
+Proof. by apply/rowP=> j; rewrite !mxE coef_map. Qed.
+
+Lemma map_char_poly_mx : map_mx fp (char_poly_mx A) = char_poly_mx A^f.
+Proof.
+rewrite raddfB /= map_scalar_mx /= map_polyX; congr (_ - _).
+by apply/matrixP=> i j; rewrite !mxE map_polyC.
+Qed.
+
+Lemma map_char_poly : fp (char_poly A) = char_poly A^f.
+Proof. by rewrite -det_map_mx map_char_poly_mx. Qed.
+
+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).
+Proof.
+move=> nz_fp nz_fq; rewrite /resultant /Sylvester_mx !size_map_poly_id0 //.
+rewrite -det_map_mx /= map_col_mx; congr (\det (col_mx _ _));
+  by apply: map_lin1_mx => v; rewrite map_poly_rV rmorphM /= map_rVpoly.
+Qed.
+
+End MapResultant.
+
+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 e : (powers_mx A e)^f = powers_mx A^f e.
+Proof. by apply/row_matrixP=> i; rewrite -map_row !rowK map_mxvec rmorphX. Qed.
+
+Lemma map_horner_mx p : (horner_mx A p)^f = horner_mx A^f (fp p).
+Proof.
+rewrite -[p](poly_rV_K (leqnn _)) map_rVpoly.
+by rewrite !horner_rVpoly map_vec_mx map_mxM map_powers_mx.
+Qed.
+
+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.
+Proof. by apply: eq_ex_minn => e; rewrite -map_powers_mx mxrank_map. Qed.
+
+Lemma mxminpoly_map : mxminpoly A^f = fp (mxminpoly A).
+Proof.
+rewrite rmorphB; congr (_ - _).
+  by rewrite /= map_polyXn degree_mxminpoly_map.
+rewrite degree_mxminpoly_map -rmorphX /=.
+apply/polyP=> i; rewrite coef_map //= !coef_rVpoly degree_mxminpoly_map.
+case/insub: i => [i|]; last by rewrite rmorph0.
+by rewrite -map_powers_mx -map_pinvmx // -map_mxvec -map_mxM // mxE.
+Qed.
+
+Lemma map_mx_inv_horner u : fp (mx_inv_horner A u) = mx_inv_horner A^f u^f.
+Proof.
+rewrite map_rVpoly map_mxM map_mxvec map_pinvmx map_powers_mx.
+by rewrite /mx_inv_horner degree_mxminpoly_map.
+Qed.
+
+End MapField.
+
+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 := forall 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).
+Proof. by case=> p; exists p; rewrite // map_poly_comp rmorph_root. Qed.
+
+Lemma integral_id x : integralOver RtoK (RtoK x).
+Proof. by exists ('X - x%:P); rewrite ?monicXsubC ?rmorph_root ?root_XsubC. Qed.
+
+Lemma integral_nat n : integralOver RtoK n%:R.
+Proof. by rewrite -(rmorph_nat RtoK); apply: integral_id. Qed.
+
+Lemma integral0 : integralOver RtoK 0. Proof. exact: (integral_nat 0). Qed.
+
+Lemma integral1 : integralOver RtoK 1. Proof. exact: (integral_nat 1). Qed.
+
+Lemma integral_poly (p : {poly K}) :
+  (forall i, integralOver RtoK p`_i) <-> {in p : seq K, integralRange RtoK}.
+Proof.
+split=> intRp => [_ /(nthP 0)[i _ <-] // | i]; rewrite -[p]coefK coef_poly.
+by case: ifP => [ltip | _]; [apply/intRp/mem_nth | apply: integral0].
+Qed.
+
+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].
+Proof.
+move=> mon_p pw0 intRp intRq.
+pose memR y := exists x, y = RtoK x.
+have memRid x: memR (RtoK x) by exists x.
+have memR_nat n: memR n%:R by rewrite -(rmorph_nat RtoK).
+have [memR0 memR1]: memR 0 * memR 1 := (memR_nat 0%N, memR_nat 1%N).
+have memRN1: memR (- 1) by exists (- 1); rewrite rmorphN1.
+pose rVin (E : K -> Prop) n (a : 'rV[K]_n) := forall i, E (a 0 i).
+pose pXin (E : K -> Prop) (r : {poly K}) := forall i, E r`_i.
+pose memM E n (X : 'rV_n) y := exists a, rVin E n a /\ y = (a *m X^T) 0 0.
+pose finM E S := exists n, exists X, forall y, memM E n X y <-> S y.
+have tensorM E n1 n2 X Y: finM E (memM (memM E n2 Y) n1 X).
+  exists (n1 * n2)%N, (mxvec (X^T *m Y)) => y.
+  split=> [[a [Ea Dy]] | [a1 [/fin_all_exists[a /all_and2[Ea Da1]] ->]]].
+    exists (Y *m (vec_mx a)^T); split=> [i|].
+      exists (row i (vec_mx a)); split=> [j|]; first by rewrite !mxE; apply: Ea.
+      by rewrite -row_mul -{1}[Y]trmxK -trmx_mul !mxE.
+    by rewrite -[Y]trmxK -!trmx_mul mulmxA -mxvec_dotmul trmx_mul trmxK vec_mxK.
+  exists (mxvec (\matrix_i a i)); split.
+    by case/mxvec_indexP=> i j; rewrite mxvecE mxE; apply: Ea.
+  rewrite -[mxvec _]trmxK -trmx_mul mxvec_dotmul -mulmxA trmx_mul !mxE.
+  apply: eq_bigr => i _; rewrite Da1 !mxE; congr (_ * _).
+  by apply: eq_bigr => j _; rewrite !mxE.
+suffices [m [X [[u [_ Du]] idealM]]]: exists m,
+  exists X, let M := memM memR m X in M 1 /\ forall y, M y -> M (q.[w] * y).
+- do [set M := memM _ m X; move: q.[w] => z] in idealM *.
+  have MX i: M (X 0 i).
+    by exists (delta_mx 0 i); split=> [j|]; rewrite -?rowE !mxE.
+  have /fin_all_exists[a /all_and2[Fa Da1]] i := idealM _ (MX i).
+  have /fin_all_exists[r Dr] i := fin_all_exists (Fa i).
+  pose A := \matrix_(i, j) r j i; pose B := z%:M - map_mx RtoK A.
+  have XB0: X *m B = 0.
+    apply/eqP; rewrite mulmxBr mul_mx_scalar subr_eq0; apply/eqP/rowP=> i.
+    by rewrite !mxE Da1 mxE; apply: eq_bigr=> j _; rewrite !mxE mulrC Dr.
+  exists (char_poly A); first exact: char_poly_monic.
+  have: (\det B *: (u *m X^T)) 0 0 == 0.
+    rewrite scalemxAr -linearZ -mul_mx_scalar -mul_mx_adj mulmxA XB0 /=.
+    by rewrite mul0mx trmx0 mulmx0 mxE.
+  rewrite mxE -Du mulr1 rootE -horner_evalE -!det_map_mx; congr (\det _ == 0).
+  rewrite !raddfB /= !map_scalar_mx /= map_polyX horner_evalE hornerX.
+  by apply/matrixP=> i j; rewrite !mxE map_polyC /horner_eval hornerC.
+pose gen1 x E y := exists2 r, pXin E r & y = r.[x]; pose gen := foldr gen1 memR.
+have gen1S (E : K -> Prop) x y: E 0 -> E y -> gen1 x E y.
+  by exists y%:P => [i|]; rewrite ?hornerC ?coefC //; case: ifP.
+have genR S y: memR y -> gen S y.
+  by elim: S => //= x S IH in y * => /IH; apply: gen1S; apply: IH.
+have gen0 := genR _ 0 memR0; have gen_1 := genR _ 1 memR1.
+have{gen1S} genS S y: y \in S -> gen S y.
+  elim: S => //= x S IH /predU1P[-> | /IH//]; last exact: gen1S.
+  by exists 'X => [i|]; rewrite ?hornerX // coefX; apply: genR.
+pose propD (R : K -> Prop) := forall x y, R x -> R y -> R (x + y).
+have memRD: propD memR.
+  by move=> _ _ [a ->] [b ->]; exists (a + b); rewrite rmorphD.
+have genD S: propD (gen S).
+  elim: S => //= x S IH _ _ [r1 Sr1 ->] [r2 Sr2 ->]; rewrite -hornerD.
+  by exists (r1 + r2) => // i; rewrite coefD; apply: IH.
+have gen_sum S := big_ind _ (gen0 S) (genD S).
+pose propM (R : K -> Prop) := forall x y, R x -> R y -> R (x * y).
+have memRM: propM memR.
+  by move=> _ _ [a ->] [b ->]; exists (a * b); rewrite rmorphM.
+have genM S: propM (gen S).
+  elim: S => //= x S IH _ _ [r1 Sr1 ->] [r2 Sr2 ->]; rewrite -hornerM.
+  by exists (r1 * r2) => // i; rewrite coefM; apply: gen_sum => j _; apply: IH.
+have gen_horner S r y: pXin (gen S) r -> gen S y -> gen S r.[y].
+  move=> Sq Sy; rewrite horner_coef; apply: gen_sum => [[i _] /= _].
+  by elim: {2}i => [|n IHn]; rewrite ?mulr1 // exprSr mulrA; apply: genM.
+pose S := w :: q ++ p; suffices [m [X defX]]: finM memR (gen S).
+  exists m, X => M; split=> [|y /defX Xy]; first exact/defX.
+  apply/defX/genM => //; apply: gen_horner => // [i|]; last exact/genS/mem_head.
+  rewrite -[q]coefK coef_poly; case: ifP => // lt_i_q.
+  by apply: genS; rewrite inE mem_cat mem_nth ?orbT.
+pose intR R y := exists r, [/\ r \is monic, root r y & pXin R r].
+pose fix genI s := if s is y :: s1 then intR (gen s1) y /\ genI s1 else True.
+have{mon_p pw0 intRp intRq}: genI S.
+  split; set S1 := _ ++ _; first exists p.
+    split=> // i; rewrite -[p]coefK coef_poly; case: ifP => // lt_i_p.
+    by apply: genS; rewrite mem_cat orbC mem_nth.
+  have: all (mem S1) S1 by exact/allP.
+  elim: {-1}S1 => //= y S2 IH /andP[S1y S12]; split; last exact: IH.
+  have{q S S1 IH S1y S12 intRp intRq} [q mon_q qx0]: integralOver RtoK y.
+    by move: S1y; rewrite mem_cat => /orP[]; [apply: intRq | apply: intRp].
+  exists (map_poly RtoK q); split=> // [|i]; first exact: monic_map.
+  by rewrite coef_map /=; apply: genR.
+elim: {w p q}S => /= [_|x S IH [[p [mon_p px0 Sp]] /IH{IH}[m2 [X2 defS]]]].
+  exists 1%N, 1 => y; split=> [[a [Fa ->]] | Fy].
+    by rewrite tr_scalar_mx mulmx1; apply: Fa.
+  by exists y%:M; split=> [i|]; rewrite 1?ord1 ?tr_scalar_mx ?mulmx1 mxE.
+pose m1 := (size p).-1; pose X1 := \row_(i < m1) x ^+ i.
+have [m [X defM]] := tensorM memR m1 m2 X1 X2; set M := memM _ _ _ in defM.
+exists m, X => y; rewrite -/M; split=> [/defM[a [M2a]] | [q Sq]] -> {y}.
+  exists (rVpoly a) => [i|].
+    by rewrite coef_rVpoly; case/insub: i => // i; apply/defS/M2a.
+  rewrite mxE (horner_coef_wide _ (size_poly _ _)) -/(rVpoly a).
+  by apply: eq_bigr => i _; rewrite coef_rVpoly_ord !mxE.
+have M_0: M 0 by exists 0; split=> [i|]; rewrite ?mul0mx mxE.
+have M_D: propD M.
+  move=> _ _ [a [Fa ->]] [b [Fb ->]]; exists (a + b).
+  by rewrite mulmxDl !mxE; split=> // i; rewrite mxE; apply: memRD.
+have{M_0 M_D} Msum := big_ind _ M_0 M_D.
+rewrite horner_coef; apply: (Msum) => i _; case: i q`_i {Sq}(Sq i) => /=.
+elim: {q}(size q) => // n IHn i i_le_n y Sy.
+have [i_lt_m1 | m1_le_i] := ltnP i m1.
+  apply/defM; exists (y *: delta_mx 0 (Ordinal i_lt_m1)); split=> [j|].
+    by apply/defS; rewrite !mxE /= mulr_natr; case: eqP.
+  by rewrite -scalemxAl -rowE !mxE.
+rewrite -(subnK m1_le_i) exprD -[x ^+ m1]subr0 -(rootP px0) horner_coef.
+rewrite polySpred ?monic_neq0 // -/m1 big_ord_recr /= -lead_coefE.
+rewrite opprD addrC (monicP mon_p) mul1r subrK !mulrN -mulNr !mulr_sumr.
+apply: Msum => j _; rewrite mulrA mulrACA -exprD; apply: IHn.
+  by rewrite -addnS addnC addnBA // leq_subLR leq_add.
+by rewrite -mulN1r; do 2!apply: (genM) => //; apply: genR. 
+Qed.
+
+Lemma integral_root_monic u p :
+    p \is monic -> root p u -> {in p : seq K, integralRange RtoK} -> 
+  integralOver RtoK u.
+Proof.
+move=> mon_p pu0 intRp; rewrite -[u]hornerX.
+apply: integral_horner_root mon_p pu0 intRp _.
+by apply/integral_poly => i; rewrite coefX; apply: integral_nat.
+Qed.
+
+Hint Resolve (integral0 RtoK) (integral1 RtoK) (@monicXsubC K).
+
+Let XsubC0 (u : K) : root ('X - u%:P) u. Proof. by rewrite root_XsubC. Qed.
+Let intR_XsubC u :
+  integralOver RtoK (- u) -> {in 'X - u%:P : seq K, integralRange RtoK}.
+Proof. by move=> intRu v; rewrite polyseqXsubC !inE => /pred2P[]->. Qed.
+
+Lemma integral_opp u : integralOver RtoK u -> integralOver RtoK (- u).
+Proof. by rewrite -{1}[u]opprK => /intR_XsubC/integral_root_monic; apply. Qed.
+
+Lemma integral_horner (p : {poly K}) u :
+    {in p : seq K, integralRange RtoK} -> integralOver RtoK u -> 
+  integralOver RtoK p.[u].
+Proof. by move=> ? /integral_opp/intR_XsubC/integral_horner_root; apply. Qed.
+
+Lemma integral_sub u v :
+  integralOver RtoK u -> integralOver RtoK v -> integralOver RtoK (u - v).
+Proof.
+move=> intRu /integral_opp/intR_XsubC/integral_horner/(_ intRu).
+by rewrite !hornerE.
+Qed.
+
+Lemma integral_add u v :
+  integralOver RtoK u -> integralOver RtoK v -> integralOver RtoK (u + v).
+Proof. by rewrite -{2}[v]opprK => intRu /integral_opp; apply: integral_sub. Qed.
+
+Lemma integral_mul u v :
+  integralOver RtoK u -> integralOver RtoK v -> integralOver RtoK (u * v).
+Proof.
+rewrite -{2}[v]hornerX -hornerZ => intRu; apply: integral_horner.
+by apply/integral_poly=> i; rewrite coefZ coefX mulr_natr mulrb; case: ifP.
+Qed.
+
+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.
+Proof.
+split=> [] [p p_nz pu0]; last by exists p; rewrite ?monic_neq0.
+exists (mk_mon p); first by rewrite monicE lead_coefZ mulVf ?lead_coef_eq0.
+by rewrite linearZ rootE hornerZ (rootP pu0) mulr0.
+Qed.
+
+Lemma algebraic_id a : algebraicOver FtoE (FtoE a).
+Proof. exact/integral_algebraic/integral_id. Qed.
+
+Lemma algebraic0 : algebraicOver FtoE 0.
+Proof. exact/integral_algebraic/integral0. Qed.
+
+Lemma algebraic1 : algebraicOver FtoE 1.
+Proof. exact/integral_algebraic/integral1. Qed.
+
+Lemma algebraic_opp x : algebraicOver FtoE x -> algebraicOver FtoE (- x).
+Proof. by move/integral_algebraic/integral_opp/integral_algebraic. Qed.
+
+Lemma algebraic_add x y :
+  algebraicOver FtoE x -> algebraicOver FtoE y -> algebraicOver FtoE (x + y).
+Proof.
+move/integral_algebraic=> intFx /integral_algebraic intFy.
+exact/integral_algebraic/integral_add.
+Qed.
+
+Lemma algebraic_sub x y :
+  algebraicOver FtoE x -> algebraicOver FtoE y -> algebraicOver FtoE (x - y).
+Proof. by move=> algFx /algebraic_opp; apply: algebraic_add. Qed.
+
+Lemma algebraic_mul x y :
+  algebraicOver FtoE x -> algebraicOver FtoE y -> algebraicOver FtoE (x * y).
+Proof.
+move/integral_algebraic=> intFx /integral_algebraic intFy.
+exact/integral_algebraic/integral_mul.
+Qed.
+
+Lemma algebraic_inv u : algebraicOver FtoE u -> algebraicOver FtoE u^-1.
+Proof.
+have [-> | /expf_neq0 nz_u_n] := eqVneq u 0; first by rewrite invr0.
+case=> p nz_p pu0; exists (Poly (rev p)).
+  apply/eqP=> /polyP/(_ 0%N); rewrite coef_Poly coef0 nth_rev ?size_poly_gt0 //.
+  by apply/eqP; rewrite subn1 lead_coef_eq0.
+apply/eqP/(mulfI (nz_u_n (size p).-1)); rewrite mulr0 -(rootP pu0).
+rewrite (@horner_coef_wide _ (size p)); last first.
+  by rewrite size_map_poly -(size_rev p) size_Poly.
+rewrite horner_coef mulr_sumr size_map_poly.
+rewrite [rhs in _ = rhs](reindex_inj rev_ord_inj) /=.
+apply: eq_bigr => i _; rewrite !coef_map coef_Poly nth_rev // mulrCA.
+by congr (_ * _); rewrite -{1}(subnKC (valP i)) addSn addnC exprD exprVn ?mulfK.
+Qed.
+
+Lemma algebraic_div x y :
+  algebraicOver FtoE x -> algebraicOver FtoE y -> algebraicOver FtoE (x / y).
+Proof. by move=> algFx /algebraic_inv; apply: algebraic_mul. Qed.
+
+Lemma integral_inv x : integralOver FtoE x -> integralOver FtoE x^-1.
+Proof. by move/integral_algebraic/algebraic_inv/integral_algebraic. Qed.
+
+Lemma integral_div x y :
+  integralOver FtoE x -> integralOver FtoE y -> integralOver FtoE (x / y).
+Proof. by move=> algFx /integral_inv; apply: integral_mul. Qed.
+
+Lemma integral_root p u :
+    p != 0 -> root p u -> {in p : seq E, integralRange FtoE} ->
+  integralOver FtoE u.
+Proof.
+move=> nz_p pu0 algFp.
+have mon_p1: mk_mon p \is monic.
+  by rewrite monicE lead_coefZ mulVf ?lead_coef_eq0.
+have p1u0: root (mk_mon p) u by rewrite rootE hornerZ (rootP pu0) mulr0.
+apply: integral_root_monic mon_p1 p1u0 _ => _ /(nthP 0)[i ltip <-].
+rewrite coefZ mulrC; rewrite size_scale ?invr_eq0 ?lead_coef_eq0 // in ltip.
+by apply: integral_div; apply/algFp/mem_nth; rewrite -?polySpred.
+Qed.
+
+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.
+
+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 e m n (A : 'M_(m, n)) : eval_mx e (mx_term A) = A.
+Proof. by apply/matrixP=> i j; rewrite !mxE. Qed.
+
+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.
+Proof.
+apply/matrixP=> i k; rewrite !mxE /= ((big_morph (eval e)) 0 +%R) //=.
+by apply: eq_bigr => j _; rewrite /= !mxE.
+Qed.
+
+Local Notation morphAnd f := ((big_morph f) 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 r m n (A : 'M_(m, n)) : qf_form (mxrank_form r A).
+Proof.
+by elim: m r n A => [|m IHm] r [|n] A //=; rewrite GRing.Pick_form_qf /=.
+Qed.
+
+Lemma eval_mxrank e r m n (A : 'M_(m, n)) :
+  qf_eval e (mxrank_form r A) = (\rank (eval_mx e A) == r).
+Proof.
+elim: m r n A => [|m IHm] r [|n] A /=; try by case r.
+rewrite GRing.eval_Pick /mxrank unlock /=; set pf := fun _ => _.
+rewrite -(@eq_pick _ pf) => [|k]; rewrite {}/pf ?mxE // eq_sym.
+case: pick => [[i j]|] //=; set B := _ - _; have:= mxrankE B.
+case: (Gaussian_elimination B) r => [[_ _] _] [|r] //= <-; rewrite {}IHm eqSS.
+by congr (\rank _ == r); apply/matrixP=> k l; rewrite !(mxE, big_ord1) !tpermR.
+Qed.
+
+Lemma eval_vec_mx e m n (u : 'rV_(m * n)) :
+  eval_mx e (vec_mx u) = vec_mx (eval_mx e u).
+Proof. by apply/matrixP=> i j; rewrite !mxE. Qed.
+
+Lemma eval_mxvec e m n (A : 'M_(m, n)) :
+  eval_mx e (mxvec A) = mxvec (eval_mx e A).
+Proof. by rewrite -{2}[A]mxvecK eval_vec_mx vec_mxK. Qed.
+
+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).
+Proof. by apply/matrixP=> i j; do 2![rewrite !mxE //; case: split => ?]. Qed.
+
+Lemma submx_form_qf : qf_form submx_form.
+Proof.
+by rewrite (morphAnd (@qf_form _)) ?big1 //= => r _; rewrite !mxrank_form_qf.
+Qed.
+
+Lemma eval_submx e : qf_eval e submx_form = (eval_mx e A <= eval_mx e B)%MS.
+Proof.
+rewrite (morphAnd (qf_eval e)) //= big_andE /=.
+apply/forallP/idP=> /= [|sAB d]; last first.
+  rewrite !eval_mxrank eval_col_mx -addsmxE; apply/implyP=> /eqP <-.
+  by rewrite mxrank_leqif_sup ?addsmxSr // addsmx_sub sAB /=.
+move/(_ (inord (\rank (eval_mx e (col_mx A B))))).
+rewrite inordK ?ltnS ?rank_leq_col // !eval_mxrank eqxx /= eval_col_mx.
+by rewrite -addsmxE mxrank_leqif_sup ?addsmxSr // addsmx_sub; case/andP.
+Qed.
+
+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.
+Proof. by rewrite tuple.size_tuple card_ord. Qed.
+
+Lemma nth_seq_of_rV x0 v (i : 'I_d) : nth x0 (seq_of_rV v) i = v 0 i.
+Proof. by rewrite nth_fgraph_ord ffunE. Qed.
+
+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.
+Proof.
+elim: e k => [|v e IHe] k; first by rewrite !nth_nil mxE.
+rewrite /row_env /= nth_cat size_seq_of_rV.
+case: k => [|k]; first by rewrite (valP i) nth_seq_of_rV.
+by rewrite mulSn -addnA -if_neg -leqNgt leq_addr addKn IHe.
+Qed.
+
+Lemma eval_row_var e k : eval_mx (row_env e) (row_var k) = e`_k :> 'rV_d.
+Proof. by apply/rowP=> i; rewrite !mxE /= nth_row_env. Qed.
+
+Definition Exists_row_form k (f : form) :=
+  foldr GRing.Exists f (codom (fun i : 'I_d => k * d + i)%N).
+
+Lemma Exists_rowP 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)).
+Proof.
+move=> d_gt0; pose i_ j := Ordinal (ltn_pmod j d_gt0).
+have d_eq j: (j = j %/ d * d + i_ j)%N := divn_eq j d.
+split=> [[v f_v] | ]; last case/GRing.foldExistsP=> e' ee' f_e'.
+  apply/GRing.foldExistsP; exists (row_env (set_nth 0 e k v)) => {f f_v}// j.
+  rewrite [j]d_eq !nth_row_env nth_set_nth /=; case: eqP => // ->.
+  by case/imageP; exists (i_ j).
+exists (\row_i e'`_(k * d + i)); apply: eq_holds f_e' => j /=.
+move/(_ j): ee'; rewrite [j]d_eq !nth_row_env nth_set_nth /=.
+case: eqP => [-> | ne_j_k -> //]; first by rewrite mxE.
+apply/mapP=> [[r lt_r_d]]; rewrite -d_eq => def_j; case: ne_j_k.
+by rewrite def_j divnMDl // divn_small ?addn0.
+Qed.
+
+End Env.
+
+End MatrixFormula.
+
+End MatrixFormula.