(Joint Center)Library qe_rcf_th

(* (c) Copyright Microsoft Corporation and Inria. All rights reserved. *)
Require Import ssreflect ssrfun ssrbool eqtype ssrnat seq choice path fintype.
Require Import div bigop ssralg poly polydiv ssrnum perm zmodp ssrint.
Require Import polyorder polyrcf interval matrix mxtens.

Set Implicit Arguments.

Import GRing.Theory Num.Theory Num.Def Pdiv.Ring Pdiv.ComRing.

Local Open Scope nat_scope.
Local Open Scope ring_scope.

Section extra.

Variable R : rcfType.
Implicit Types (p q : {poly R}).

Proof. move=> sq; rewrite comp_polyE; case hp: (size p) => [|n]. move/eqP: hp; rewrite size_poly_eq0 => /eqP ->. by rewrite !big_ord0 mulr1 lead_coef0. rewrite big_ord_recr /= addrC lead_coefDl. by rewrite lead_coefZ lead_coef_exp // !lead_coefE hp. rewrite (leq_ltn_trans (size_sum _ _ _)) // size_scale; last first. rewrite - [n]/(n.+1.-1) -hp -lead_coefE ?lead_coef_eq0 //. by rewrite -size_poly_eq0 hp. rewrite polySpred ?ltnS ?expf_eq0; last first. by rewrite andbC -size_poly_eq0 gtn_eqF // ltnW. apply/bigmax_leqP => i _; rewrite size_exp. have [->|/size_scale-> ] := eqVneq p`_i 0; first by rewrite scale0r size_poly0. by rewrite (leq_trans (size_exp_leq _ _)) // ltn_mul2l -subn1 subn_gt0 sq /=. Qed.

Lemma mul2n n : (2 × n = n + n)%N.
Lemma mul3n n : (3 × n = n + (n + n))%N.
Lemma exp3n n : (3 ^ n)%N = (3 ^ n).-1.+1.

Definition exp3S n : (3 ^ n.+1 = 3 ^ n + (3 ^ n + 3 ^ n))%N
  := etrans (expnS 3 n) (mul3n (3 ^ n)).

