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.

Fixpoint Dnormalize : diadic := Cases (Dnum x) of (POS (xO y)) => (Dnormalize (Diadic (POS y) `(Dexp x)+1`)) | (NEG (xO y)) => (Dnormalize (Diadic (NEG y) `(Dexp x)+1`)) | _ => x end.
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_supRounding_inf
| Rounding_infRounding_sup
| Rounding_nearestRounding_nearest
| Rounding_zeroRounding_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.EqTrue
| Datatypes.LtFalse
| Datatypes.GtFalse
end.
Definition Dneq (x y : diadic) :=
match Dcompare x y with
| Datatypes.EqFalse
| Datatypes.LtTrue
| Datatypes.GtTrue
end.
Definition Dle (x y : diadic) :=
match Dcompare x y with
| Datatypes.EqTrue
| Datatypes.LtTrue
| Datatypes.GtFalse
end.
Definition Dlt (x y : diadic) :=
match Dcompare x y with
| Datatypes.EqFalse
| Datatypes.LtTrue
| Datatypes.GtFalse
end.
Definition Dge (x y : diadic) :=
match Dcompare x y with
| Datatypes.EqTrue
| Datatypes.LtFalse
| Datatypes.GtTrue
end.
Definition Dgt (x y : diadic) :=
match Dcompare x y with
| Datatypes.EqFalse
| Datatypes.LtFalse
| Datatypes.GtTrue
end.

Definition Deq_bool (x y : diadic) :=
match Dcompare x y with
| Datatypes.Eqtrue
| Datatypes.Ltfalse
| Datatypes.Gtfalse
end.
Definition Dneq_bool (x y : diadic) :=
match Dcompare x y with
| Datatypes.Eqfalse
| Datatypes.Lttrue
| Datatypes.Gttrue
end.
Definition Dle_bool (x y : diadic) :=
match Dcompare x y with
| Datatypes.Eqtrue
| Datatypes.Lttrue
| Datatypes.Gtfalse
end.
Definition Dlt_bool (x y : diadic) :=
match Dcompare x y with
| Datatypes.Eqfalse
| Datatypes.Lttrue
| Datatypes.Gtfalse
end.
Definition Dge_bool (x y : diadic) :=
match Dcompare x y with
| Datatypes.Eqtrue
| Datatypes.Ltfalse
| Datatypes.Gttrue
end.
Definition Dgt_bool (x y : diadic) :=
match Dcompare x y with
| Datatypes.Eqfalse
| Datatypes.Ltfalse
| Datatypes.Gttrue
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
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
| Z0y = 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_infDle y x Dlt x (Dsucc y)
| Rounding_supDlt (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 qZROUND m q nx
| Zneg q ⇒ (nx × two_power_pos q)%Z
| Z0nx
end.