Library mathcomp.algebra.mxpoly
Require Import mathcomp.ssreflect.ssreflect.
From mathcomp
Require Import ssrfun ssrbool eqtype ssrnat seq div fintype tuple.
From mathcomp
Require Import finfun bigop fingroup perm ssralg zmodp matrix mxalgebra.
From mathcomp
Require Import poly polydiv.
Set Implicit Arguments.
Unset Strict Implicit.
Unset Printing Implicit Defensive.
Import GRing.Theory.
Import Monoid.Theory.
Local Open Scope ring_scope.
Import Pdiv.Idomain.
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.
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 //.
by apply: leq_trans (leq_trans (leq_bigmax i) le_m_k); apply: (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; apply: 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; apply: 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.
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 apply/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.
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.
From mathcomp
Require Import ssrfun ssrbool eqtype ssrnat seq div fintype tuple.
From mathcomp
Require Import finfun bigop fingroup perm ssralg zmodp matrix mxalgebra.
From mathcomp
Require Import poly polydiv.
Set Implicit Arguments.
Unset Strict Implicit.
Unset Printing Implicit Defensive.
Import GRing.Theory.
Import Monoid.Theory.
Local Open Scope ring_scope.
Import Pdiv.Idomain.
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.
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 //.
by apply: leq_trans (leq_trans (leq_bigmax i) le_m_k); apply: (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; apply: 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; apply: 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.
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 apply/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.
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.