(Joint Center)Library closed_field

(* (c) Copyright Microsoft Corporation and Inria. All rights reserved. *)
Require Import ssreflect ssrfun ssrbool eqtype ssrnat seq.
Require Import bigop ssralg poly polydiv.

Set Implicit Arguments.

Import GRing.
Open Local Scope ring_scope.

Import Pdiv.Ring.
Import PreClosedField.

Section TermEqType.
Variable R : UnitRing.type.
Fixpoint term_eq (t t' : term R) := match t, t' with | Var x, Var y => x == y | Const r, Const s => r == s | NatConst n, NatConst m => n == m | Add t t', Add s s' => term_eq t s && term_eq t' s' | Opp t, Opp s => term_eq t s | NatMul t n, NatMul s m => term_eq t s && (n == m) | Mul t t', Mul s s' => term_eq t s && term_eq t' s' | Inv t, Inv s => term_eq t s | Exp t n, Exp s m => term_eq t s && (n == m) | _, _ => false end.
Lemma term_eq_axiom : Equality.axiom term_eq. Proof. elim; do ? [by move=> ? [ ] *; apply: (iffP idP)=> //=; [move/eqP->|case=> -> ]#].
  • move=> ? P ? P' [ ] /= *; apply: (iffP idP)=> //=. by case/andP; move/P->; move/P'->. by case=> <- <-; apply/andP; split; [apply/P|apply/P' ].
  • move=> ? P [ ] /= *; apply: (iffP idP)=> //=; first by move/P->. by case=> <-; apply/P.
  • move=> ? P ? [ ] /= *; apply: (iffP idP)=> //=. by case/andP; move/P->; move/eqP->. by case=> <- <-; apply/andP; split; do 1?apply/P.
  • move=> ? P ? P' [ ] /= *; apply: (iffP idP)=> //=. by case/andP; move/P->; move/P'->. by case=> <- <-; apply/andP; split; [apply/P|apply/P' ].
  • move=> ? P [ ] /= *; apply: (iffP idP)=> //=; first by move/P->. by case=> <-; apply/P.
  • move=> ? P ? [ ] /= *; apply: (iffP idP)=> //=. by case/andP; move/P->; move/eqP->. by case=> <- <-; apply/andP; split; do 1?apply/P.
Canonical Structure term_eqType := EqType (term R) (EqMixin term_eq_axiom).
End TermEqType.
Section ClosedFieldQE.

Variable F : Field.type.

Variable axiom : ClosedField.axiom F.

Notation fF := (formula F).
Notation qf f := (qf_form f && rformula f).

Definition ifF (th el f: fF) : fF := ((f /\ th) \/ ((~ f) /\ el))%T.
Lemma ifFP : forall th el f e, qf_eval e (ifF th el f) = (fun e f => if f then qf_eval e th else qf_eval e el) e (qf_eval e f). Proof. move=> th el f e; rewrite /ifF /=. case: (qf_eval e f); rewrite //=. by case: (qf_eval _ _). Qed. Lemma ifF_qf : forall th el f & qf th & qf el & qf f, qf (ifF th el f). Proof. by move=> ? ? ? /=; do ? [case/andP=> -> -> ]. Qed.
Definition polyF := seq (term F).
Fixpoint eval_poly (e:seq F) pf :=
  if pf is c::qf then (eval_poly e qf)*'X + (eval e c)%:P else 0.

Definition sizeT (k : nat -> fF) (p : polyF) := Pick (fun i : 'I_(size p) => nth 0 p i != 0 /\ \big[And/True]#_(j < size p | j > i) (nth 0 p j == 0))%T (fun i => k i.+1) (k 0%N).

Definition rpoly p := (all (@rterm F) p).

Fixpoint sizeT (k : natfF) (p:polyF) :=
  if p is c::q then
    sizeT (fun n
      if n is m.+1 then k m.+2
        else GRing.If (c == 0) (k 0%N) (k 1%N)) q
    else k O%N.

Lemma sizeTP : k,
   p e, qf_eval e (sizeT k p) = qf_eval e (k (size (eval_poly e p))).

Lemma sizeT_qf : k p, ( n, qf (k n))
  → rpoly pqf (sizeT k p).

Definition isnull (k : boolfF) (p: polyF) :=
  sizeT (fun nk (n == 0%N)) p.

Lemma isnullP : k,
   p e, qf_eval e (isnull k p) = qf_eval e (k (eval_poly e p == 0)).

Lemma isnull_qf : k p, ( b, qf (k b))
  → rpoly pqf (isnull k p).

Definition lt_sizeT (k : boolfF) (p q : polyF) : fF :=
  sizeT (fun nsizeT (fun mk (n<m)) q) p.

Definition lift (p : {poly F}) := let: q := p in map Const q.

Lemma eval_lift : e p, eval_poly e (lift p) = p.

Fixpoint lead_coefT (k : term FfF) p :=
  if p is c::q then
    lead_coefT (fun lGRing.If (l == 0) (k c) (k l)) q
    else k (Const 0).

Lemma lead_coefTP : k,
  ( x e, qf_eval e (k x) = qf_eval e (k (Const (eval e x))))
  → p e, qf_eval e (lead_coefT k p)
    = qf_eval e (k (Const (lead_coef (eval_poly e p)))).

Lemma lead_coefT_qf : k p, ( c, rterm cqf (k c))
  → rpoly pqf (lead_coefT k p).

Fixpoint amulXnT (a:term F) (n:nat) : polyF:=
  if n is n'.+1 then (Const 0)::(amulXnT a n') else [::a].

Lemma eval_amulXnT : a n e,
  eval_poly e (amulXnT a n) = (eval e a)%:P × 'X^n.

Lemma ramulXnT: a n, rterm arpoly (amulXnT a n).

Fixpoint sumpT (p q : polyF) :=
  if p is a::p' then
    if q is b::q' then (Add a b)::(sumpT p' q')
      else p
    else q.

Lemma eval_sumpT : p q e,
  eval_poly e (sumpT p q) = (eval_poly e p) + (eval_poly e q).

Lemma rsumpT: p q, rpoly prpoly qrpoly (sumpT p q).

Fixpoint mulpT (p q : polyF) :=
  if p is a::p' then sumpT (map (Mul a) q) (Const 0::(mulpT p' q)) else [::].

Lemma eval_mulpT : p q e,
  eval_poly e (mulpT p q) = (eval_poly e p) × (eval_poly e q).

Lemma rpoly_map_mul : t p, rterm trpoly (map (Mul t) p) = rpoly p.

Lemma rmulpT: p q, rpoly prpoly qrpoly (mulpT p q).

Definition opppT := map (Mul (@Const F (-1))).
Lemma eval_opppT : p e, eval_poly e (opppT p) = - eval_poly e p.

Definition natmulpT n := map (Mul (@NatConst F n)).
Lemma eval_natmulpT : p n e,
  eval_poly e (natmulpT n p) = (eval_poly e p) *+ n.

Fixpoint redivp_rec_loopT (q : polyF) sq cq (k : nat × polyF × polyFfF)
  (c : nat) (qq r : polyF) (n : nat) {struct n}:=
  sizeT (fun sr
    if sr < sq then k (c, qq, r) else
      lead_coefT (fun lr
        let m := amulXnT lr (sr - sq) in
        let qq1 := sumpT (mulpT qq [::cq]) m in
        let r1 := sumpT (mulpT r ([::cq])) (opppT (mulpT m q)) in
        if n is n1.+1 then redivp_rec_loopT q sq cq k c.+1 qq1 r1 n1
          else k (c.+1, qq1, r1)
      ) r
  ) r.

Fixpoint redivp_rec_loop (q : {poly F}) sq cq
   (k : nat) (qq r : {poly F})(n : nat) {struct n} :=
  (* (n : nat) (k : nat) (qq r : {poly F}) {struct n} := *)
    if size r < sq then (k, qq, r) else
    let m := (lead_coef r) *: 'X^(size r - sq) in
    (* let m := (lead_coef r):P + m in
    let r1 := r * cqN, ::Const 0, p) else
      sizeT (fun sq =>
        sizeT (fun sp =>
          lead_coefT (fun lq =>
            redivp_rec_loopT q sq lq k 0 ::Const 0 p sp
          ) q
        ) p
      ) q
  ) q.

Lemma redivp_rec_loopP : forall q c qq r n, redivp_rec q c qq r n 
    (* = redivp_rec_loop q (size q) (lead_coef q) n c qq r. *)
    = redivp_rec_loop q (size q) (lead_coef q) c qq r n.
move=> q c qq r n.
by elim: n c qq r => | n Pn c qq r //; rewrite /= Pn.

Lemma redivpTP : forall k,
  (forall c qq r e,  qf_eval e (k (c,qq,r)) 
    = qf_eval e (k (c, lift (eval_poly e qq), lift (eval_poly e r))))
  -> forall p q e (d := (redivp (eval_poly e p) (eval_poly e q))),
    qf_eval e (redivpT p k q) = qf_eval e (k (d.1.1, lift d.1.2, lift d.2)).
move=> k Pk.
move=> p q e /=.
rewrite isnullP unlock.
case q0 : (_==_); first by rewrite Pk /= mul0r add0r polyC0.
rewrite !sizeTP lead_coefTP /=; last by move=> *; rewrite !redivp_rec_loopTP.
rewrite redivp_rec_loopTP /=; last by move=> *; rewrite Pk.
rewrite mul0r add0r polyC0.
by rewrite redivp_rec_loopP.

Lemma redivpT_qf : forall p k q,
  (forall r, && rpoly r.1.2 & rpoly r.2 -> qf (k r))
  -> rpoly p -> rpoly q -> qf (redivpT p k q).
move=> p k q kP rp rq; rewrite /redivpT.
apply: isnull_qf=> // b.
case b; first by apply: kP=> /=.
apply: sizeT_qf => // sq.
apply: sizeT_qf=> // sp.
apply: lead_coefT_qf=> // lq rlq.
exact: redivp_rec_loopT_qf.

Definition rmodpT (p : polyF) (k:polyF -> fF) (q : polyF) : fF :=
  redivpT p (fun d => k d.2) q.
Definition rdivpT (p : polyF) (k:polyF -> fF) (q : polyF) : fF :=
  redivpT p (fun d => k d.1.2) q.
Definition rscalpT (p : polyF) (k: nat -> fF) (q : polyF) : fF :=
  redivpT p (fun d => k d.1.1) q.
Definition rdvdpT (p : polyF) (k:bool -> fF) (q : polyF) : fF :=
  rmodpT p (isnull k) q.

Fixpoint rgcdp_loop n (pp qq : {poly F}) {struct n} :=
  if rmodp pp qq == 0 then qq 
    else if n is n1.+1 then rgcdp_loop n1 qq (rmodp pp qq)
        else rmodp pp qq.

Fixpoint rgcdp_loopT pp k n qq {struct n} :=
  rmodpT pp (isnull 
    (fun b => if b 
      then (k qq) 
      else (if n is n1.+1 
        then rmodpT pp (rgcdp_loopT qq k n1) qq 
        else rmodpT pp k qq)
  ) qq.


Lemma rgcdp_loopP: forall k,
  (forall p e, qf_eval e (k p) = qf_eval e (k (lift (eval_poly e p))))
  -> forall n p q e, qf_eval e (rgcdp_loopT p k n q) = 
    qf_eval e (k (lift (rgcdp_loop n (eval_poly e p) (eval_poly e q)))).
move=> k Pk n p q e.
elim: n p q e => /=.
  move=> p q e.
  rewrite redivpTP; last by move=>*; rewrite !isnullP eval_lift.
  rewrite isnullP eval_lift.
  case: (_ == 0); first by rewrite Pk.
  by rewrite redivpTP; last by move=>*; rewrite Pk.
move=> m Pm p q e.
rewrite redivpTP; last by move=>*; rewrite !isnullP eval_lift.
rewrite isnullP eval_lift.
case: (_ == 0); first by rewrite Pk.
by rewrite redivpTP; move=>*; rewrite ?Pm !eval_lift.

Lemma rgcdp_loopT_qf : forall p k q n,
  (forall r, rpoly r -> qf (k r))
  -> rpoly p -> rpoly q -> qf (rgcdp_loopT p k n q).
move=> p k q n; move: p k q.
elim: n=> |n ihn p k q kP rp rq.
  apply: redivpT_qf=> // r; case/andP=> _ rr.
  apply: isnull_qf=> // []; first exact: kP.
  by apply: redivpT_qf=> // r'; case/andP=> _ rr'; apply: kP.
apply: redivpT_qf=> // r; case/andP=> _ rr.
apply: isnull_qf=> // []; first exact: kP.
by apply: redivpT_qf=> // r'; case/andP=> _ rr'; apply: ihn.

Definition rgcdpT (p:polyF) k (q:polyF) : fF :=
  let aux p1 k q1 := isnull 
    (fun b => if b 
      then (k q1) 
      else (sizeT (fun n => (rgcdp_loopT p1 k n q1)) p1)) p1
    in (lt_sizeT (fun b => if b then (aux q k p) else (aux p k q)) p q). 

Lemma rgcdpTP : forall k,
  (forall p e, qf_eval e (k p) = qf_eval e (k (lift (eval_poly e p))))
    -> forall p q e, qf_eval e (rgcdpT p k q) = qf_eval e (k (lift (rgcdp (eval_poly e p) (eval_poly e q)))).
move=> k Pk p q e.
rewrite /rgcdpT !sizeTP.
case lqp: (_ < _).  
  rewrite isnullP.
  case q0: (_ == _); first by rewrite Pk (eqP q0) rgcdp0.
  rewrite sizeTP rgcdp_loopP; first by rewrite /rgcdp lqp q0.
  by move=> e' p'; rewrite Pk.
rewrite isnullP.
case p0: (_ == _); first by rewrite Pk (eqP p0) rgcd0p.
rewrite sizeTP rgcdp_loopP; first by rewrite /rgcdp lqp p0.
by move=> e' q'; rewrite Pk.

Lemma rgcdpT_qf :  forall p k q,  (forall r, rpoly r -> qf (k r))
  -> rpoly p -> rpoly q -> qf (rgcdpT p k q).
move=> p k q kP rp rq.
apply: sizeT_qf=> // n; apply: sizeT_qf=> // m.
by case:(_ < _); apply: isnull_qf=> //; 
  case; do ?apply: kP=> //;
  apply: sizeT_qf=> // n';
  apply: rgcdp_loopT_qf.

Fixpoint rgcdpTs k (ps : seq polyF) : fF :=
  if ps is p::pr then rgcdpTs (rgcdpT p k) pr else k ::Const 0.

Lemma rgcdpTsP : forall k,
  (forall p e, qf_eval e (k p) = qf_eval e (k (lift (eval_poly e p))))
  -> forall ps e, qf_eval e (rgcdpTs k ps) = qf_eval e (k (lift (\big@rgcdp _/0%:P_(i <- ps)(eval_poly e i)))).
move=> k Pk ps e.
elim: ps k Pk; first by move=> p Pk; rewrite /= big_nil Pk /= mul0r add0r.
move=> p ps Pps /= k Pk /=.
rewrite big_cons Pps.
  by rewrite rgcdpTP // eval_lift.
by move=>  p' e'; rewrite !rgcdpTP; first by rewrite Pk !eval_lift .

Definition rseq_poly ps := all rpoly ps.

Lemma rgcdpTs_qf : forall k ps,  (forall r, rpoly r -> qf (k r))
  -> rseq_poly ps -> qf (rgcdpTs k ps).
move=> k p; elim: p k=> |c p ihp k kP rps=> /=; first exact: kP.
move: rps; case/andP=> rc rp.
by apply: ihp=> // r rr; apply: rgcdpT_qf.

Fixpoint rgdcop_recT (q: polyF) k  (p : polyF) n :=
  if n is m.+1 then
    rgcdpT p (sizeT (fun sd =>
      if sd == 1:R]) q.

Lemma rgdcop_recTP : forall k,
  (forall p e, qf_eval e (k p) = qf_eval e (k (lift (eval_poly e p))))
  -> forall p q n e, qf_eval e (rgdcop_recT p k q n) 
    = qf_eval e (k (lift (rgdcop_rec (eval_poly e p) (eval_poly e q) n))).
move=> k Pk p q n e.
elim: n k Pk p q e => |n Pn k Pk p q e /=.
  rewrite isnullP /=.
  by case: (_==_); rewrite Pk /= mul0r add0r ?(polyC0,polyC1).
rewrite rgcdpTP ?sizeTP ?eval_lift //.
  rewrite /rcoprimep; case se : (_==_); rewrite Pk //.
  do ?rewrite (rgcdpTP,Pn,eval_lift,redivpTP) | move × //=.
by do ?rewrite (sizeTP,eval_lift) | move × //=.

Lemma rgdcop_recT_qf :  forall p k q n,  (forall r, rpoly r -> qf (k r))
  -> rpoly p -> rpoly q -> qf (rgdcop_recT p k q n).
move=> p k q n; elim: n p k q=> |n ihn p k q kP rp rq /=.
apply: isnull_qf=> //; first by case; rewrite kP.
apply: rgcdpT_qf=> // g rg.
apply: sizeT_qf=> // n'.
case:(_ == _); first exact: kP.
apply: rgcdpT_qf=> // g' rg'.
apply: redivpT_qf=> // r; case/andP=> rr _.
exact: ihn.

Definition rgdcopT q k p := sizeT (rgdcop_recT q k p) p.

Lemma rgdcopTP : forall k,
  (forall p e, qf_eval e (k p) = qf_eval e (k (lift (eval_poly e p))))
    -> forall p q e, qf_eval e (rgdcopT p k q)
      = qf_eval e (k (lift (rgdcop (eval_poly e p) (eval_poly e q)))).
Proof. by move=> *; rewrite sizeTP rgdcop_recTP 1?Pk. Qed.

Lemma rgdcopT_qf : forall p k q, (forall r, rpoly r -> qf (k r))
  -> rpoly p -> rpoly q -> qf (rgdcopT p k q).
by move=> p k q kP rp rq; apply: sizeT_qf => // n; apply: rgdcop_recT_qf.

Definition ex_elim_seq (ps : seq polyF) (q : polyF) :=
  (rgcdpTs (rgdcopT q (sizeT (fun n => Bool (n != 1:P]_(p <- ps)(eval_poly e p)) in
      qf_eval e (ex_elim_seq ps q) = (size (rgdcop (eval_poly e q) gp) != 1:P]_(j <- ps) eval_poly e (abstrX i j)
    =  \big@rgcdp _/0%:P_(j <- (map (eval_poly e) (map (abstrX i) (ps)))) j.
  by rewrite !big_map.
rewrite -!map_comp.
  have aux I (l : seq I) (P : I -> {poly F}) :
    \big(@gcdp F)/0_(j <- l) P j :P]_(j <- map (eval_poly e \o abstrX i) ps) j == 0).
  rewrite (eqP g0) rgdcop0.
  case m0 : (_ == 0)=> //=; rewrite ?(size_poly1,size_poly0) //=.
    rewrite abstrX_bigmul eval_bigmul -bigmap_id in m0.
    constructor=> [x] // [] //.
    case=> _; move/holds_conjn=> hc; move/hc:rqs.
    by rewrite -root_bigmul //= (eqP m0) root0.
  constructor; move/negP:m0; move/negP=>m0.
  case: (closed_nonrootP axiom _ m0) => x {m0}.
  rewrite abstrX_bigmul eval_bigmul -bigmap_id.
  rewrite root_bigmul=> m0.
  exists x; do 2?constructor=> //.
    apply/holds_conj; rewrite //= -root_biggcd.
    by rewrite (eqp_root (aux _ _ _ )) (eqP g0) root0.
  by apply/holds_conjn.
apply:(iffP (closed_rootP axiom _)); case=> x Px; exists x; move:Px => //=.
  rewrite (eqp_root (eqp_rgdco_gdco _ _)).
  rewrite root_gdco ?g0 //.
  rewrite -(eqp_root (aux _ _ _ )) root_biggcd.
  rewrite abstrX_bigmul eval_bigmul -bigmap_id root_bigmul.
  case/andP=> psr qsr.
  do 2?constructor.
    by apply/holds_conj.
  by apply/holds_conjn.
rewrite (eqp_root (eqp_rgdco_gdco _ _)).
rewrite root_gdco ?g0 // -(eqp_root (aux _ _ _ )) root_biggcd.
rewrite abstrX_bigmul eval_bigmul -bigmap_id root_bigmul=> [] // [hps hqs].
apply/andP; constructor.
  by apply/holds_conj.
by apply/holds_conjn.

Lemma wf_ex_elim : GRing.wf_QE_proj ex_elim.
Proof. by move=> i bc /= rbc; apply: ex_elim_qf. Qed.

Definition closed_fields_QEMixin := 
  QEdecFieldMixin wf_ex_elim holds_ex_elim.

End ClosedFieldQE.