more on FPDivision

1 files changed, 133 insertions, 3 deletions

diff --git a/kvx/FPDivision.v b/kvx/FPDivision.v
index 4e7fce31..ebf9136d 100644
--- a/kvx/FPDivision.v
+++ b/kvx/FPDivision.v
@@ -1,5 +1,5 @@
 From Flocq Require Import Core Digits Operations Round Bracket Sterbenz
-                          Binary Round_odd.
+                          Binary Round_odd Bits.
 Require Archi.
 Require Import Coqlib.
 Require Import Compopts.
@@ -9,12 +9,14 @@ Require Import Floats.
 Require Import Op.
 Require Import CminorSel.
 Require Import OpHelpers.
-Require Import ExtValues ExtFloats.
+Require Import ExtFloats.
 Require Import DecBoolOps.
 Require Import Chunks.
 Require Import Builtins.
 Require Import Values Globalenvs.
 Require Compopts.
+Require Import Psatz.
+Require Import IEEE754_extra.
 
 Local Open Scope cminorsel_scope.
 
@@ -31,13 +33,15 @@ Definition approx_inv b :=
   let x := Eop Ofmaddf (invb_d_var ::: alpha ::: invb_d_var ::: Enil) in
   Elet b (Elet invb_d x).
 
+Definition approx_inv_thresh := (0.000000000000001)%R.
+
 Theorem approx_inv_correct :
   forall (ge : genv)  (sp: val) cmenv memenv
          (le : letenv) (expr_b : expr) b
          (EVAL_b : eval_expr ge sp cmenv memenv le expr_b (Vint b)),
     exists f : float,
       eval_expr ge sp cmenv memenv le (approx_inv expr_b) (Vfloat f) /\
-        is_finite _ _ f = true /\ (Rabs((B2R _ _ f) - (1 / IZR (Int.unsigned b))) <= 0.1)%R.
+        is_finite _ _ f = true /\ (Rabs((B2R _ _ f) - (1 / IZR (Int.unsigned b))) <= approx_inv_thresh)%R.
 Proof.
   intros. unfold approx_inv.
   repeat econstructor.
@@ -47,3 +51,129 @@ Proof.
   all: set (b_d := (Float.of_longu (Int64.repr (Int.unsigned b)))).
   all: cbn.
 Admitted.
+
+Definition approx_div a b :=
+  let a_var := Eletvar (1%nat) in
+  let b_var := Eletvar (0%nat) in
+  let a_d := Eop Ofloatoflongu ((Eop Ocast32unsigned (a_var ::: Enil)) ::: Enil) in
+  let qr := Eop Olonguoffloat_ne ((Eop Omulf (a_d:::(approx_inv b_var):::Enil)):::Enil) in
+  let qr_var := Eletvar 0%nat in
+  let rem := Eop Omsubl ((Eop Ocast32unsigned ((Eletvar (2%nat)):::Enil)):::
+                         qr_var :::
+                         (Eop Ocast32unsigned ((Eletvar (1%nat)):::Enil)):::Enil) in
+  let qr_m1 := Eop (Oaddlimm (Int64.repr (-1)%Z)) (qr_var:::Enil) in
+  let cases := Eop (Osel (Ccompl0 Clt) Tlong)
+                   (qr_m1 ::: qr_var ::: rem ::: Enil) in (* (Elet qr cases) *)
+  Eop Olowlong ((Elet a (Elet (lift b) (Elet qr cases))) ::: Enil).
+
+
+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.
+
+Opaque approx_inv.
+
+Theorem approx_div_correct :
+  forall (ge : genv)  (sp: val) cmenv memenv
+         (le : letenv) (expr_a expr_b : expr) (a b : int)
+         (EVAL_a : eval_expr ge sp cmenv memenv le expr_a (Vint a))
+         (EVAL_b : eval_expr ge sp cmenv memenv le expr_b (Vint b))
+         (b_nz : (Int.unsigned b > 0)%Z),
+  eval_expr ge sp cmenv memenv le (approx_div expr_a expr_b)
+            (Vint (Int.divu a b)).
+Proof.
+  intros.
+  assert (eval_expr ge sp cmenv memenv (Vint b :: Vint a :: le)
+                    (Eletvar 0) (Vint b)) as EVAL_b'.
+  { constructor. reflexivity. }
+  destruct (approx_inv_correct ge sp cmenv memenv (Vint b :: Vint a :: le) (Eletvar 0) b EVAL_b') as (inv_b & inv_b_eval & inv_b_finite & inv_b_correct).
+  unfold approx_div.
+  repeat econstructor.
+  { exact EVAL_a. }
+  { apply eval_lift. exact EVAL_b. }
+  exact inv_b_eval.
+  cbn. f_equal.
+  rewrite <- Float.of_intu_of_longu.
+  unfold Float.to_longu_ne.
+  rewrite ZofB_ne_range_correct by lia.
+  set (prod := (Float.mul (Float.of_intu a) inv_b)).
+  set (a' := Int.unsigned a) in *.
+  set (b' := Int.unsigned b) in *.
+  set (prod' := (B2R 53 1024 prod)).
+  set (prod_r := ZnearestE prod') in *.
+  replace (_ && _ && _) with true by admit.
+  cbn.
+  f_equal.
+  unfold Int.divu.
+  rewrite <- div_approx_reals_correct with (x := B2R _ _ prod).
+  2: lia.
+  {
+    unfold div_approx_reals.
+    fold a' b' prod_r.
+    unfold Int64.mul.
+    rewrite Int64.unsigned_repr by admit.
+    rewrite Int64.unsigned_repr by admit.
+    unfold Int64.sub.
+    rewrite Int64.unsigned_repr by admit.
+    rewrite Int64.unsigned_repr by admit.
+    case Z.ltb_spec; intro CMP.
+    - replace (Int64.lt (Int64.repr (a' - prod_r * b')) Int64.zero) with true by admit.
+      cbn.
+      f_equal.
+      admit.
+    - replace (Int64.lt (Int64.repr (a' - prod_r * b')) Int64.zero) with false by admit.
+      cbn.
+      f_equal.
+      admit.
+  }
+  fold prod'.
+  admit.
+Admitted.
+    
+