Lemma tens_I3_mx (cR : comRingType) m n (M : 'M[cR]_(m,n)) :
  1%:M ×t M = castmx (esym (mul3n _ ), esym (mul3n _ ))
               (block_mx M 0
                         0 (block_mx M 0
                                     0 M : 'M_(m+m,n+n)%N)).

Lemma mul_1tensmx (cR : comRingType) (m n p: nat)
  (e3n : (n + (n + n) = 3 × n)%N)
  (A B C : 'M[cR]_(m, n)) (M : 'M[cR]_(n, p)) :
  castmx (erefl _, e3n)
         (row_mx A (row_mx B C)) ×m (1%:M ×t M)
  = castmx (erefl _, esym (mul3n _))
         (row_mx (A ×m M) (row_mx (B ×m M) (C ×m M))).

:TODO: backport to polydiv
Lemma coprimep_rdiv_gcd p q : (p != 0) || (q != 0)
  coprimep (rdivp p (gcdp p q)) (rdivp q (gcdp p q)).

:TODO: generalize to non idomainTypes and backport to polydiv
Lemma rgcdp_eq0 p q : rgcdp p q == 0 = (p == 0) && (q == 0).

:TODO: : move in polyorder
Lemma mu_eq0 : p x, p != 0 → (\mu_x p == 0%N) = (~~ root p x).

Notation lcn_neq0 := lc_expn_rscalp_neq0.

:TODO: : move to polyorder
Lemma mu_mod p q x : (\mu_x p < \mu_x q)%N
 \mu_x (rmodp p q) = \mu_x p.

:TODO: : move to polyorder
Lemma mu_add p q x : p + q != 0 →
  (minn (\mu_x p) (\mu_x q) \mu_x (p + q)%R)%N .

:TODO: : move to polydiv
Lemma mu_mod_leq : p q x, ~~ (q %| p)
  (\mu_x q \mu_x p)%N
  (\mu_x q \mu_x (rmodp p q)%R)%N.

Lemma sgp_right0 : forall (x : R), sgp_right 0 x = 0. Proof. by move=> x; rewrite /sgp_right size_poly0. Qed.

End extra.

Section ctmat.

Variable R : numFieldType.

Definition ctmat1 := \matrix_(i < 3, j < 3)
  (nth [::] [:: [:: 1%:Z ; 1 ; 1 ]
              ; [:: -1 ; 1 ; 1 ]
              ; [:: 0 ; 0 ; 1 ] ] i)`_j.

Lemma det_ctmat1 : \det ctmat1 = 2.

Notation zmxR := ((map_mx ((intmul 1) : intR)) _ _).

Lemma ctmat1_unit : zmxR ctmat1 \in unitmx.

Definition ctmat n := (ctmat1 ^t n).

Lemma ctmat_unit : n, zmxR (ctmat n) \in unitmx.

Lemma ctmat1_blocks : ctmat1 = (block_mx
         1 (row_mx 1 1)
  (col_mx (-1) 0) (block_mx 1 1 0 1 : 'M_(1+1)%N)).

Lemma tvec_sub n : (3 × (3 ^ n).-1.+1 = 3 ^ (n.+1) )%N.

Lemma tens_ctmat1_mx n (M : 'M_n) :
  ctmat1 ×t M = castmx (esym (mul3n _ ), esym (mul3n _ ))
         (block_mx M (row_mx M M)
                   (col_mx (-M) 0) (block_mx M M
                                             0 M : 'M_(n+n)%N)).

Definition coefs n i :=
  [seq (castmx (erefl _, exp3n _) (invmx (zmxR (ctmat n)))) i ord0
     | i <- enum 'I__]`_i.

End ctmat.

Section QeRcfTh.

Variable R : rcfType.
Implicit Types a b : R.
Implicit Types p q : {poly R}.

Notation zmxR := ((map_mx ((intmul 1) : intR)) _ _).
Notation midf a b := ((a + b) / 2%:~R).

Constraints and Tarski queries

Local Notation sgp_is q s := (fun x ⇒ (sgr q.[x] == s)).

Definition constraints (z : seq R) (sq : seq {poly R}) (sigma : seq int) :=
  (\sum_(x <- z) \prod_(i < size sq) (sgz (sq`_i).[x] == sigma`_i))%N.

Definition taq (z : seq R) (q : {poly R}) : int := \sum_(x <- z) (sgz q.[x]).

Lemma taq_constraint1 z q :
  taq z q = (constraints z [::q] [::1])%:~R - (constraints z [::q] [::-1])%:~R.

Lemma taq_constraint0 z q :
  taq z 1 = (constraints z [::q] [:: 0])%:~R
          + (constraints z [::q] [:: 1])%:~R
          + (constraints z [::q] [::-1])%:~R.

Lemma taq_no_constraint z : taq z 1 = (constraints z [::] [::])%:~R.

Lemma taq_constraint2 z q :
  taq z (q ^+ 2) = (constraints z [::q] [:: 1])%:~R
                 + (constraints z [::q] [::-1])%:~R.

Fixpoint sg_tab n : seq (seq int) :=
  if n is m.+1
    then flatten (map (fun xmap (fun lx :: l) (sg_tab m)) [::1;-1;0])
    else [::[::]].

Lemma sg_tab_nil n : (sg_tab n == [::]) = false.

Lemma size_sg_tab n : size (sg_tab n) = (3 ^ n)%N.

Lemma size_sg_tab_neq0 n : size (sg_tab n) != 0%N.

Definition comb_exp (R : realDomainType) (s : R) :=
  match sgz s with Posz 1 ⇒ 1%N | Negz 0 ⇒ 2 | _ ⇒ 0%N end.

Definition poly_comb (sq : seq {poly R}) (sc : seq int) : {poly R} :=
  \prod_(i < size sq) ((sq`_i) ^+ (comb_exp sc`_i)).

Eval compute in sg_tab 4.

Definition cvec z sq := let sg_tab := sg_tab (size sq) in
  \row_(i < 3 ^ size sq) ((constraints z sq (nth [::] sg_tab i))%:~R : int).
Definition tvec z sq := let sg_tab := sg_tab (size sq) in
  \row_(i < 3 ^ size sq) (taq z (map (poly_comb sq) sg_tab)`_i).

Lemma tvec_cvec1 z q : tvec z [::q] = (cvec z [::q]) ×m ctmat1.

Lemma cvec_rec z q sq :
  cvec z (q :: sq) = castmx (erefl _, esym (exp3S _))
       (row_mx (cvec (filter (sgp_is q 1) z) (sq))
               (row_mx (cvec (filter (sgp_is q (-1)) z) (sq))
                       (cvec (filter (sgp_is q 0) z) (sq)))).

Lemma poly_comb_cons q sq s ss :
  poly_comb (q :: sq) (s :: ss) = (q ^ (comb_exp s)) × poly_comb sq ss.

Lemma comb_expE (rR : realDomainType):
  (comb_exp (1 : rR) = 1%N) × (comb_exp (-1 : rR) = 2%N) × (comb_exp (0 : rR) = 0%N).

Lemma tvec_rec z q sq :
  tvec z (q :: sq) =
  castmx (erefl _, esym (exp3S _)) (
    (row_mx (tvec (filter (sgp_is q 1) z) (sq))
      (row_mx (tvec (filter (sgp_is q (-1)) z) (sq))
        (tvec (filter (sgp_is q 0) z) (sq)))) ×m
    (castmx (mul3n _, mul3n _) (ctmat1 ×t 1%:M))).

Lemma tvec_cvec z sq :
  tvec z sq = (cvec z sq) ×m (ctmat (size sq)).

Lemma cvec_tvec z sq :
  zmxR (cvec z (sq)) = (zmxR (tvec z (sq))) ×m (invmx (zmxR (ctmat (size (sq))))).

Lemma constraints1_tvec : z sq,
  (constraints z (sq) (nseq (size (sq)) 1))%:~R = (castmx (erefl _, exp3n _)
    (zmxR (tvec z (sq)) ×m (invmx (zmxR (ctmat (size (sq))))))) ord0 ord0.

Cauchy Index, relation with Tarski query

Local Notation seq_mids a s b := (pairmap (fun x ymidf x y) a (rcons s b)).
Local Notation noroot p := ( x, ~~ root p x).
Notation lcn_neq0 := lc_expn_rscalp_neq0.

Definition jump q p x: int :=
  let non_null := (q != 0) && odd (\mu_x p - \mu_x q) in
  let sign := (sgp_right (q × p) x < 0)%R in
    (-1) ^+ sign *+ non_null.

Definition cindex (a b : R) (q p : {poly R}) : int :=
  \sum_(x <- roots p a b) jump q p x.

Definition cindexR q p := \sum_(x <- rootsR p) jump q p x.

Definition sjump p x : int :=
  ((-1) ^+ (sgp_right p x < 0)%R) *+ odd (\mu_x p).

Definition variation (x y : R) : int := (sgz y) × (x × y < 0).

Definition cross p a b := variation p.[a] p.[b].

Definition crossR p := variation (sgp_minfty p) (sgp_pinfty p).

Definition sum_var (s : seq R) := \sum_(n <- pairmap variation 0 s) n.

Lemma cindexEba a b : b a p q, cindex a b p q = 0.

Lemma jump0p q x : jump 0 q x = 0.

Lemma taq_cindex a b p q : taq (roots p a b) q = cindex a b (p^` × q) p.

Lemma sum_varP s : 0 \notin ssum_var s = variation (head 0 s) (last 0 s).

Lemma jump_coprime p q : p != 0 → coprimep p q
  → x, root p xjump q p x = sjump (q × p) x.

Lemma sjump_neigh a b p x : p != 0 →
  {in neighpl p a x & neighpr p x b,
    yl yr, sjump p x = cross p yl yr}.

Lemma jump_neigh a b p q x : q × p != 0 →
  {in neighpl (q × p) a x & neighpr (q × p) x b, yl yr,
   jump q p x = cross (q × p) yl yr *+ ((q != 0) && (\mu_x p > \mu_x q)%N)}.

Lemma jump_mul2l (p q r : {poly R}) :
  p != 0 → jump (p × q) (p × r) =1 jump q r.

Lemma jump_mul2r (p q r : {poly R}) :
  p != 0 → jump (q × p) (r × p) =1 jump q r.

Lemma jumppc p c x : jump p c%:P x = 0.

Lemma noroot_jump q p x : ~~ root p xjump q p x = 0.

Lemma jump_mulCp c p q x : jump (c *: p) q x = (sgz c) × jump p q x.

Lemma jump_pmulC c p q x : jump p (c *: q) x = (sgz c) × jump p q x.

Lemma jump_mod p q x :
  jump p q x = sgz (lead_coef q) ^+ (rscalp p q) × jump (rmodp p q) q x.

Lemma cindexRP q p a b :
  {in `]-oo, a], noroot p}{in `[b , +oo[, noroot p}
  cindex a b q p = cindexR q p.

Lemma cindex0p a b q : cindex a b 0 q = 0.

Lemma cindexR0p p : cindexR 0 p = 0.

Lemma cindexpC a b p c : cindex a b p (c%:P) = 0.

Lemma cindexRpC q c : cindexR q c%:P = 0.

Lemma cindex_mul2r a b p q r : r != 0 →
  cindex a b (p × r) (q × r) = cindex a b p q.

Lemma cindex_mulCp a b p q c :
  cindex a b (c *: p) q = (sgz c) × cindex a b p q.

Lemma cindex_pmulC a b p q c :
  cindex a b p (c *: q) = (sgz c) × cindex a b p q.

Lemma cindex_mod a b p q :
  cindex a b p q =
  (sgz (lead_coef q) ^+ rscalp p q) × cindex a b (rmodp p q) q.

Lemma variation0r b : variation 0 b = 0.

Lemma variationC a b : variation a b = - variation b a.

Lemma variationr0 a : variation a 0 = 0.

Lemma variation_pmull a b c : c > 0 → variation (a × c) (b) = variation a b.

Lemma variation_pmulr a b c : c > 0 → variation a (b × c) = variation a b.

Lemma congr_variation a b a' b' : sgz a = sgz a'sgz b = sgz b'
  variation a b = variation a' b'.

Lemma crossRP p a b :
  {in `]-oo, a], noroot p}{in `[b , +oo[, noroot p}
  cross p a b = crossR p.

Lemma noroot_cross p a b : a b
  {in `]a, b[, noroot p}cross p a b = 0.

Lemma cross_pmul p q a b : q.[a] > 0 → q.[b] > 0 →
  cross (p × q) a b = cross p a b.

Lemma cross0 a b : cross 0 a b = 0.

Lemma crossR0 : crossR 0 = 0.

Lemma cindex_seq_mids a b : a < b
   p q, p != 0 → q != 0 → coprimep p q
  cindex a b q p + cindex a b p q =
  sum_var (map (horner (p × q)) (seq_mids a (roots (p × q) a b) b)).

Lemma cindex_inv a b : a < b p q,
  ~~ root (p × q) a~~ root (p × q) b
  cindex a b q p + cindex a b p q = cross (p × q) a b.

Definition next_mod p q := - (lead_coef q ^+ rscalp p q) *: rmodp p q.

Lemma next_mod0p q : next_mod 0 q = 0.

Lemma cindex_rec a b : a < b p q,
  ~~ root (p × q) a~~ root (p × q) b
  cindex a b q p = cross (p × q) a b + cindex a b (next_mod p q) q.

Lemma cindexR_rec p q :
  cindexR q p = crossR (p × q) + cindexR (next_mod p q) q.

Computation of cindex through changes_mods

Definition mods p q :=
  let fix aux p q n :=
    if n is m.+1
      then if p == 0 then [::] else p :: (aux q (next_mod p q) m)
      else [::] in aux p q (maxn (size p) (size q).+1).

Lemma mods_rec p q : mods p q =
  if p == 0 then [::] else p :: (mods q (next_mod p q)).

Lemma mods_eq0 p q : (mods p q == [::]) = (p == 0).

Lemma neq0_mods_rec p q : p != 0 → mods p q = p :: mods q (next_mod p q).

Lemma mods0p q : mods 0 q = [::].

Lemma modsp0 p : mods p 0 = if p == 0 then [::] else [::p].

Fixpoint changes (s : seq R) : nat :=
  (if s is a :: q then (a × (head 0 q) < 0)%R + changes q else 0)%N.

Definition changes_pinfty (p : seq {poly R}) := changes (map lead_coef p).
Definition changes_minfty (p : seq {poly R}) :=
  changes (map (fun p : {poly _}(-1) ^+ (~~ odd (size p)) × lead_coef p) p).

Definition changes_poly (p : seq {poly R}) :=
  (changes_minfty p)%:Z - (changes_pinfty p)%:Z.
Definition changes_mods p q := changes_poly (mods p q).

Lemma changes_mods0p q : changes_mods 0 q = 0.

Lemma changes_modsp0 p : changes_mods p 0 = 0.

Lemma changes_mods_rec p q :
  changes_mods p q = crossR (p × q) + changes_mods q (next_mod p q).

Lemma changes_mods_cindex p q : changes_mods p q = cindexR q p.

Definition taqR p q := changes_mods p (p^` × q).

Lemma taq_taqR p q : taq (rootsR p) q = taqR p q.

Section ChangesItvMod_USELESS.
Not used anymore, but the content of this section is used in the LMCS 2012 paper and in Cyril's thesis

Definition changes_horner (p : seq {poly R}) x :=
  changes (map (fun pp.[x]) p).
Definition changes_itv_poly a b (p : seq {poly R}) :=
  (changes_horner p a)%:Z - (changes_horner p b)%:Z.

Definition changes_itv_mods a b p q := changes_itv_poly a b (mods p q).

Lemma changes_itv_mods0p a b q : changes_itv_mods a b 0 q = 0.

Lemma changes_itv_modsp0 a b p : changes_itv_mods a b p 0 = 0.

Lemma changes_itv_mods_rec a b : a < b p q,
  ~~ root (p × q) a~~ root (p × q) b
  changes_itv_mods a b p q = cross (p × q) a b
                          + changes_itv_mods a b q (next_mod p q).

Lemma changes_itv_mods_cindex a b : a < b p q,
  all (fun p~~ root p a) (mods p q) →
  all (fun p~~ root p b) (mods p q) →
  changes_itv_mods a b p q = cindex a b q p.

Definition taq_itv a b p q := changes_itv_mods a b p (p^` × q).

Lemma taq_taq_itv a b : a < b p q,
  all (fun pp.[a] != 0) (mods p (p^` × q)) →
  all (fun pp.[b] != 0) (mods p (p^` × q)) →
  taq (roots p a b) q = taq_itv a b p q.

End ChangesItvMod_USELESS.

Definition tvecR p sq := let sg_tab := sg_tab (size sq) in
  \row_(i < 3^size sq) (taqR p (map (poly_comb sq) sg_tab)`_i).

Lemma tvec_tvecR sq p : tvec (rootsR p) sq = tvecR p sq.

Lemma all_prodn_gt0 : (I : finType) r (P : pred I) (F : Inat),
  (\prod_(i <- r | P i) F i > 0)%N
   i, i \in rP i → (F i > 0)%N.

Definition taqsR p sq i : R :=
  (taqR p (map (poly_comb sq) (sg_tab (size sq)))`_i)%:~R.

Definition ccount_weak p sq : R :=
  let fix aux s (i : nat) := if i is i'.+1
    then aux (taqsR p sq i' × coefs R (size sq) i' + s) i'
    else s in aux 0 (3 ^ size sq)%N.

Lemma constraints1P (p : {poly R}) (sq : seq {poly R}) :
  (constraints (rootsR p) (sq) (nseq (size (sq)) 1))%:~R
  = ccount_weak p sq.

Lemma ccount_weakP p sq : p != 0 →
  reflect ( x, (p.[x] == 0) && \big[andb/true]_(q <- sq) (q.[x] > 0))
  (ccount_weak p sq > 0).

Lemma myprodf_eq0 (S : idomainType)(I : eqType) (r : seq I) P (F : IS) :
  reflect (exists2 i, ((i \in r) && (P i)) & (F i == 0))
  (\prod_(i <- r| P i) F i == 0).

Definition bounding_poly (sq : seq {poly R}) := (\prod_(q <- sq) q)^`.

Lemma bounding_polyP (sq : seq {poly R}) :
  [\/ \big[andb/true]_(q <- sq) (lead_coef q > 0),
      \big[andb/true]_(q <- sq) ((-1)^+(size q).-1 × (lead_coef q) > 0) |
    x,
   ((bounding_poly sq).[x] == 0) && \big[andb/true]_(q <- sq) (q.[x] > 0)]
     x, \big[andb/true]_(q <- sq) (q.[x] > 0).

Lemma size_prod_eq1 (sq : seq {poly R}) :
  reflect ( q, q \in sqsize q = 1%N) (size (\prod_(q0 <- sq) q0) == 1%N).

Definition ccount_gt0 (sp sq : seq {poly R}) :=
  let p := \big[@rgcdp _/0%R]_(p <- sp) p in
    if p != 0 then 0 < ccount_weak p sq
    else let bq := bounding_poly sq in
        [|| \big[andb/true]_(q <- sq)(lead_coef q > 0) ,
            \big[andb/true]_(q <- sq)((-1)^+(size q).-1 *(lead_coef q) > 0) |
         0 < ccount_weak bq sq].

Lemma ccount_gt0P (sp sq : seq {poly R}) :
  reflect ( x, \big[andb/true]_(p <- sp) (p.[x] == 0)
                  && \big[andb/true]_(q <- sq) (q.[x] > 0))
    (ccount_gt0 sp sq).

End QeRcfTh.