Library closedfields

Require Import ssreflect ssrfun ssrbool eqtype ssrnat seq.
Require Import bigops ssralg poly polydiv.

Import GRing.

Import Prenex Implicits.

Open Scope ring_scope.

Section SeqExtension.

Lemma all_map : forall (A R : Type) p (f: A -> R) s,
  all p (map f s) = all (p \o f) s.

End SeqExtension.

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.

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).
Lemma ifF_qf : forall th el f & qf th & qf el & qf f, qf (ifF th el f).

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 rpoly p := (all (@rterm F) p).

Fixpoint sizeT (k : nat -> fF) (p:polyF) :=
  if p is c::q then
    sizeT (fun n =>
      if n is m.+1 then k m.+2
        else ifF (k 0%N) (k 1%N) (Equal c (Const 0))) q
    else k O%N.

Lemma sizeTP : forall k,
  forall p e, qf_eval e (sizeT k p) = qf_eval e (k (size (eval_poly e p))).
Lemma sizeT_qf : forall k p, (forall n, qf (k n))
  -> rpoly p -> qf (sizeT k p).

Definition isnull (k : bool -> fF) (p: polyF) := sizeT (fun n => k (n == 0%N)) p.
Lemma isnullP : forall k,
  forall p e, qf_eval e (isnull k p) = qf_eval e (k (eval_poly e p == 0)).

Lemma isnull_qf : forall k p, (forall b, qf (k b))
  -> rpoly p -> qf (isnull k p).

Definition lt_sizeT (k : bool -> fF) (p q : polyF) : fF :=
  sizeT (fun n => sizeT (fun m => k (n<m)) q) p.

Definition lift (p : {poly F}) := let: q := p in map Const q.
Lemma eval_lift : forall e p, eval_poly e (lift p) = p.

Fixpoint lead_coefT (k : term F -> fF) p :=
  if p is c::q then
    lead_coefT (fun l =>
      ifF (k c) (k l) (Equal l (Const 0))
    ) q
    else k (Const 0).

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

