Require Import Omega.
Require Import Zcomplements.
Require Import Zpower.
Require Import Zlogarithm.
Require Import ZArithRing.

Section definitions.


Record diadic : Set := Diadic {Dnum : Z; Dexp : Z}.

Definition Dshift (n : Z) (x : diadic) := Diadic (Dnum x) (Dexp x + n).

Definition Dzero (x : Z) := Diadic 0 x.

Definition is_Dzero (x : diadic) := Dnum x = 0%Z.

Inductive rounding_mode : Set :=
  | Rounding_sup : rounding_mode
  | Rounding_inf : rounding_mode
  | Rounding_nearest : rounding_mode
  | Rounding_zero : rounding_mode.

Definition Rounding_mode_opp (m : rounding_mode) :=
  match m with
  | Rounding_sup ⇒ Rounding_inf
  | Rounding_inf ⇒ Rounding_sup
  | Rounding_nearest ⇒ Rounding_nearest
  | Rounding_zero ⇒ Rounding_zero
  end.



End definitions.

Section comparisons.


Definition Dcompare (x y : diadic) : Datatypes.comparison :=
  let nx := Dnum x in
  let ny := Dnum y in
  let ex := Dexp x in
  let ey := Dexp y in
  (two_p (ex - Zmin ex ey) × nx ?= two_p (ey - Zmin ex ey) × ny)%Z.

Definition Deq (x y : diadic) :=
  match Dcompare x y with
  | Datatypes.Eq ⇒ True
  | Datatypes.Lt ⇒ False
  | Datatypes.Gt ⇒ False
  end.
Definition Dneq (x y : diadic) :=
  match Dcompare x y with
  | Datatypes.Eq ⇒ False
  | Datatypes.Lt ⇒ True
  | Datatypes.Gt ⇒ True
  end.
Definition Dle (x y : diadic) :=
  match Dcompare x y with
  | Datatypes.Eq ⇒ True
  | Datatypes.Lt ⇒ True
  | Datatypes.Gt ⇒ False
  end.
Definition Dlt (x y : diadic) :=
  match Dcompare x y with
  | Datatypes.Eq ⇒ False
  | Datatypes.Lt ⇒ True
  | Datatypes.Gt ⇒ False
  end.
Definition Dge (x y : diadic) :=
  match Dcompare x y with
  | Datatypes.Eq ⇒ True
  | Datatypes.Lt ⇒ False
  | Datatypes.Gt ⇒ True
  end.
Definition Dgt (x y : diadic) :=
  match Dcompare x y with
  | Datatypes.Eq ⇒ False
  | Datatypes.Lt ⇒ False
  | Datatypes.Gt ⇒ True
  end.

Definition Deq_bool (x y : diadic) :=
  match Dcompare x y with
  | Datatypes.Eq ⇒ true
  | Datatypes.Lt ⇒ false
  | Datatypes.Gt ⇒ false
  end.
Definition Dneq_bool (x y : diadic) :=
  match Dcompare x y with
  | Datatypes.Eq ⇒ false
  | Datatypes.Lt ⇒ true
  | Datatypes.Gt ⇒ true
  end.
Definition Dle_bool (x y : diadic) :=
  match Dcompare x y with
  | Datatypes.Eq ⇒ true
  | Datatypes.Lt ⇒ true
  | Datatypes.Gt ⇒ false
  end.
Definition Dlt_bool (x y : diadic) :=
  match Dcompare x y with
  | Datatypes.Eq ⇒ false
  | Datatypes.Lt ⇒ true
  | Datatypes.Gt ⇒ false
  end.
Definition Dge_bool (x y : diadic) :=
  match Dcompare x y with
  | Datatypes.Eq ⇒ true
  | Datatypes.Lt ⇒ false
  | Datatypes.Gt ⇒ true
  end.
Definition Dgt_bool (x y : diadic) :=
  match Dcompare x y with
  | Datatypes.Eq ⇒ false
  | Datatypes.Lt ⇒ false
  | Datatypes.Gt ⇒ true
  end.

Lemma Dcompare_shift :
 ∀ (x y : diadic) (n : Z),
 Dcompare (Dshift n x) (Dshift n y) = Dcompare x y.
