From 30c2e58daf62d53530505bea2401e646df79bd31 Mon Sep 17 00:00:00 2001
From: David Monniaux <David.Monniaux@univ-grenoble-alpes.fr>
Date: Mon, 24 Jan 2022 21:06:06 +0100
Subject: progress

---
 kvx/FPDivision64.v | 190 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 190 insertions(+)

(limited to 'kvx/FPDivision64.v')

diff --git a/kvx/FPDivision64.v b/kvx/FPDivision64.v
index 758061b5..580f4dd9 100644
--- a/kvx/FPDivision64.v
+++ b/kvx/FPDivision64.v
@@ -1467,3 +1467,193 @@ Proof.
   }
   reflexivity.
 Qed.
+
+
+Lemma step2_real_div_long_bigb_correct :
+    forall (a b : int64)
+           (b_BIG : ((Int64.unsigned b) > smallb_thresh)%Z),
+    exists (q : float),
+      (step2_real_div_long (Vlong a) (Vlong b)) = Vfloat q /\
+        (Rabs ((B2R _ _ q) - (IZR (Int64.signed a)) / (IZR (Int64.unsigned b))) <= (32767/65536))%R /\
+      is_finite _ _ q = true.
+Proof.
+  intros.
+  unfold step2_real_div_long.
+  assert (0 < Int64.unsigned b)%Z as b_NOT0 by (unfold smallb_thresh in *; lia).
+  destruct (approx_inv_longu_correct_rel b b_NOT0) as (f & C0E & C0F & C0R).
+  rewrite C0E.
+  econstructor.
+  split. reflexivity.
+  Local Transparent Float.of_long.
+  unfold Float.mul, Float.of_long.
+  set (re := (@eq_refl Datatypes.comparison Lt)) in *.
+  pose proof (Int64.unsigned_range b) as b_RANGE.
+  change Int64.modulus with 18446744073709551616%Z in b_RANGE.
+  pose proof (Int64.signed_range a) as a_RANGE.
+  set (a' := Int64.signed a) in *.
+  set (b' := Int64.unsigned b) in *.
+  assert (IZR (1 + smallb_thresh) <= IZR b' <= 18446744073709551615)%R as b_RANGE'.
+  { split; apply IZR_le; lia.
+  }
+  assert(IZR Int64.min_signed <= IZR a' <= IZR Int64.max_signed)%R as a_RANGE'.
+  { split; apply IZR_le; lia.
+  }
+  change Int64.min_signed with (-9223372036854775808)%Z in a_RANGE'.
+  change Int64.max_signed with (9223372036854775807)%Z in a_RANGE'.
+  pose proof (BofZ_correct 53 1024 re re a') as C1.
+  rewrite Rlt_bool_true in C1 ; cycle 1.
+  { clear C1.
+    apply Rabs_relax with (b := bpow radix2 64).
+    { apply bpow_lt; lia. }
+    cbn.
+    gappa.
+  }
+  cbn in C1.
+  destruct C1 as (C1R & C1F & C1S).
+
+  unfold smallb_thresh in b_RANGE'; cbn in b_RANGE'.
+
+  pose proof (Bmult_correct 53 1024 re re Float.binop_nan mode_NE (BofZ 53 1024 re re a') f) as C2.
+  rewrite Rlt_bool_true in C2 ; cycle 1.
+  { clear C2.
+    apply Rabs_relax with (b := bpow radix2 53).
+    { apply bpow_lt. lia. }
+    cbn.
+    rewrite C1R.
+    unfold approx_inv_rel_thresh in C0R.
+    replace (B2R 53 1024 f) with
+      ((1/IZR b') * ((IZR b' * B2R 53 1024 f - 1) + 1))%R ; cycle 1.
+    { field.  lra. }
+    gappa.
+  }
+  rewrite C0F in C2.
+  rewrite C1R in C2.
+  rewrite C1F in C2.
+  rewrite C1S in C2.
+  cbn in C2.
+  destruct C2 as (C2R & C2F & _).
+  split.
+  2: exact C2F.
+  rewrite C2R.
+  set (f' := (B2R 53 1024 f)) in *.  
+  replace (rd(rd (IZR a') * f') - IZR a' / IZR b')%R with
+    ((rd(rd (IZR a') * f') - IZR a' * f') + IZR a' / IZR b' * (IZR b' * f' - 1))%R ; cycle 1.
+  { field. lra. }
+  unfold approx_inv_rel_thresh in *.
+  gappa.
+Qed.
+
+Lemma step2_real_div_long_bigb_correct :
+    forall (a b : int64)
+           (b_BIG : ((Int64.unsigned b) > smallb_thresh)%Z),
+    exists (q : float),
+      (step2_real_div_long (Vlong a) (Vlong b)) = Vfloat q /\
+        (Rabs ((B2R _ _ q) - (IZR (Int64.signed a)) / (IZR (Int64.unsigned b))) <= (32767/65536))%R /\
+      is_finite _ _ q = true.
+Proof.
+  intros.
+  pose proof (Int64.unsigned_range b) as b_RANGE.
+  change Int64.modulus with 18446744073709551616%Z in b_RANGE.
+  set (a' := (Int64.signed a)) in *.
+  set (b' := (Int64.unsigned b)) in *.
+  assert (IZR (1 + smallb_thresh) <= IZR b' <= 18446744073709551615)%R as b_RANGE'.
+  { split; apply IZR_le; lia.
+  } 
+  destruct (step2_real_div_long_smalla_correct a b a_SMALL b_NOT0) as (q & C1R & C1E & C1F).
+  fold a' b' in C1E.
+  assert ((Int64.min_signed <=? ZnearestE (B2R 53 1024 q))=true)%Z as q_LOW.
+  { apply Zle_imp_le_bool.
+    change Int64.min_signed with (-9223372036854775808)%Z.
+    apply Znearest_lub.
+    set (q' :=  B2R 53 1024 q) in *.
+    replace q' with (IZR a' / IZR b' + (q' - IZR a' / IZR b'))%R by ring.
+    unfold smalla_thresh in a_RANGE'.
+    gappa.
+  }
+  assert ((ZnearestE (B2R 53 1024 q) <=? Int64.max_signed)=true)%Z as q_HIGH.
+  { apply Zle_imp_le_bool.
+    change Int64.max_signed with (9223372036854775807)%Z.
+    apply Znearest_glb.
+    set (q' :=  B2R 53 1024 q) in *.
+    replace q' with (IZR a' / IZR b' + (q' - IZR a' / IZR b'))%R by ring.
+    unfold smalla_thresh in a_RANGE'.
+    gappa.
+  }
+  unfold step2_div_long, step2_div_long'.
+  rewrite C1R.
+  cbn.
+  unfold Float.to_long_ne.
+  rewrite (ZofB_ne_range_correct _ _ q Int64.min_signed Int64.max_signed).
+  rewrite C1F.
+  rewrite q_LOW.
+  rewrite q_HIGH.
+  cbn.
+  rewrite normalize_ite.
+  cbn.
+  rewrite <- (function_ite Vlong).
+  f_equal.
+  unfold Int64.lt.
+  set (q' :=  B2R 53 1024 q) in *.
+  fold a'.
+  assert (Int64.signed (Int64.repr (ZnearestE q')) = ZnearestE q') as q_SIGNED.
+  { apply Int64.signed_repr.
+    split; lia.
+  }
+  rewrite Int64.add_signed.
+  rewrite q_SIGNED.
+  rewrite Int64.signed_mone.
+  rewrite Int64.signed_zero.
+  rewrite <- (function_ite Int64.repr).
+  f_equal.
+  replace  (ZnearestE q' + -1)%Z with (ZnearestE q' - 1)%Z by ring.
+  
+  set (q'' :=  (ZnearestE q')) in *.
+  Check Int64.sub_signed.
+  fold a'.
+  rewrite int64_mul_signed_unsigned.
+  rewrite q_SIGNED.
+  fold b'.
+
+  rewrite Int64.sub_signed.
+  fold a'.
+  rewrite signed_repr_sub.
+
+  assert ((Rabs (IZR a' / IZR b' - q') < / 2)%R) as HALF.
+  { replace (IZR a' / IZR b' - q')%R with
+      (-(q' - IZR a' / IZR b'))%R by ring.
+    rewrite Rabs_Ropp.
+    lra.
+  }
+  pose proof (find_quotient a' b' b_NOT0 q' HALF) as QUOTIENT.
+  fold q'' in QUOTIENT.
+  cbn zeta in QUOTIENT.
+
+  assert (b' <> 0)%Z as NONZ by lia.
+  pose proof (Zmod_eq_full a' b' NONZ) as MOD.
+  assert (b' > 0)%Z as b_GT0 by lia.
+  pose proof (Z_mod_lt a' b' b_GT0) as MOD_LT.
+  destruct (Z_lt_dec a' (b' * q'')) as [LT | GE].
+  { replace (b' * q'' >? a')%Z with true in QUOTIENT by lia.
+    replace q'' with (1 + (a' / b'))%Z by lia.
+    replace ((1 + a' / b') * b' - a')%Z
+      with (-(a' - a' / b' * b')+b')%Z by ring.
+    rewrite <- MOD.
+    rewrite Int64.signed_repr; cycle 1.
+    { change Int64.min_signed with (-9223372036854775808)%Z in *.
+      change Int64.max_signed with (9223372036854775807)%Z in *.
+      lia.
+    }
+    rewrite zlt_true by lia.
+    ring.
+  }
+  replace (b' * q'' >? a')%Z with false in QUOTIENT by lia.
+  rewrite <- QUOTIENT.
+  replace (a' / b' * b' - a')%Z with (-(a' - a' / b' * b'))%Z by ring.
+  rewrite <- MOD.
+  rewrite Int64.signed_repr ; cycle 1.
+  { change Int64.min_signed with (-9223372036854775808)%Z in *.
+    change Int64.max_signed with (9223372036854775807)%Z in *.
+    lia.
+  }
+  rewrite zlt_false by lia.
+  reflexivity.
-- 
cgit