Library mathcomp.algebra.mxpoly

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

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

This file provides basic support for formal computation with matrices, mainly results combining matrices and univariate polynomials, such as the Cayley-Hamilton theorem; it also contains an extension of the first order representation of algebra introduced in ssralg (GRing.term/formula). rVpoly v == the little-endian decoding of the row vector v as a polynomial p = \sum_i (v 0 i)%:P * 'X^i. poly_rV p == the partial inverse to rVpoly, for polynomials of degree less than d to 'rV_d (d is inferred from the context). Sylvester_mx p q == the Sylvester matrix of p and q. resultant p q == the resultant of p and q, i.e., \det (Sylvester_mx p q). horner_mx A == the morphism from {poly R} to 'M_n (n of the form n'.+1) mapping a (scalar) polynomial p to the value of its scalar matrix interpretation at A (this is an instance of the generic horner_morph construct defined in poly). powers_mx A d == the d x (n ^ 2) matrix whose rows are the mxvec encodings of the first d powers of A (n of the form n'.+1). Thus, vec_mx (v *m powers_mx A d) = horner_mx A (rVpoly v). char_poly A == the characteristic polynomial of A. char_poly_mx A == a matrix whose determinant is char_poly A. companionmx p == a matrix whose char_poly is p mxminpoly A == the minimal polynomial of A, i.e., the smallest monic polynomial that annihilates A (A must be nontrivial). degree_mxminpoly A == the (positive) degree of mxminpoly A. mx_inv_horner A == the inverse of horner_mx A for polynomials of degree smaller than degree_mxminpoly A. 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.

Local Open Scope ring_scope.

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

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

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

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

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

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

Lemma rVpolyK : cancel rVpoly poly_rV.

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

Lemma poly_rV_is_linear : linear poly_rV.
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.


Section Resultant.

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

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

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

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

Definition resultant := \det Sylvester_mx.

End Resultant.


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

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

Section HornerMx.

Variables (R : comRingType) (n' : nat).
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 <- index_enum _ & true].
Let size_diagA : size diagA = n.

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

Lemma size_char_poly : size char_poly = n.+1.

Lemma char_poly_monic : char_poly \is monic.

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

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

End CharPoly.


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

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

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

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

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

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

Section MinPoly.

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

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

Lemma mxminpoly_nonconstant : d > 0.

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

Lemma minpoly_mx_free : row_free Ad.

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

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

Lemma mx_inv_horner0 : mx_inv_horner 0 = 0.

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

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

Lemma minpoly_mx_ring : mxring Ad.

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

Lemma size_mxminpoly : size p_A = d.+1.

Lemma mxminpoly_monic : p_A \is monic.

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

Lemma mx_root_minpoly : horner_mx A p_A = 0.

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

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

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

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

Lemma mxminpoly_linear_is_scalar : (d 1) = is_scalar_mx A.

Lemma mxminpoly_dvd_char : p_A %| char_poly A.

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

End MinPoly.



Parametricity.
Section MapRingMatrix.

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

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

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

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

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

End MapRingMatrix.

Section MapResultant.

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

End MapResultant.

Section MapComRing.

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

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

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

End MapComRing.

Section MapField.

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

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

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

Lemma degree_mxminpoly_map : degree_mxminpoly A^f = degree_mxminpoly A.

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

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

End MapField.

Section IntegralOverRing.

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

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

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

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

Lemma integral_id x : integralOver RtoK (RtoK x).

Lemma integral_nat n : integralOver RtoK n%:R.

Lemma integral0 : integralOver RtoK 0.

Lemma integral1 : integralOver RtoK 1.

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

End IntegralOverRing.

Section IntegralOverComRing.

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

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

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

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

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

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

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

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

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

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

End IntegralOverComRing.

Section IntegralOverField.

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

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

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

Lemma integral_algebraic u : algebraicOver FtoE u integralOver FtoE u.

Lemma algebraic_id a : algebraicOver FtoE (FtoE a).

Lemma algebraic0 : algebraicOver FtoE 0.

Lemma algebraic1 : algebraicOver FtoE 1.

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

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

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

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

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

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

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

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

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

End IntegralOverField.

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

Section MatrixFormula.

Variable F : fieldType.


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

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

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

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

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


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

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

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

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

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

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

Section Subsetmx.

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

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

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

Lemma submx_form_qf : qf_form submx_form.

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

End Subsetmx.

Section Env.

Variable d : nat.

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

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

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

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

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

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

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

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

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

End Env.

End MatrixFormula.

End MatrixFormula.