unfold Dcompare in |- *; simpl in |- *; intros;
 rewrite (Zmin.Zmin_plus (Dexp x) (Dexp y) n).
do 2 rewrite BinInt.Zminus_plus_simpl_r.
reflexivity.
Qed.

Lemma eq_Deq : ∀ x y : diadic, x = y → Deq x y.
intros; rewrite H; unfold Deq in |- *; unfold Dcompare in |- *;
 apply Zcompare.Zcompare_eq_case; trivial.
Qed.


Lemma Dcompare_zero :
 ∀ (x : diadic) (n : Z), Dcompare x (Dzero n) = (Dnum x ?= 0)%Z.
intros (nx, ex) n.
unfold Dcompare in |- *; simpl in |- ×.
symmetry in |- ×.
replace (two_p (n - Zmin ex n) × 0)%Z with (two_p (ex - Zmin ex n) × 0)%Z.
apply Zmult_compare_compat_l.
apply two_p_gt_ZERO.
generalize (Zle_min_l ex n); generalize (Zmin ex n); intro; omega.
do 2 rewrite Zmult_0_r; reflexivity.
Qed.

Lemma Deq_shift :
 ∀ (x y : diadic) (n : Z), Deq x y → Deq (Dshift n x) (Dshift n y).
unfold Deq in |- *; intros; rewrite (Dcompare_shift x y n); trivial.
Qed.

Lemma Deq_x_shift_x :
 ∀ (x : diadic) (n : Z),
 (0 ≤ n)%Z → Deq x (Diadic (Dnum x × two_p n) (Dexp x - n)).
intros (nx, ex) n Hn; unfold Deq in |- *; unfold Dcompare in |- *;
 simpl in |- ×.
cut
 ((two_p (ex - Zmin ex (ex - n)) × nx)%Z =
  (two_p (ex - n - Zmin ex (ex - n)) × (nx × two_p n))%Z).
intro H; rewrite H.
generalize (two_p (ex - n - Zmin ex (ex - n)) × (nx × two_p n))%Z.
intros. generalize (Zcompare.Zcompare_refl z).
elim (z ?= z)%Z; discriminate || trivial.
rewrite (Zmult_comm nx (two_p n)).
rewrite <- Zmult_assoc_reverse.
rewrite <- two_p_is_exp.
ring_simplify (ex - n - Zmin ex (ex - n) + n)%Z (ex - Zmin ex (ex - n))%Z.
reflexivity.
generalize (Zle_min_l ex (ex - n)) (Zle_min_r ex (ex - n)).
omega.
assumption.
Qed.

Lemma Dle_Zle :
 ∀ n1 n2 d : Z, (n1 ≤ n2)%Z → Dle (Diadic n1 d) (Diadic n2 d).
intros; unfold Dle in |- *; unfold Dcompare in |- *; simpl in |- ×.
rewrite (Zmin_n_n d).
rewrite <- (Zminus_diag_reverse d).
unfold two_p in |- ×.
rewrite (Zcompare_mult_compat 1 n1 n2).
apply Zcompare.Zle_compare; assumption.
Qed.

Lemma Dlt_Zlt :
 ∀ n1 n2 d : Z, (n1 < n2)%Z → Dlt (Diadic n1 d) (Diadic n2 d).
intros; unfold Dlt in |- *; unfold Dcompare in |- *; simpl in |- ×.
rewrite (Zmin_n_n d).
rewrite <- (Zminus_diag_reverse d).
unfold two_p in |- ×.
rewrite (Zcompare_mult_compat 1 n1 n2).
apply Zcompare.Zlt_compare; assumption.
Qed.

Lemma Dge_Zge :
 ∀ n1 n2 d : Z, (n1 ≥ n2)%Z → Dge (Diadic n1 d) (Diadic n2 d).
intros; unfold Dge in |- *; unfold Dcompare in |- *; simpl in |- ×.
rewrite (Zmin_n_n d).
rewrite <- (Zminus_diag_reverse d).
unfold two_p in |- ×.
rewrite (Zcompare_mult_compat 1 n1 n2).
apply Zcompare.Zge_compare; assumption.
Qed.

Lemma Dgt_Zgt :
 ∀ n1 n2 d : Z, (n1 > n2)%Z → Dgt (Diadic n1 d) (Diadic n2 d).
