(**************************************************************)
(*   Copyright Dominique Larchey-Wendling [*]                 *)
(*                                                            *)
(*                             [*] Affiliation LORIA -- CNRS  *)
(**************************************************************)
(*      This file is distributed under the terms of the       *)
(*         CeCILL v2 FREE SOFTWARE LICENSE AGREEMENT          *)
(**************************************************************)

(** ** Single Diophantine equations *)

Require Import List Arith Omega Nat.

Require Import utils_tac utils_list sums pos vec. 
Require Import dio_logic dio_elem.

Set Implicit Arguments.

Local Notation "âˆ‘" := (msum plus 0).

Section convexity.

  Let convex_1 x p : 2*(x*(x+p)) <= x*x+(x+p)*(x+p).
  Proof.
    rewrite mult_assoc.
    repeat rewrite Nat.mul_add_distr_r.
    repeat rewrite Nat.mul_add_distr_l.
    rewrite (mult_comm p x).
    repeat rewrite <- mult_assoc.
    generalize (x*x) (x*p) (p*p); intros; omega.
  Qed.

  Let convex_2 x p : 2*(x*(x+p)) = x*x+(x+p)*(x+p) -> p = 0.
  Proof.
    rewrite mult_assoc.
    intros H.
    cut (p*p = 0).
    { destruct p; simpl; auto; discriminate. }
    revert H.
    repeat rewrite Nat.mul_add_distr_r.
    repeat rewrite Nat.mul_add_distr_l.
    rewrite (mult_comm p x).
    repeat rewrite <- mult_assoc.
    generalize (x*x) (x*p) (p*p); intros; omega.
  Qed.

  Fact convex_le x y : 2*(x*y) <= x*x+y*y.
  Proof.
    destruct (le_lt_dec x y).
    + replace y with (x+(y-x)) by omega.
      apply convex_1.
    + rewrite (mult_comm x y), plus_comm.
      replace x with (y+(x-y)) by omega.
      apply convex_1.
  Qed.

  Fact convex_eq x y : 2*(x*y) = x*x+y*y -> x = y.
  Proof.
    destruct (le_lt_dec x y).
    + replace y with (x+(y-x)) by omega.
      intros H; apply convex_2 in H; omega.
    + rewrite (mult_comm x y), plus_comm.
      replace x with (y+(x-y)) by omega.
      intros H; apply convex_2 in H; omega.
  Qed.

  Let convex_3 a t x y : 0 < t -> a*x+(a+t)*y = a*y+(a+t)*x -> x = y.
  Proof.
    intros H.
    repeat rewrite Nat.mul_add_distr_r.
    intros H1.
    apply Nat.mul_cancel_l with t; omega.
  Qed.
   
  Fact convex_neq a b x y : a < b -> a*x+b*y = a*y+b*x -> x = y.
  Proof.
    intros H.
    replace b with (a+(b-a)) by omega.
    apply convex_3; omega.
  Qed.

  Hint Resolve convex_le.

  Fact convex_n_le n (f g : nat -> nat) :  âˆ‘ n (fun i => 2*(f i*g i)) 
                                        <= âˆ‘ n (fun i => f i*f i + g i*g i).
  Proof.
    revert f g; induction n as [ | n IHn ]; intros f g.
    + rewrite msum_0; auto.
    + do 2 rewrite msum_S.
      apply plus_le_compat; auto.
  Qed.

  Hint Resolve convex_n_le.

  Let nat_le_sum a b c d : a <= b -> c <= d -> a+c = b+d -> a = b /\ c = d.
  Proof. intros; omega. Qed.


  (* This one can be used to encode a list of equations
     x1 = y1 ... xn = yn into a single equation of size
     linear in n if all the xi, yi have a bounded size 

     Hence one can transform a system elementary diophantine
     equations into a single polynomial equation of linear
     size !! *)

  Fact convex_n_eq n (f g : nat -> nat) : âˆ‘ n (fun i => 2*(f i*g i)) 
                                        = âˆ‘ n (fun i => f i*f i + g i*g i)
                                    <-> forall i, i < n -> f i = g i.
  Proof.
    split.
    + revert f g; induction n as [ | n IHn ]; intros f g.
      * intros; omega.
      * do 2 rewrite msum_S; intros H.
        apply nat_le_sum in H; auto.
        destruct H as (H1 & H2).
        apply convex_eq in H1.
        specialize (IHn _ _ H2).
        intros [ | ] ?; auto; apply IHn; omega.
    + intros Hfg.
      apply msum_ext.
      intros i Hi; rewrite Hfg; auto; ring.
  Qed.

