From ba7a3de74288f8721699291d8fc2464c452498e1 Mon Sep 17 00:00:00 2001
From: David Monniaux <David.Monniaux@univ-grenoble-alpes.fr>
Date: Fri, 7 Jan 2022 15:02:58 +0100
Subject: progress in proofs

---
 kvx/FPDivision64.v | 183 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 183 insertions(+)
 create mode 100644 kvx/FPDivision64.v

(limited to 'kvx/FPDivision64.v')

diff --git a/kvx/FPDivision64.v b/kvx/FPDivision64.v
new file mode 100644
index 00000000..33a980c6
--- /dev/null
+++ b/kvx/FPDivision64.v
@@ -0,0 +1,183 @@
+From Flocq Require Import Core Digits Operations Round Bracket Sterbenz
+                          Binary Round_odd Bits.
+Require Archi.
+Require Import Coqlib.
+Require Import Compopts.
+Require Import AST.
+Require Import Integers.
+Require Import Floats.
+Require Import Op.
+Require Import CminorSel.
+Require Import OpHelpers.
+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.
+
+From Gappa Require Import Gappa_tactic.
+
+Definition approx_inv_longu b :=
+  let invb_s := ExtValues.invfs (Val.maketotal (Val.singleoflongu b)) in
+  let invb_d := Val.floatofsingle invb_s in
+  let b_d := Val.maketotal (Val.floatoflongu b) in
+  let one := Vfloat (ExtFloat.one) in
+  let alpha := ExtValues.fmsubf one invb_d b_d in
+  ExtValues.fmaddf invb_d alpha invb_d.
+
+Definition approx_inv_thresh := (1/70368744177664)%R.
+
+Lemma Rabs_relax:
+  forall b b' (INEQ : (b < b')%R) x,
+    (-b <= x <= b)%R -> (Rabs x < b')%R.
+Proof.
+  intros.
+  apply Rabs_lt.
+  lra.
+Qed.
+
+Theorem approx_inv_longu_correct :
+  forall b,
+    (1 <= Int64.unsigned b <= 9223372036854775808)%Z ->
+    exists (f : float),
+      (approx_inv_longu (Vlong b)) = Vfloat f /\
+      is_finite _ _ f = true /\ (Rabs((B2R _ _ f) - (1 / IZR (Int64.unsigned b))) <= approx_inv_thresh)%R.
+Proof.
+  intros b NONZ.
+  unfold 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 *.
+  assert(1 <= IZR b' <= 9223372036854775808)%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 9223372036854775808)).
+    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 <= 9223372036854775808)%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/9223372036854775808 <= 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/9223372036854775808 <= 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.
+  }
+
+  pose proof (is_finite_Bopp 53 1024 Float.neg_nan invb_d) as opp_finite.
+  rewrite C3F in opp_finite.
+
+  pose proof (BofZ_correct 53 1024 re re 1) as C4.
+  rewrite Rlt_bool_true in C4; cycle 1.
+  { clear C4.
+    cbn.
+    eapply (Rabs_relax (IZR 9223372036854775808)).
+    lra.
+    set (b'' := IZR b') in *.
+    gappa.
+  }
+  destruct C4 as (C4R & C4F & _).
+  
+  pose proof (BofZ_correct 53 1024 re re b') as C5.
+  cbn in C5.
+  rewrite Rlt_bool_true in C5; cycle 1.
+  { clear C5.
+    eapply (Rabs_relax (IZR 9223372036854775808)).
+    lra.
+    set (b'' := IZR b') in *.
+    gappa.
+  }
+  destruct C5 as (C5R & C5F & _).
+  set (b_d :=  (BofZ 53 1024 re re b')) in *.
+  
+  assert(1 <= B2R 53 1024 b_d <= 9223372036854775808)%R as b_d_RANGE.
+  { rewrite C5R.
+    gappa.
+  }
+
+  pose proof (Bfma_correct 53 1024 re re Float.fma_nan mode_NE
+          (Bopp 53 1024 Float.neg_nan invb_d) (BofZ 53 1024 re re b')
+          (BofZ 53 1024 re re 1) opp_finite C5F C4F) as C6.
+  rewrite Rlt_bool_true in C6; cycle 1.
+  { clear C6.
+    rewrite C4R.
+    rewrite B2R_Bopp.
+    cbn.
+    eapply (Rabs_relax (IZR 9223372036854775808)).
+    { lra. }
+    fold invb_d.
+    fold b_d.
+    set (r_invb_d := B2R 53 1024 invb_d) in *.
+    set (r_b_d := B2R 53 1024 b_d) in *.
+    
+    unfold invb_d.
+    unfold invb_s.
+    gappa.
+  }
+Admitted.
-- 
cgit