intros; unfold Dgt in |- *; unfold Dcompare in |- *; simpl in |- ×.
rewrite (Zmin_n_n d).
rewrite <- (Zminus_diag_reverse d).
unfold two_p in |- ×.
rewrite (Zcompare_mult_compat 1 n1 n2).
apply Zcompare.Zgt_compare; assumption.
Qed.

Lemma Dle_refl : ∀ x y : diadic, Deq x y → Dle x y.
unfold Deq in |- *; unfold Dle in |- *; intros x y; elim (Dcompare x y);
 trivial.
Qed.


End comparisons.

Section operations.


Definition Dsucc (x : diadic) := Diadic (Dnum x + 1) (Dexp x).
Definition Dpred (x : diadic) := Diadic (Dnum x - 1) (Dexp x).

Definition Dadd (x y : diadic) :=
  let nx := Dnum x in
  let ny := Dnum y in
  let ex := Dexp x in
  let ey := Dexp y in
  Diadic (two_p (ex - Zmin ex ey) × nx + two_p (ey - Zmin ex ey) × ny)
    (Zmin ex ey).

Definition Dopp (x : diadic) := Diadic (- Dnum x) (Dexp x).

Definition Dabs (x : diadic) := Diadic (Zabs (Dnum x)) (Dexp x).

Definition Dminus (x y : diadic) := Dadd x (Dopp y).

Definition Dmult (x y : diadic) := Diadic (Dnum x × Dnum y) (Dexp x + Dexp y).



Definition Dproj (m : rounding_mode) (x : diadic) (P : diadic → Prop)
  (y : diadic) :=
  P y ∧
  match m with
  | Rounding_sup ⇒ ∀ z : diadic, P z → Dle x z → Dle y z
  | Rounding_inf ⇒ ∀ z : diadic, P z → Dle z x → Dle z y
  | Rounding_nearest ⇒
      ∀ z : diadic, P z → Dle (Dabs (Dminus x y)) (Dabs (Dminus x z))
  | Rounding_zero ⇒
      ∀ z : diadic,
      P z →
      IF Dle (Dzero 0) x then Dle z x → Dle z y else Dle z x → Dle z y
  end.


Lemma ZROUND_inf_spec :
 ∀ (p : positive) (x : Z),
 {y : Z | (y × two_power_pos p ≤ x < Zsucc y × two_power_pos p)%Z}.
intros; elim (Zdiv_rest_correct x p); intros q r Hx Hr1 Hr2; ∃ q;
 rewrite (Zplus_0_r_reverse (q × two_power_pos p));
 rewrite Hx; split;
 [ apply Zplus_le_compat_l; assumption
 | unfold Zsucc in |- *; rewrite Zmult_plus_distr_l; apply Zplus_lt_compat_l;
    rewrite Zmult_1_l; assumption ].
Qed.

