From 120446974d1a83e6f9ef99f049d597166be09271 Mon Sep 17 00:00:00 2001
From: David Monniaux <David.Monniaux@univ-grenoble-alpes.fr>
Date: Mon, 10 Jan 2022 14:18:07 +0100
Subject: rough approx

---
 kvx/FPDivision64.v | 116 +++++++++++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 113 insertions(+), 3 deletions(-)

(limited to 'kvx')

diff --git a/kvx/FPDivision64.v b/kvx/FPDivision64.v
index ab593591..6c053d7c 100644
--- a/kvx/FPDivision64.v
+++ b/kvx/FPDivision64.v
@@ -37,9 +37,6 @@ Definition approx_inv_longu b :=
   let alpha := ExtValues.fmsubf one invb_d b_d in
   ExtValues.fmaddf invb_d alpha invb_d.
 
-Definition approx_inv_thresh := (25/2251799813685248)%R.
-(* 1.11022302462516e-14 *)
-
 Lemma Rabs_relax:
   forall b b' (INEQ : (b < b')%R) x,
     (-b <= x <= b)%R -> (Rabs x < b')%R.
@@ -48,6 +45,9 @@ Proof.
   apply Rabs_lt.
   lra.
 Qed.
+
+Definition approx_inv_thresh := (25/2251799813685248)%R.
+(* 1.11022302462516e-14 *)
     
 Theorem approx_inv_longu_correct :
   forall b,
@@ -282,3 +282,113 @@ Proof.
   unfold approx_inv_thresh.
   gappa.
 Qed.
+
+Definition rough_approx_inv_longu b :=
+  let invb_s := ExtValues.invfs (Val.maketotal (Val.singleoflongu b)) in
+  Val.floatofsingle invb_s.
+
+Definition rough_approx_inv_thresh := (3/33554432)%R.
+(* 8.94069671630859e-8 *)
+    
+Theorem rough_approx_inv_longu_correct :
+  forall b,
+    (0 < Int64.unsigned b)%Z ->
+    exists (f : float),
+      (rough_approx_inv_longu (Vlong b)) = Vfloat f /\
+      is_finite _ _ f = true /\ (Rabs((B2R _ _ f) - (1 / IZR (Int64.unsigned b))) <= rough_approx_inv_thresh)%R.
+Proof.
+  intros b NONZ.
+  unfold rough_approx_inv_longu.
+  cbn.
+  econstructor.
+  split.
+  reflexivity.
+  Local Transparent Float.neg Float.of_single Float32.of_longu Float32.div Float.of_longu Float32.of_int Float.of_int.
+  unfold Float.fma, Float.neg, Float.of_single, Float32.of_longu, ExtFloat32.inv, Float32.div, Float.of_longu, ExtFloat32.one, Float32.of_int, ExtFloat.one, Float.of_int.
+  set (re := (@eq_refl Datatypes.comparison Lt)).
+  change (Int.signed (Int.repr 1)) with 1%Z.
+  set (b' := Int64.unsigned b) in *.
+  pose proof (Int64.unsigned_range b) as RANGE.
+  change Int64.modulus with 18446744073709551616%Z in RANGE.                                              
+  assert(1 <= IZR b' <= 18446744073709551616)%R as RANGE'.
+  { split; apply IZR_le; lia.
+  }
+
+  assert (-16777216 <= 1 <= 16777216)%Z as SILLY by lia.
+  destruct (BofZ_exact 24 128 re re 1 SILLY) as (C0R & C0F & _).
+  clear SILLY.
+  set (one_s := (BofZ 24 128 re re 1)) in *.
+  
+  pose proof (BofZ_correct 24 128 re re b') as C1.
+  cbn in C1.
+  rewrite Rlt_bool_true in C1; cycle 1.
+  { clear C1.
+    eapply (Rabs_relax (IZR 18446744073709551616)).
+    lra.
+    set (b'' := IZR b') in *.
+    gappa.
+  }
+  destruct C1 as (C1R & C1F & _).
+  set (b_s :=  (BofZ 24 128 re re b')) in *.
+
+  assert(1 <= B2R 24 128 b_s <= 18446744073709551616)%R as b_s_RANGE.
+  { rewrite C1R.
+    gappa.
+  }
+  assert(B2R 24 128 b_s <> 0)%R as b_s_NONZ by lra.
+  
+  pose proof (Bdiv_correct 24 128 re re Float32.binop_nan mode_NE one_s b_s b_s_NONZ) as C2.
+  rewrite Rlt_bool_true in C2; cycle 1.
+  { clear C2.
+    apply Rabs_relax with (b := 1%R).
+    { cbn; lra. }
+    rewrite C0R.
+    set (r_b_s := B2R 24 128 b_s) in *.
+    cbn.
+    gappa.
+  }
+  
+  destruct C2 as (C2R & C2F & _).
+  set (invb_s := (Bdiv 24 128 re re Float32.binop_nan mode_NE one_s b_s)) in *.
+  rewrite C0F in C2F.
+
+  assert ((1/18446744073709551616 <= B2R 24 128 invb_s <= 1)%R) as invb_s_RANGE.
+  { rewrite C2R.
+    set (r_b_s := B2R 24 128 b_s) in *.
+    rewrite C0R.
+    cbn.
+    gappa.
+  }
+  
+  pose proof (Bconv_correct 24 128 53 1024 re re Float.of_single_nan mode_NE invb_s C2F) as C3.
+  rewrite Rlt_bool_true in C3; cycle 1.
+  { clear C3.
+    set (r_invb_s := (B2R 24 128 invb_s)) in *.
+    apply Rabs_relax with (b := 1%R).
+    { replace 1%R with (bpow radix2 0)%R by reflexivity.
+      apply bpow_lt.
+      lia.
+    }
+    cbn.
+    gappa.
+  }
+  
+  destruct C3 as (C3R & C3F & _).
+  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.
+  { 
+    rewrite C3R.
+    set (r_invb_s := B2R 24 128 invb_s) in *.
+    cbn.
+    gappa.
+  }
+
+  split. assumption.
+  unfold rough_approx_inv_thresh.
+  rewrite C3R.
+  rewrite C2R.
+  rewrite C1R.
+  rewrite C0R.
+  cbn.
+  gappa.
+Qed.
-- 
cgit