Lemma lead_coefT_qf : forall k p, (forall c, rterm c -> qf (k c))
  -> rpoly p -> qf (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 : forall a n e,
  eval_poly e (amulXnT a n) = (eval e a)%:P * 'X^n.

Lemma ramulXnT: forall a n, rterm a -> rpoly (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 : forall p q e,
  eval_poly e (sumpT p q) = (eval_poly e p) + (eval_poly e q).

Lemma rsumpT: forall p q, rpoly p -> rpoly q -> rpoly (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 : forall p q e,
  eval_poly e (mulpT p q) = (eval_poly e p) * (eval_poly e q).

Lemma rpoly_map_mul : forall t p, rterm t -> rpoly (map (Mul t) p) = rpoly p.

Lemma rmulpT: forall p q, rpoly p -> rpoly q -> rpoly (mulpT p q).

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

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

Fixpoint edivp_rec_loopT (q : polyF) sq cq (k : term F * polyF * polyF -> fF)
  (c : term F) (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 c1 := Mul cq c 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 edivp_rec_loopT q sq cq k c1 qq1 r1 n1
          else k (c1, qq1, r1)
      ) r
  ) r.

Fixpoint edivp_rec_loop (q : {poly F}) sq cq
  (n : nat) (c : F) (qq r : {poly F}) {struct n} :=
    if size r < sq then (c, qq, r) else
    let m := (lead_coef r)%:P * 'X^(size r - sq) in
    let c1 := cq * c in
    let qq1 := qq * cq%:P + m in
    let r1 := r * cq%:P - m * q in
    if n is n1.+1 then edivp_rec_loop q sq cq n1 c1 qq1 r1 else (c1, qq1, r1).

Lemma edivp_rec_loopTP : forall k,
  (forall c qq r e, qf_eval e (k (c,qq,r))
    = qf_eval e (k (Const (eval e c), lift (eval_poly e qq), lift (eval_poly e r))))
  -> forall q sq cq c qq r n e
    (d := edivp_rec_loop (eval_poly e q) sq (eval e cq) n
      (eval e c) (eval_poly e qq) (eval_poly e r)),
    qf_eval e (edivp_rec_loopT q sq cq k c qq r n)
    = qf_eval e (k (Const d.1.1, lift d.1.2, lift d.2)).

Lemma edivp_rec_loopT_qf : forall q sq cq k c qq r n,
  (forall r, [&& rterm r.1.1, rpoly r.1.2 & rpoly r.2] -> qf (k r))
  -> rpoly q -> rterm cq -> rterm c -> rpoly qq -> rpoly r
    -> qf (edivp_rec_loopT q sq cq k c qq r n).

Definition edivpT (p : polyF) (k : term F * polyF * polyF -> fF) (q : polyF) : fF :=
  isnull (fun b =>
    if b then k (Const 1, [::Const 0], p) else
      sizeT (fun sq =>
        sizeT (fun sp =>
          lead_coefT (fun lq =>
            edivp_rec_loopT q sq lq k 1 [::Const 0] p sp
          ) q
        ) p
      ) q
  ) q.

Lemma edivp_rec_loopP : forall q c qq r n, edivp_rec q n c qq r
    = edivp_rec_loop q (size q) (lead_coef q) n c qq r.

Lemma edivpTP : forall k,
  (forall c qq r e, qf_eval e (k (c,qq,r))
    = qf_eval e (k (Const (eval e c), lift (eval_poly e qq), lift (eval_poly e r))))
  -> forall p q e (d := (edivp (eval_poly e p) (eval_poly e q))),
    qf_eval e (edivpT p k q) = qf_eval e (k (Const d.1.1, lift d.1.2, lift d.2)).

Lemma edivpT_qf : forall p k q,
  (forall r, [&& rterm r.1.1, rpoly r.1.2 & rpoly r.2] -> qf (k r))
  -> rpoly p -> rpoly q -> qf (edivpT p k q).

Definition modpT (p : polyF) (k:polyF -> fF) (q : polyF) : fF :=
  edivpT p (fun d => k d.2) q.
Definition divpT (p : polyF) (k:polyF -> fF) (q : polyF) : fF :=
  edivpT p (fun d => k d.1.2) q.
Definition scalpT (p : polyF) (k:term F -> fF) (q : polyF) : fF :=
  edivpT p (fun d => k d.1.1) q.
Definition dvdpT (p : polyF) (k:bool -> fF) (q : polyF) : fF :=
  modpT p (isnull k) q.

Fixpoint gcdp_loop n (pp qq : {poly F}) {struct n} :=
  if pp %% qq == 0 then qq
    else if n is n1.+1 then gcdp_loop n1 qq (pp %% qq)
        else pp %% qq.
Fixpoint gcdp_loopT pp k n qq {struct n} :=
  modpT pp (isnull
    (fun b => if b
      then (k qq)
      else (if n is n1.+1
        then modpT pp (gcdp_loopT qq k n1) qq
        else modpT pp k qq)
    )
  ) qq.

Lemma gcdp_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 (gcdp_loopT p k n q) =
    qf_eval e (k (lift (gcdp_loop n (eval_poly e p) (eval_poly e q)))).

Lemma gcdp_loopT_qf : forall p k q n,
  (forall r, rpoly r -> qf (k r))
  -> rpoly p -> rpoly q -> qf (gcdp_loopT p k n q).


Definition gcdpT (p:polyF) k (q:polyF) : fF :=
  let aux p1 k q1 := isnull
    (fun b => if b
      then (k q1)
      else (sizeT (fun n => (gcdp_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 gcdpTP : 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 (gcdpT p k q) = qf_eval e (k (lift (gcdp (eval_poly e p) (eval_poly e q)))).

Lemma gcdpT_qf : forall p k q, (forall r, rpoly r -> qf (k r))
  -> rpoly p -> rpoly q -> qf (gcdpT p k q).

Fixpoint gcdpTs k (ps : seq polyF) : fF :=
  if ps is p::pr then gcdpTs (gcdpT p k) pr else k [::Const 0].

Lemma gcdpTsP : 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 (gcdpTs k ps) = qf_eval e (k (lift (\big[@gcdp _/0%:P]_(i <- ps)(eval_poly e i)))).

Definition rseq_poly ps := all rpoly ps.

Lemma gcdpTs_qf : forall k ps, (forall r, rpoly r -> qf (k r))
  -> rseq_poly ps -> qf (gcdpTs k ps).

Fixpoint gdcop_recT (q: polyF) k (p : polyF) n :=
  if n is m.+1 then
    gcdpT p (sizeT (fun sd =>
      if sd == 1%N then k p
        else gcdpT p (divpT p (fun r => gdcop_recT q k r m)) q
    )) q
    else isnull (fun b => k [::Const b%:R]) q.
Lemma gdcop_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 (gdcop_recT p k q n)
    = qf_eval e (k (lift (gdcop_rec (eval_poly e p) (eval_poly e q) n))).

Lemma gdcop_recT_qf : forall p k q n, (forall r, rpoly r -> qf (k r))
  -> rpoly p -> rpoly q -> qf (gdcop_recT p k q n).

Definition gdcopT q k p := sizeT (gdcop_recT q k p) p.
Lemma gdcopTP : 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 (gdcopT p k q)
      = qf_eval e (k (lift (gdcop (eval_poly e p) (eval_poly e q)))).
Lemma gdcopT_qf : forall p k q, (forall r, rpoly r -> qf (k r))
  -> rpoly p -> rpoly q -> qf (gdcopT p k q).

Definition ex_elim_seq (ps : seq polyF) (q : polyF) :=
  (gcdpTs (gdcopT q (sizeT (fun n => Bool (n != 1%N)))) ps).
Lemma ex_elim_seqP :
  forall ps q e,
    let gp := (\big[@gcdp _/0%:P]_(p <- ps)(eval_poly e p)) in
      qf_eval e (ex_elim_seq ps q) = (size (gdcop (eval_poly e q) gp) != 1%N).

Lemma ex_elim_seq_qf : forall ps q, rseq_poly ps -> rpoly q
  -> qf (ex_elim_seq ps q).

Fixpoint abstrX (i : nat) (t : term F) :=
  match t with
    | (Var n) => if n == i then [::Const 0; Const 1] else [::t]
    | (Opp x) => opppT (abstrX i x)
    | (Add x y) => sumpT (abstrX i x) (abstrX i y)
    | (Mul x y) => mulpT (abstrX i x) (abstrX i y)
    | (NatMul x n) => natmulpT n (abstrX i x)
    | (Exp x n) => let ax := (abstrX i x) in
      iter n (mulpT ax) [::Const 1]
    | _ => [::t]
  end.

Lemma abstrXP : forall i t e x,
  rterm t -> (eval_poly e (abstrX i t)).[x] = eval (set_nth 0 e i x) t.

Lemma rabstrX : forall i t, rterm t -> rpoly (abstrX i t).

Implicit Types tx ty : term F.

Lemma abstrX_mulM : forall i, {morph abstrX i : x y / Mul x y >-> mulpT x y}.
Lemma abstrX1 : forall i, abstrX i (Const 1) = [::Const 1].

Lemma eval_poly_mulM : forall e, {morph eval_poly e : x y / mulpT x y >-> mul x y}.
Lemma eval_poly1 : forall e, eval_poly e [::Const 1] = 1.

Notation abstrX_bigmul := (big_morph _ (abstrX_mulM _) (abstrX1 _)).
Notation eval_bigmul := (big_morph _ (eval_poly_mulM _) (eval_poly1 _)).
Notation bigmap_id := (big_map _ (fun _ => true) id).

Lemma rseq_poly_map : forall x ts,
  all (@rterm _) ts -> rseq_poly (map (abstrX x) ts).

Definition ex_elim (x : nat) (pqs : seq (term F) * seq (term F)) :=
  ex_elim_seq (map (abstrX x) pqs.1)
  (abstrX x (\big[Mul/Const 1]_(q <- pqs.2) q)).
Lemma ex_elim_qf : forall x pqs,
  dnf_rterm pqs -> qf (ex_elim x pqs).

Lemma holds_conj : forall e i x ps, all (@rterm _) ps ->
  (holds (set_nth 0 e i x) (foldr (fun t : term F => And (t == 0)) True ps)
  <-> all (fun p => root p x) (map (eval_poly e \o abstrX i) ps)).

Lemma holds_conjn : forall e i x ps, all (@rterm _) ps ->
  (holds (set_nth 0 e i x) (foldr (fun t : term F => And (t != 0)) True ps)
  <-> all (fun p => ~~root p x) (map (eval_poly e \o abstrX i) ps)).

Lemma holds_ex_elim : QE.holds_proj_axiom ex_elim.

Lemma wf_ex_elim : QE.wf_proj_axiom ex_elim.

Definition closed_fields_QEMixin :=
  QE.Mixin wf_ex_elim holds_ex_elim.

End ClosedFieldQE.