1 files changed, 61 insertions, 0 deletions
diff --git a/kvx/ExtValues.v b/kvx/ExtValues.v
index c493f708..f08c1157 100644
--- a/kvx/ExtValues.v
+++ b/kvx/ExtValues.v
@@ -755,6 +755,67 @@ Definition fmaddfs := triple_op_single (fun f1 f2 f3 => Float32.fma f2 f3 f1).
 Definition fmsubf := triple_op_float (fun f1 f2 f3 => Float.fma (Float.neg f2) f3 f1).
 Definition fmsubfs := triple_op_single (fun f1 f2 f3 => Float32.fma (Float32.neg f2) f3 f1).
 
+From Flocq Require Import Core Digits Operations Round Bracket Sterbenz
+                          Binary Round_odd.
+Require Import  IEEE754_extra Zdiv Psatz Floats ExtFloats.
+
+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.
+    ring_simplify (y-1+1).
+    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.
+
+Definition my_div (a b : val) :=
+  let b_d := Val.maketotal (Val.floatofintu b) in
+  let invb_d := Val.floatofsingle (invfs (Val.maketotal (Val.singleofintu b))) in
+  let alpha := fmsubf (Vfloat ExtFloat.one) invb_d b_d in
+  let x := fmaddf invb_d alpha invb_d in
+  Val.mulf (Val.maketotal (Val.floatofintu a)) x.
+
 Definition int_abs i1 := Int.repr (Z.abs (Int.signed i1)).
 Definition long_abs i1 := Int64.repr (Z.abs (Int64.signed i1)).