(Joint Center)Library MathComp.mxpoly

(* (c) Copyright Microsoft Corporation and Inria.                       
 You may distribute this file under the terms of the CeCILL-B license *)

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.

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.

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

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

Lemma rVpolyK : cancel rVpoly poly_rV.

Lemma poly_rV_K p : size p drVpoly (poly_rV p) = p.

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

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

End RowPoly.

Implicit Arguments poly_rV [R d].

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 kS_ p k | inr kS_ q k end.

Definition resultant := \det Sylvester_mx.

End Resultant.

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

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

Section HornerMx.

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

Definition horner_mx := horner_morph (fun ascalar_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.

Lemma horner_mx_X : horner_mx 'X = A.

Lemma horner_mxZ : scalable horner_mx.

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

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.

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

Lemma size_char_poly : size char_poly = n.+1.

Lemma char_poly_monic : char_poly \is monic.

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

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

End CharPoly.

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

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

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

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 : d, \rank (powers_mx A d.+1) d.
Definition degree_mxminpoly := ex_minn degree_mxminpoly_proof.
Local Notation d := degree_mxminpoly.
Local Notation Ad := (powers_mx A d).

Lemma mxminpoly_nonconstant : d > 0.

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

Lemma minpoly_mx_free : row_free Ad.

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

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

Lemma mx_inv_horner0 : mx_inv_horner 0 = 0.

Lemma mx_inv_hornerK B : (B \in Ad)%MShorner_mx A (mx_inv_horner B) = B.

Lemma minpoly_mxM B C : (B \in AdC \in AdB × C \in Ad)%MS.

Lemma minpoly_mx_ring : mxring Ad.

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

Lemma size_mxminpoly : size p_A = d.+1.

Lemma mxminpoly_monic : p_A \is monic.

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

Lemma mx_root_minpoly : horner_mx A p_A = 0.

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

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

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

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

Lemma mxminpoly_linear_is_scalar : (d 1) = is_scalar_mx A.

Lemma mxminpoly_dvd_char : p_A %| char_poly A.

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

End MinPoly.

Parametricity.
Section MapRingMatrix.

Variables (aR rR : ringType) (f : {rmorphism aRrR}).
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.

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

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

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

End MapRingMatrix.

Section MapResultant.

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

End MapResultant.

Section MapComRing.

Variables (aR rR : comRingType) (f : {rmorphism aRrR}).
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.

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

End MapComRing.

Section MapField.

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

Lemma degree_mxminpoly_map : degree_mxminpoly A^f = degree_mxminpoly A.

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

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

End MapField.

Section IntegralOverRing.

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

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

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

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

Lemma integral_id x : integralOver RtoK (RtoK x).

Lemma integral_nat n : integralOver RtoK n%:R.

Lemma integral0 : integralOver RtoK 0.

Lemma integral1 : integralOver RtoK 1.

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

End IntegralOverRing.

Section IntegralOverComRing.

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

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

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

Hint Resolve (integral0 RtoK) (integral1 RtoK) (@monicXsubC K).

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

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

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

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

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

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

End IntegralOverComRing.

Section IntegralOverField.

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

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

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

Lemma integral_algebraic u : algebraicOver FtoE u integralOver FtoE u.

Lemma algebraic_id a : algebraicOver FtoE (FtoE a).

Lemma algebraic0 : algebraicOver FtoE 0.

Lemma algebraic1 : algebraicOver FtoE 1.

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

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

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

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

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

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

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

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

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

End IntegralOverField.

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

Section MatrixFormula.

Variable F : fieldType.

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

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

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

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

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

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

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'.+1fun A : 'M_(1 + m', 1 + n')
    let nzA k := A k.1 k.2 != 0 in
    let xSchur k := Schur (xrow k.1 0%R (xcol k.2 0%R A)) in
    let recf k := Bool (r > 0) mxrank_form r.-1 (xSchur k) in
    GRing.Pick nzA recf (Bool (r == 0%N))
  | _, _fun _Bool (r == 0%N)
  end%T.

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

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

Lemma eval_vec_mx e m n (u : 'rV_(m × n)) :
  eval_mx e (vec_mx u) = vec_mx (eval_mx e u).

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

Section Subsetmx.

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

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

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

Lemma submx_form_qf : qf_form submx_form.

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

End Subsetmx.

Section Env.

Variable d : nat.

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

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

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

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

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

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

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

Definition Exists_row_form k (f : form) :=
  foldr GRing.Exists f (codom (fun i : 'I_dk × d + i)%N).

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

End Env.

End MatrixFormula.

End MatrixFormula.