From 185846ddf99f5c9f17b74e40fc9ad1deb91cb1ff Mon Sep 17 00:00:00 2001
From: David Monniaux <David.Monniaux@univ-grenoble-alpes.fr>
Date: Wed, 12 Jan 2022 11:35:05 +0100
Subject: finer proof

---
 kvx/FPDivision64.v | 105 +++++++++++++++++++++++++++++++++++++++++++++++------
 1 file changed, 94 insertions(+), 11 deletions(-)

(limited to 'kvx')

diff --git a/kvx/FPDivision64.v b/kvx/FPDivision64.v
index f20efa4e..425e11da 100644
--- a/kvx/FPDivision64.v
+++ b/kvx/FPDivision64.v
@@ -289,14 +289,17 @@ Definition rough_approx_inv_longu b :=
 
 Definition rough_approx_inv_thresh := (3/33554432)%R.
 (* 8.94069671630859e-8 *)
-    
+
+Local Notation "'rd'" := (round radix2 (FLT_exp (-1074) 53) ZnearestE).
+Local Notation "'rs'" := (round radix2 (FLT_exp (-149) 24) ZnearestE).
+
 Theorem rough_approx_inv_longu_correct :
   forall b,
     (0 < Int64.unsigned b)%Z ->
     exists (f : float),
       (rough_approx_inv_longu (Vlong b)) = Vfloat f /\
+        (B2R _ _ f) = (rd (rs (1 / rs (IZR (Int64.unsigned b))))) /\
         is_finite _ _ f = true /\
-        (Rabs((B2R _ _ f) - (1 / IZR (Int64.unsigned b))) <= rough_approx_inv_thresh)%R /\
         Bsign _ _ f = false.
 Proof.
   intros b NONZ.
@@ -350,6 +353,7 @@ Proof.
     cbn.
     gappa.
   }
+  rewrite C1R in C2.
   destruct C2 as (C2R & C2F & C2Sz).
   rewrite C1S in C2Sz.
   change (xorb _ _) with false in C2Sz.
@@ -382,7 +386,6 @@ Proof.
     cbn.
     gappa.
   }
-  
   destruct C3 as (C3R & C3F & C3S).
   set (invb_d :=  (Bconv 24 128 53 1024 re re Float.of_single_nan mode_NE invb_s)) in *.
   assert ((1/18446744073709551616 <= B2R 53 1024 invb_d <= 1)%R) as invb_d_RANGE.
@@ -393,15 +396,10 @@ Proof.
     gappa.
   }
   rewrite C2S in C3S.
+  rewrite C2R in C3R.
+  rewrite C0R in C3R.
 
-  split. assumption. split. 2: assumption.
-  unfold rough_approx_inv_thresh.
-  rewrite C3R.
-  rewrite C2R.
-  rewrite C1R.
-  rewrite C0R.
-  cbn.
-  gappa.
+  auto.
 Qed.
 
 Theorem rough_approx_inv_longu_correct1 :
@@ -611,3 +609,88 @@ Proof.
   set (delta3 := (f_r - 1 / IZR b')%R) in *.
   gappa.
 Qed.
+
+Definition rough_approx_euclid_div_longu a b :=
+  Val.maketotal (Val.longuoffloat_ne (rough_approx_div_longu a b)).
+
+Definition smallb_threshold := 4398046511104%Z.
+Definition rough_approx_div_longu_smallb_delta := 3917010173952%Z.
+
+Theorem rough_approx_div_longu_smallb :
+  forall a b,
+    (2 <= Int64.unsigned b <= smallb_threshold)%Z ->
+    exists (qr : int64),
+      (rough_approx_euclid_div_longu (Vlong a) (Vlong b)) = Vlong qr /\
+        (Z.abs (Int64.unsigned a - Int64.unsigned b* Int64.unsigned qr) <= rough_approx_div_longu_smallb_delta)%Z /\
+        (0 <= Int64.unsigned qr <= 11529215046068469760)%Z.
+Proof.
+  intros a b b_RANGE.
+  unfold rough_approx_euclid_div_longu.
+  destruct (rough_approx_div_longu_correct a b) as (q & C1R & C1D & C1F & C1S).
+  lia.
+  rewrite C1R.
+  cbn.
+  set (qr := B2R 53 1024 q) in *.
+  pose proof (Int64.unsigned_range a) as a_RANGE.
+  change Int64.modulus with 18446744073709551616%Z in a_RANGE.
+  set (a' := Int64.unsigned a) in *.
+  assert (0 <= IZR a' <= 18446744073709551615)%R as IZR_a_RANGE.
+  { split ; apply IZR_le; lia.
+  }
+
+  set (b' := Int64.unsigned b) in *.
+  assert (2 <= IZR b' <= IZR smallb_threshold)%R as IZR_b_RANGE.
+  { split ; apply IZR_le; lia.
+  }
+
+  assert (0 <= qr <= 11529215046068469760)%R as q_RANGE.
+  { split.
+    { unfold qr.
+      apply Bsign_false_nonneg.
+      assumption.
+    }
+    unfold rough_approx_div_thresh in C1D.
+    replace qr with ((qr - IZR a' / IZR b') + (IZR a' / IZR b'))%R by (field ; lra).
+    set (delta := (qr - IZR a' / IZR b')%R) in *.
+    unfold smallb_threshold in *.
+    gappa.
+  }
+  unfold Float.to_longu_ne.
+  pose proof (ZofB_ne_range_correct _ _ q 0 Int64.max_unsigned) as C2.
+  rewrite C1F in C2.
+  cbn in C2.
+  fold qr in C2.
+  assert (0 <= ZnearestE qr <= 11529215046068469760)%Z as round_RANGE.
+  { split.
+    { apply Znearest_IZR_le.
+      lra.
+    }
+    apply Znearest_le_IZR.
+    lra.
+  }
+  rewrite Zle_bool_true in C2 by lia.
+  rewrite Zle_bool_true in C2 by lia.
+  change Int64.max_unsigned with 18446744073709551615%Z.
+  cbn in C2.
+  rewrite C2.
+  cbn.
+  econstructor. split. reflexivity.
+  rewrite Int64.unsigned_repr; cycle 1.
+  { change Int64.max_unsigned with 18446744073709551615%Z.
+    lia.
+  }
+  split.
+  2: assumption.
+  unfold rough_approx_div_thresh, rough_approx_div_longu_smallb_delta in *.
+  apply le_IZR.
+  rewrite <- Rabs_Zabs.
+  rewrite minus_IZR.
+  rewrite mult_IZR.
+  replace  (IZR a' - IZR b' * IZR (ZnearestE qr))%R with
+    (- IZR b' * (IZR (ZnearestE qr) - qr)
+     - IZR b' * (qr - IZR a' / IZR b'))%R by (field ; lra).
+  pose proof (Znearest_imp2  (fun x : Z => negb (Z.even x)) qr) as ROUND.
+  set (delta1 := (IZR (ZnearestE qr) - qr)%R) in *.
+  set (delta2 := (qr - IZR a' / IZR b')%R) in *.
+  unfold smallb_threshold in *.
+ 
-- 
cgit