End convexity.

Section diophantine_polynomial.

  Variable (V P : Set).

  Inductive dio_polynomial : Set :=
    | dp_nat : nat -> dio_polynomial                  (* natural number constant *)
    | dp_var : V   -> dio_polynomial                  (* existentially quantified variable *)
    | dp_par : P   -> dio_polynomial                  (* parameter *)
    | dp_comp : dio_op -> dio_polynomial -> dio_polynomial -> dio_polynomial.

  Notation dp_add := (dp_comp do_add).
  Notation dp_mul := (dp_comp do_mul).

  Fixpoint dp_var_list p :=
    match p with
      | dp_nat _      => nil
      | dp_var v      => v::nil
      | dp_par _      => nil
      | dp_comp _ p q => dp_var_list p ++ dp_var_list q
    end.

  Fixpoint dp_par_list p :=
    match p with
      | dp_nat _      => nil
      | dp_var _      => nil
      | dp_par x      => x::nil
      | dp_comp _ p q => dp_par_list p ++ dp_par_list q
    end.

  (* Ï� Ïƒ Î½ Ï† *)

  Fixpoint dp_eval Ï† Î½ p := 
    match p with
      | dp_nat n => n
      | dp_var v => Ï† v
      | dp_par i => Î½ i
      | dp_comp do_add p q => dp_eval Ï† Î½ p + dp_eval Ï† Î½ q 
      | dp_comp do_mul p q => dp_eval Ï† Î½ p * dp_eval Ï† Î½ q 
    end.

  Fact dp_eval_ext Ï† Î½ Ï†' Î½' p :
        (forall v, In v (dp_var_list p) -> Ï† v = Ï†' v) 
     -> (forall i, In i (dp_par_list p) -> Î½ i = Î½' i) 
     -> dp_eval Ï† Î½ p = dp_eval Ï†' Î½' p.
  Proof.
    induction p as [ | | | [] p Hp q Hq ]; simpl; intros H1 H2; f_equal; auto;
      ((apply Hp || apply Hq); intros; [ apply H1 | apply H2 ]; apply in_or_app; auto).
  Qed.

  Fact dp_eval_fix_add Ï† Î½ p q : dp_eval Ï† Î½ (dp_add p q) = dp_eval Ï† Î½ p + dp_eval Ï† Î½ q.
  Proof. trivial. Qed.

  Fact dp_eval_fix_mul Ï† Î½ p q : dp_eval Ï† Î½ (dp_mul p q) = dp_eval Ï† Î½ p * dp_eval Ï† Î½ q.
  Proof. trivial. Qed.

  Fixpoint dp_size p :=
    match p with
      | dp_nat n => 1
      | dp_var v => 1
      | dp_par i => 1
      | dp_comp _ p q => 1 + dp_size p + dp_size q 
    end.

  Fact dp_size_fix_comp o p q : dp_size (dp_comp o p q) = 1 + dp_size p + dp_size q.
  Proof. auto. Qed.

  Definition dio_single := (dio_polynomial * dio_polynomial)%type.
  Definition dio_single_size (e : dio_single) := dp_size (fst e) + dp_size (snd e).

  Definition dio_single_pred e Î½ := exists Ï†, dp_eval Ï† Î½ (fst e) = dp_eval Ï† Î½ (snd e).

End diophantine_polynomial.

Arguments dp_nat {V P}.
Arguments dp_var {V P}.
Arguments dp_par {V P}.
Arguments dp_comp {V P}.

Notation dp_add := (dp_comp do_add).
Notation dp_mul := (dp_comp do_mul).

Section dio_elem_dio_poly.

  Let dp_2xy u v : dio_polynomial nat nat := dp_mul (dp_nat 2) (dp_mul u v).
  Let dp_x2y2 u v : dio_polynomial nat nat := dp_add (dp_mul u u) (dp_mul v v).

  Let dp_2xy_size u v : dp_size (dp_2xy u v) = 3+dp_size u+dp_size v.
  Proof. auto. Qed.

  Let dp_x2y2_size u v : dp_size (dp_x2y2 u v) = 3+2*dp_size u+2*dp_size v.
  Proof. simpl; omega. Qed.

  Let dp_common e : dio_polynomial nat nat :=
    match e with
      | dee_nat c      => dp_nat c
      | dee_var v      => dp_var v
      | dee_par p      => dp_par p
      | dee_comp o v w => dp_comp o (dp_var v) (dp_var w)
    end.

  Let dp_common_size e : dp_size (dp_common e) <= 3.
  Proof. destruct e as [ | | | [] ]; simpl; auto. Qed.

  Let dp_eval_common Ï† Î½ e : dp_eval Ï† Î½ (dp_common e) = dee_eval Ï† Î½ e.
  Proof. destruct e as [ | | | [] ]; auto. Qed.  

  Let dp_left (c : dio_constraint) := dp_2xy (dp_var (fst c)) (dp_common (snd c)).
  Let dp_right (c : dio_constraint) := dp_x2y2 (dp_var (fst c)) (dp_common (snd c)).
 
  Let dee2dp_1 l := fold_right dp_add (dp_nat 0) (map dp_left l).
  Let dee2dp_2 l := fold_right dp_add (dp_nat 0) (map dp_right l).

  Let dee2dp_1_size l : dp_size (dee2dp_1 l) <= 1+8*length l.
  Proof.
    induction l as [ | (x,e) l IHl ].
    + simpl; auto.
    + unfold dee2dp_1; simpl fold_right; fold (dee2dp_1 l).
      rewrite dp_size_fix_comp.
      unfold dp_left; rewrite dp_2xy_size.
      unfold fst, snd.
      generalize (dp_common_size e); intros.
      simpl dp_size at 1; simpl length.
      rewrite Nat.mul_succ_r; omega.
  Qed.

  Let dee2dp_2_size l : dp_size (dee2dp_2 l) <= 1+12*length l.
  Proof.
    induction l as [ | (x,e) l IHl ].
    + simpl; auto.
    + unfold dee2dp_2; simpl fold_right; fold (dee2dp_2 l).
      rewrite dp_size_fix_comp.
      unfold dp_right; rewrite dp_x2y2_size.
      unfold fst, snd.
      generalize (dp_common_size e); intros.
      simpl dp_size at 1; simpl length.
      rewrite Nat.mul_succ_r; omega.
  Qed.

  Let dc_value_1 Ï† Î½ (c : dio_constraint) := 2*(Ï† (fst c)*dee_eval Ï† Î½ (snd c)).
  Let dc_value_2 Ï† Î½ (c : dio_constraint) := (Ï† (fst c)*Ï† (fst c)) + (dee_eval Ï† Î½ (snd c)*dee_eval Ï† Î½ (snd c)).

  Let dee2dp_1_eval Ï† Î½ l : dp_eval Ï† Î½ (dee2dp_1 l) = fold_right plus 0 (map (dc_value_1 Ï† Î½) l).
  Proof.
    induction l as [ | (u,e) l IHl ].
    + simpl; auto.
    + simpl fold_right; rewrite <- IHl.
      unfold dee2dp_1; simpl fold_right.
      rewrite dp_eval_fix_add; f_equal; auto.
      unfold dc_value_1, dp_left, dp_2xy.
      repeat rewrite dp_eval_fix_mul.
      unfold fst, snd; do 2 f_equal.
      apply dp_eval_common.
  Qed.

  Let dee2dp_2_eval Ï† Î½ l : dp_eval Ï† Î½ (dee2dp_2 l) = fold_right plus 0 (map (dc_value_2 Ï† Î½) l).
  Proof.
    induction l as [ | (u,e) l IHl ].
    + simpl; auto.
    + simpl fold_right; rewrite <- IHl.
      unfold dee2dp_2; simpl fold_right.
      rewrite dp_eval_fix_add; f_equal; auto.
      unfold dc_value_2, dp_right, dp_x2y2.
      rewrite dp_eval_fix_add.
      repeat rewrite dp_eval_fix_mul.
      unfold fst, snd; do 2 f_equal;
      apply dp_eval_common.
  Qed.

  Let dee2dp_spec Ï† Î½ l : dp_eval Ï† Î½ (dee2dp_1 l) = dp_eval Ï† Î½ (dee2dp_2 l)
                      <-> Forall (dc_eval Ï† Î½) l.
  Proof.
    rewrite dee2dp_1_eval, dee2dp_2_eval.
    destruct (list_fun_inv l (0,dee_nat 0)) as (f & Hf).
    rewrite Hf at 1 2.
    do 2 rewrite map_map.
    do 2 rewrite <- sum_fold_map.
    unfold dc_value_1, dc_value_2.
    rewrite convex_n_eq.
    unfold dc_eval.
    apply Forall_forall_map with (P := fun i =>  Ï† (fst i) = dee_eval Ï† Î½ (snd i)); auto.
  Qed.

  Theorem dio_elem_single l : { E : dio_single nat nat | dio_single_size E <= 2+20*length l
              /\ forall Î½ Ï†, dp_eval Ï† Î½ (fst E) = dp_eval Ï† Î½ (snd E) <-> Forall (dc_eval Ï† Î½) l }.
  Proof.
    exists (dee2dp_1 l,dee2dp_2 l); split.
    + unfold dio_single_size, fst, snd.
      generalize (dee2dp_1_size l) (dee2dp_2_size l); intros; omega.
    + unfold dio_single_pred, fst, snd; split; apply dee2dp_spec.
  Defined.

  Theorem dio_elem_equation l : { E : dio_single nat nat | dio_single_size E <= 2+20*length l
                                            /\ forall Î½, dio_single_pred E Î½ <-> exists Ï†, Forall (dc_eval Ï† Î½) l }.
  Proof.
    destruct (dio_elem_single l) as (p & H1 & H2); exists p; split; auto.
    split; intros (phi & H); exists phi; revert H; apply H2.
  Defined.

End dio_elem_dio_poly.

Corollary dio_rel_single R : 
      ð�”»R R -> { E : dio_single nat nat | forall Î½, R Î½ <-> dio_single_pred E Î½}.
Proof.
  intros (A & HA).
  destruct dio_formula_elem with (f := A) as (l & _ & _ & Hl).
  destruct dio_elem_equation with (l := l) as (E & _ & HE).
  exists E; intro; rewrite HE, <- Hl, HA; tauto.
Qed.

Section dio_poly_pos.

  Variable P : Set.

  Implicit Type (p : dio_polynomial nat P).

  Definition dio_poly_pos m p : (forall x, In x (dp_var_list p) -> x < m) -> { q | forall Ï† Î½, dp_eval Ï† Î½ p = dp_eval (vec_pos (fun2vec 0 m Ï†)) Î½ q }.
  Proof.
    induction p as [ n | v | i | o p Hp q Hq ]; intros H.
    + exists (dp_nat n); auto.
    + specialize (H v); spec in H; simpl; auto.
      exists (dp_var (nat2pos H)); intros phi psi; simpl.
      rewrite vec_pos_fun2vec, pos2nat_nat2pos; auto.
    + exists (dp_par i); auto.
    + simpl in H.
      destruct Hp as (p1 & H1). { intros; apply H, in_or_app; auto. }
      destruct Hq as (q1 & H2). { intros; apply H, in_or_app; auto. }
      exists (dp_comp o p1 q1); intros phi psi; simpl.
      destruct o; f_equal; auto.
  Qed.

  Theorem dio_poly_eq_pos (e : dio_single nat P) : { m : nat 
                                                 & { p' : dio_polynomial (pos m) P 
                                                 & { q' | forall Î½, dio_single_pred e Î½ <-> dio_single_pred (p',q') Î½ } } }.
  Proof.
    destruct e as (p,q).
    destruct (list_upper_bound (dp_var_list p++dp_var_list q)) as (m & Hm).
    destruct (@dio_poly_pos m p) as (p1 & H1). { intros; apply Hm, in_or_app; auto. }
    destruct (@dio_poly_pos m q) as (q1 & H2). { intros; apply Hm, in_or_app; auto. }
    exists m, p1, q1; intros psi; unfold dio_single_pred.
    split; intros (phi & Hphi).
    + exists (vec_pos (fun2vec 0 m phi)); simpl.
      rewrite <- H1, <- H2; auto.
    + exists (vec2fun (vec_set_pos phi) 0).
      rewrite H1, H2.
      rewrite fun2vec_vec2fun.
      eq goal Hphi; f_equal; 
      apply dp_eval_ext; auto; intros j _; rewrite vec_pos_set; auto. 
  Qed.

  Fact dio_poly_eq_pos_equiv n (p q : dio_polynomial (pos n) P) Î½ : dio_single_pred (p,q) Î½ <-> exists w, dp_eval (vec_pos w) Î½ p = dp_eval (vec_pos w) Î½ q.
  Proof.
    split. 
    + intros (w & Hw); exists (vec_set_pos w); eq goal Hw; f_equal; simpl; apply dp_eval_ext; auto;
        intros; rewrite vec_pos_set; auto.
    + intros (w & Hw); exists (vec_pos w); auto.
  Qed.

End dio_poly_pos.

Check dio_poly_eq_pos.

Section dio_poly_inst_par.

  Variable (V P : Set) (Ïƒ : P -> nat).

  Fixpoint dp_inst_par (p : dio_polynomial V P) : dio_polynomial V Empty_set :=
    match p with
      | dp_nat c       => dp_nat c
      | dp_var v       => dp_var v
      | dp_par p       => dp_nat (Ïƒ p)
      | dp_comp o p q  => dp_comp o (dp_inst_par p) (dp_inst_par q)
    end.

  Fact dp_inst_par_eval Ï† Î½ p : 
    dp_eval Ï† Î½ (dp_inst_par p) = dp_eval Ï† Ïƒ p.
  Proof. induction p as [ | | | [] ]; simpl; f_equal; auto. Qed.

End dio_poly_inst_par.

Section dio_poly_ren_par.

  Variable (V P Q : Set) (f : P -> Q).

  Fixpoint dp_ren_par p : dio_polynomial V Q :=
    match p with
      | dp_nat c       => dp_nat c
      | dp_var v       => dp_var v
      | dp_par p       => dp_par (f p)
      | dp_comp o p q  => dp_comp o (dp_ren_par p) (dp_ren_par q)
    end.

  Fact dp_ren_par_eval Ï† Î½ p : 
    dp_eval Ï† Î½ (dp_ren_par p) = dp_eval Ï† (fun i => Î½ (f i)) p.
  Proof. induction p as [ | | | [] ]; simpl; f_equal; auto. Qed.

End dio_poly_ren_par.

Section dio_poly_proj_par.

  Variable (V : Set) (n : nat).

  Fixpoint dp_proj_par p : dio_polynomial V (pos n) :=
    match p with
      | dp_nat c       => dp_nat c
      | dp_var v       => dp_var v
      | dp_par p       => match le_lt_dec n p with left _ => dp_nat 0 | right H => dp_par (nat2pos H) end
      | dp_comp o p q  => dp_comp o (dp_proj_par p) (dp_proj_par q)
    end.

  Fact dp_proj_par_eval Ï† Î½ p : 
    dp_eval Ï† Î½ (dp_proj_par p) = dp_eval Ï† (fun i => match le_lt_dec n i with left _ => 0 | right H => Î½ (nat2pos H) end) p.
  Proof. 
    induction p as [ | | p | [] ]; simpl; f_equal; auto.
    destruct (le_lt_dec n p); auto.
  Qed.
      
End dio_poly_proj_par.


