(* Check the positivity of a multivariate polynomial in a fixed bound box *)
Require Import Reals Ring_polynom Psatz.
Open Scope R_scope.

(* We need to re-use many functions from the implementation. *)
Ltac get_goal_rhs := match goal with |- 0 < ?a => a end.

(* This function takes an un-normalized polynomal expression
   of type Ring_polynom.PExpr and returns a normalized
   polynomial expression (sparse Horner form) of type Ring_polynom.Pol *)
Definition norm_subst :=
  norm_subst 0%Z 1%Z Z.add Z.mul Z.sub Z.opp Zeq_bool Z.quotrem.

Definition Rphi := Pphi 0 Rplus Rmult (gen_phiZ 0 1 Rplus Rmult Ropp).

(* creates a list of length (p-1), with only 0 coefficients. *)
Definition padding_list p :=
  List.tl (iter_pos _ (fun l => (Pc 0%Z::l)%list) Datatypes.nil p).

(* From a polynomial of type Pol, returns the list of coefficients. *)
Fixpoint pol_to_coef_list (p : Ring_polynom.Pol Z) :=
   match p with
      PX a k b =>
      (b :: padding_list k ++ pol_to_coef_list a)%list
    | Pc v as it =>
      (it::Datatypes.nil)%list
    | Pinj c v => 
      (mkPinj_pred c v::Datatypes.nil)%list
   end.

(* From a list of coefficients, constructs a polynomial *)
Fixpoint coef_list_to_pol (l : list (Ring_polynom.Pol Z)) :=
  match l with
     Datatypes.nil => Pc 0%Z
  |  (Pc a as it::Datatypes.nil)%list => it
  |  (Ring_polynom.Pinj j P::Datatypes.nil)%list => 
       Ring_polynom.mkPinj (j + 1) P
  |  (a::l')%list =>
     match coef_list_to_pol l' with
     | PX P' i (Pc 0%Z) => PX P' (i + 1) a
     | any => PX any 1 a
     end
  end.

Lemma nth_S :
   forall A (a b : A) n l, List.nth (S n) (b::l) a = List.nth n l a.
 reflexivity.
Qed.

(* In this tactic, we first find the list of free variables,
   then we transform the target into a PExpr (w), called p;
   then the list of coefficients l.  Then p' is the normalized
   polynomial with the reversed coefficients.  The substition of
   (x + 1) for x is done by evaluating the polynomial, using function
   Rphi, which return an expression of Type R, which we need to 
   transform back into a polynomial: mkPol (to obtain something
   of type PExpr), then norm_subst (to normalize), and then we
   have the list Mobius coefficients (same sign as Bernstein coeffiients)
   We then show an equality that makes it easy to prove that the
   expression is positive. *)
Ltac test_rhs RNG Ih rl :=
  let  mkFV := get_RingFV RNG in let mkPol := get_RingMeta RNG in
  let rhs := get_goal_rhs in let fv := mkFV rhs (@List.nil R) in
  let x := eval unfold List.hd in (List.hd 0 fv) in
  let p := mkPol rhs fv in
  set (dummy := (p, fv)).

Ltac test_tac := ring_lookup (PackRing test_rhs) [] (0 = 0 :> R).

Ltac normalize_rhs RNG Ih rl :=
  let  mkFV := get_RingFV RNG in let mkPol := get_RingMeta RNG in
  let rhs := get_goal_rhs in let fv := mkFV rhs (@List.nil R) in
  let x := eval unfold List.hd in (List.hd 0 fv) in
  let p := mkPol rhs fv in
  let l := constr:(pol_to_coef_list (norm_subst 0 Datatypes.nil p)) in
  let deg := eval lazy in (pred (List.length l)) in
  let p'' := eval compute in 
  (* this is where the translation happens *)
      (Rphi ((x + 1)::List.tl fv)%list
      (coef_list_to_pol (List.rev l))) in
  let w := mkPol p'' fv in
  let w' := eval compute in
    (List.map (Rphi (List.tl fv))
        (pol_to_coef_list (norm_subst 0 Datatypes.nil w))) in
  let w2 :=
      constr:(sum_f_R0
                 (fun i =>
                    List.nth i w' 0 * ((1 - x) ^ i * x ^ (deg - i))) deg) in
      (idtac "working on variable" x; replace rhs with w2;
        rewrite !tech5, !nth_S; unfold sum_f_R0, List.nth);
      [repeat (apply Rplus_lt_0_compat);
            (apply Rmult_lt_0_compat;
             [ match goal with id : _ < x < _ |- _ => try clear id x end
               |apply Rmult_lt_0_compat; apply pow_lt; psatzl R]) 
       | simpl minus; try ring ].

(* This tactic picks the first variable it finds.  The interval
  for this variable must be expressed by an hypothesis of the form
  id : 0 < x < 1 (but no error message if the condition is not satisfied.
  I have yet to find the way to choose this variable. *)
Ltac eliminate_one :=
   ring_lookup (PackRing normalize_rhs) [] (0 = 0 :> R).

(* This tactic transforms a problem on an arbitrary interval
  into an equivalent problem on the (0,1) interval. *)
Ltac reframe x :=
match goal with
   id : ?a < x < ?b |- 0 < _ => 
   let u := constr:((x - a) / (b - a)) in
   replace x with (u * (b - a) + a);
   [ cut (0 < u < 1);[generalize u; clear x id; intros x id
        | match goal with |- _ < ?z < _ => field_simplify z end; try psatzl R]|
  field]
end.

(*
Lemma first_ex : forall x y, 0 < x < 1 -> 2 < y < 3 ->
    0 < x * y + 1.
Proof.
intros.
eliminate_one.
 psatzl R.
psatzl R.
Qed.

Lemma second_ex : forall x y, 1 < x < 3 -> 3 <  y < 5 -> 
           0 <  x * y ^ 2 * (y - x)^3 + 2 * x^ 4 + 1/3.
intros;
reframe x; unfold Rdiv.
(* ring_lookup (PackRing test_rhs) [] (0 = 0 :> R). *)
Time eliminate_one; reframe y; eliminate_one; try psatzl R.
Qed.
*)