Definition ZROUND_inf (p : positive) (x : Z) :=
  let (x', p) := ZROUND_inf_spec p x in x'.

Lemma ZROUND_sup_spec :
 ∀ (p : positive) (x : Z),
 {y : Z | (Zpred y × two_power_pos p < x ≤ y × two_power_pos p)%Z}.
intros; elim (Zdiv_rest_correct x p); intros q r; elim r;
 [ intros Hx Hr; ∃ q; rewrite Hx; rewrite <- Zplus_0_r_reverse; split;
    [ apply Zmult_gt_0_lt_compat_r;
       [ compute in |- *; reflexivity | apply Zlt_pred ]
    | apply Zle_refl ]
 | intros p0 Hx Hr1 Hr2; ∃ (Zsucc q); rewrite Hx; split;
    [ replace (Zpred (Zsucc q) × two_power_pos p)%Z with
       (q × two_power_pos p + 0)%Z;
       [ apply Zplus_lt_compat_l; compute in |- *; reflexivity
       | rewrite <- Zpred_succ; rewrite <- Zplus_0_r_reverse; reflexivity ]
    | unfold Zsucc in |- *; rewrite Zmult_plus_distr_l;
       apply Zplus_le_compat_l; rewrite Zmult_1_l;
       apply Zlt_le_weak; assumption ]
 | intros p0 Hx Hr1 Hr2; absurd (Datatypes.Gt = Datatypes.Gt);
    [ exact Hr1 | reflexivity ] ].
Qed.

Definition ZROUND_sup (p : positive) (x : Z) :=
  let (x', p) := ZROUND_sup_spec p x in x'.

Lemma ZROUND_correct :
 ∀ (m : rounding_mode) (p : positive) (x : Z),
 {y : Z |
 match m with
 | Rounding_inf ⇒ (y × two_power_pos p ≤ x < Zsucc y × two_power_pos p)%Z
 | Rounding_sup ⇒ (Zpred y × two_power_pos p < x ≤ y × two_power_pos p)%Z
 | Rounding_nearest ⇒
     match (x - ZROUND_inf p x ?= ZROUND_sup p x - x)%Z with
     | Datatypes.Eq ⇒
         if Zeven.Zeven_bool (ZROUND_inf p x)
         then (y × two_power_pos p ≤ x < Zsucc y × two_power_pos p)%Z
         else (Zpred y × two_power_pos p < x ≤ y × two_power_pos p)%Z
     | Datatypes.Gt ⇒
         (Zpred y × two_power_pos p < x ≤ y × two_power_pos p)%Z
     | Datatypes.Lt ⇒
         (y × two_power_pos p ≤ x < Zsucc y × two_power_pos p)%Z
     end
 | Rounding_zero ⇒
     match x with
     | Zpos _ ⇒ (y × two_power_pos p ≤ x < Zsucc y × two_power_pos p)%Z
     | Z0 ⇒ y = 0%Z
     | Zneg _ ⇒ (Zpred y × two_power_pos p < x ≤ y × two_power_pos p)%Z
     end
 end}.
simple destruct m;
 [ exact ZROUND_sup_spec
 | exact ZROUND_inf_spec
 | intros p x; elim (x - ZROUND_inf p x ?= ZROUND_sup p x - x)%Z;
    [ elim (Zeven.Zeven_bool (ZROUND_inf p x));
       [ apply ZROUND_inf_spec | apply ZROUND_sup_spec ]
    | apply ZROUND_inf_spec
    | apply ZROUND_sup_spec ]
 | simple induction x;
    [ ∃ 0%Z; reflexivity
    | intro; apply ZROUND_inf_spec
    | intro; apply ZROUND_sup_spec ] ].
Qed.

Definition ZROUND (m : rounding_mode) (p : positive)
  (x : Z) := let (x', p) := ZROUND_correct m p x in x'.

Definition POS_ROUND (m : rounding_mode) (p n : positive) :=
  BinInt.Zabs_N (ZROUND m p (Zpos n)).

Definition NEG_ROUND (m : rounding_mode) (p n : positive) :=
  BinInt.Zabs_N (- ZROUND m p (Zneg n)).



Definition Ddouble (d : diadic) := Dshift 1 d.

Axiom
  ROUND_spec :
    ∀ (m : rounding_mode) (p : Z) (x : diadic),
    {y : diadic |
    N_digits (Dexp y) = p ∧
    match m with
    | Rounding_inf ⇒ Dle y x ∧ Dlt x (Dsucc y)
    | Rounding_sup ⇒ Dlt (Dpred y) x ∧ Dle x y
    | Rounding_nearest ⇒
        Dle (Dpred (Ddouble y)) (Ddouble x) ∧
        Dle (Ddouble x) (Dsucc (Ddouble y))
    | Rounding_zero ⇒
        IF Dlt (Dzero 0) x then Dle y x ∧ Dlt x (Dsucc y)
        else Dlt (Dpred y) x ∧ Dle x y
    end}.

Definition ROUND (m : rounding_mode) (p : Z) (d : diadic) :=
  let (x, _) := ROUND_spec m p d in x.


Definition ANTIROUND (m : rounding_mode) (p : Z) (x : diadic) :=
  let nx := Dnum x in
  let ex := Dexp x in
  match (p - ex)%Z with
  | Zpos q ⇒ ZROUND m q nx
  | Zneg q ⇒ (nx × two_power_pos q)%Z
  | Z0 ⇒ nx
  end.



Parameter Ddiv : rounding_mode → Z → diadic → diadic → diadic.


Parameter Dsqrt : rounding_mode → Z → positive → Z → diadic.

End operations.