From d85683086bed8b76580616306166965079a02fb3 Mon Sep 17 00:00:00 2001
From: David Monniaux <David.Monniaux@univ-grenoble-alpes.fr>
Date: Mon, 13 Dec 2021 18:43:41 +0100
Subject: div_approx_reals_correct

---
 kvx/ExtFloats.v | 56 ++++++++++++++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 54 insertions(+), 2 deletions(-)

(limited to 'kvx/ExtFloats.v')

diff --git a/kvx/ExtFloats.v b/kvx/ExtFloats.v
index b08503a5..f6dae39e 100644
--- a/kvx/ExtFloats.v
+++ b/kvx/ExtFloats.v
@@ -13,7 +13,9 @@
 (*                                                             *)
 (* *************************************************************)
 
-Require Import Floats Integers ZArith.
+From Flocq Require Import Core Digits Operations Round Bracket Sterbenz
+                          Binary Round_odd.
+Require Import Floats Integers ZArith IEEE754_extra Zdiv Psatz.
 
 Module ExtFloat.
 (** TODO check with the actual KVX;
@@ -50,5 +52,55 @@ Definition max (x : float32) (y : float32) : float32 :=
 Definition one := Float32.of_int (Int.repr (1%Z)).
 Definition inv (x : float32) : float32 :=
   Float32.div one x.
-
 End ExtFloat32.
+
+
+Definition div_approx_reals (a b : Z) (x : R) :=
+    let q:=ZnearestE x in
+    let r:=a-q*b in
+    if r <? 0
+    then q-1
+    else q.
+
+Lemma floor_ball1:
+  forall x : R, forall y : Z,
+    (Rabs (x - IZR y) < 1)%R -> Zfloor x = (y-1)%Z \/ Zfloor x = y.
+Proof.
+  intros x y BALL.
+  apply Rabs_lt_inv in BALL.
+  case (Rcompare_spec x (IZR y)); intro CMP.
+  - left. apply Zfloor_imp.
+    replace (y-1+1) with y by ring.
+    rewrite minus_IZR. lra.
+  - subst.
+    rewrite Zfloor_IZR. right. reflexivity.
+  - right. apply Zfloor_imp.
+    rewrite plus_IZR. lra.
+Qed.
+
+Theorem div_approx_reals_correct:
+  forall a b : Z, forall x : R,
+    b > 0 ->
+    (Rabs (x - IZR a/ IZR b) < 1/2)%R ->
+    div_approx_reals a b x = (a/b)%Z.
+Proof.                    
+  intros a b x bPOS GAP.
+  assert (0 < IZR b)%R by (apply IZR_lt ; lia).
+  unfold div_approx_reals.
+  pose proof (Znearest_imp2 (fun x => negb (Z.even x)) x) as NEAR.
+  assert (Rabs (IZR (ZnearestE x) - IZR a/ IZR b) < 1)%R as BALL.
+  { pose proof (Rabs_triang (IZR (ZnearestE x) - x)
+                            (x - IZR a/ IZR b)) as TRI.
+    ring_simplify (IZR (ZnearestE x) - x + (x - IZR a / IZR b))%R in TRI.
+    lra.
+  }
+  clear GAP NEAR.
+  rewrite Rabs_minus_sym in BALL.
+  pose proof (floor_ball1 _ _ BALL) as FLOOR.
+  clear BALL.
+  rewrite Zfloor_div in FLOOR by lia.
+  pose proof (Z_div_mod_eq_full a b) as DIV_MOD.
+  assert (0 < b) as bPOS' by lia.
+  pose proof (Z.mod_pos_bound a b bPOS') as MOD_BOUNDS.
+  case Z.ltb_spec; intro; destruct FLOOR; lia.
+Qed.
-- 
cgit