some progress?

1 files changed, 185 insertions, 126 deletions

diff --git a/kvx/FPDivision64.v b/kvx/FPDivision64.v
index 3685f601..69e49f29 100644
--- a/kvx/FPDivision64.v
+++ b/kvx/FPDivision64.v
@@ -1064,96 +1064,6 @@ Proof.
   gappa.
 Qed.
 
-(*
-Lemma range_up_le :
-  forall a x b,
-    (IZR a <= IZR x <= (IZR b) - 1)%R ->
-    (a <= x < b)%Z.
-Proof.
-  intros until b. intro RANGE.
-  split.
-  { apply le_IZR. lra. }
-  assert (x <= b-1)%Z.
-  { apply le_IZR. rewrite minus_IZR. lra. }
-  lia.
-Qed.
- *)
-
-(*
-Lemma find_quotient:
-  forall (a b : Z)
-         (b_POS : (b > 0)%Z)
-         (invb : R)
-         (eps : R)
-         (invb_EPS : (Rabs (invb * IZR b - 1) <= eps)%R)
-         (SMALL : (IZR (Z.abs a)*eps < IZR b / 2)%R),
-    (a / b)%Z =
-      let q := ZnearestE ((IZR a) * invb) in
-      if (b*q >? a)%Z
-      then (q-1)%Z
-      else q.
-Proof.
-  intros.
-  rewrite <- Rabs_Zabs in SMALL.
-  set (q := ZnearestE (IZR a * invb)%R).
-  cbn zeta.
-  set (b' := IZR b) in *.
-  set (a' := IZR a) in *.
-  assert (1 <= b')%R as b_POS'.
-  { apply IZR_le.
-    lia.
-  }
-  
-  pose proof (Znearest_imp2 (fun x : Z => negb (Z.even x)) (a' * invb)) as NEAR.
-  fold q in NEAR.
-  set (q' := IZR q) in *.
-  assert ((Rabs a') * Rabs (invb * b' - 1) <= (Rabs a') * eps)%R as S1.
-  { apply Rmult_le_compat_l.
-    apply Rabs_pos.
-    assumption. }
-  rewrite <- Rabs_mult in S1.
-  assert ((Rabs b') * Rabs (q' - a' * invb) <= (Rabs b') / 2)%R as S2.
-  { apply Rmult_le_compat_l.
-    apply Rabs_pos.
-    assumption. }
-  rewrite <- Rabs_mult in S2.
-  rewrite (Rabs_right b') in S2 by lra.
-  pose proof (Rabs_triang (a' * (invb * b' - 1))
-                          (b' * (q' - a' * invb))) as TRIANGLE.
-  replace (a' * (invb * b' - 1) + b' * (q' - a' * invb))%R with
-    (b' * q' - a')%R in TRIANGLE by ring.
-  assert (Rabs (b' * q' - a') < b')%R as DELTA by lra.
-  
-  pose proof (Zgt_cases (b * q) a)%Z as CASE.
-  destruct (_ >? _)%Z.
-  { unfold b', q', a' in DELTA.
-    rewrite <- mult_IZR in DELTA.
-    rewrite <- minus_IZR in DELTA.
-    rewrite Rabs_Zabs in DELTA.
-    apply lt_IZR in DELTA.
-    rewrite Z.abs_eq in DELTA by lia.
-    
-    apply Zdiv_unique with (b := (a - (q-1)*b)%Z).
-    ring.
-    split; lia.
-  }
-  
-  rewrite <- Rabs_Ropp in DELTA.
-  unfold b', q', a' in DELTA.
-  rewrite <- mult_IZR in DELTA.
-  rewrite <- minus_IZR in DELTA.
-  rewrite <- opp_IZR in DELTA.
-  rewrite Rabs_Zabs in DELTA.
-  apply lt_IZR in DELTA.
-  rewrite Z.abs_eq in DELTA by lia.
-
-  apply Zdiv_unique with (b := (a - q*b)%Z).
-  ring.
-  
-  split; lia.
-Qed.
- *)
-
 Lemma find_quotient:
   forall (a b : Z)
          (b_POS : (0 < b)%Z)
@@ -1754,48 +1664,197 @@ Proof.
   rewrite zlt_false by lia.
   reflexivity.
 Qed.
-  
-Lemma twostep_div_longu_bigb_correct :
-    forall a b
-           (b_BIG : ((Int64.unsigned b) > smallb_thresh)%Z)
-           (b_NOT_TOO_BIG : ((Int64.unsigned b) <= Int64.max_signed)%Z),
-          (twostep_div_longu (Vlong a) (Vlong b)) =
-            (Val.maketotal (Val.divlu (Vlong a) (Vlong b))).
+
+Definition step2_real_div_longu a b :=
+  Val.mulf (Val.maketotal (Val.floatoflongu a)) (approx_inv_longu b).
+
+Definition step2_div_longu' a b :=
+  Val.maketotal (Val.longuoffloat_ne (step2_real_div_longu a b)).
+
+Definition step2_div_longu a b :=
+  let q := step2_div_longu' a b in
+  Val.select (Val.cmpl_bool Cgt (Val.subl (Val.mull q b) a) (Vlong Int64.zero))
+             (Val.addl q (Vlong Int64.mone)) q Tlong.
+
+Lemma step2_real_div_longu_bigb_correct :
+    forall (a b : int64)
+           (b_BIG : ((Int64.unsigned b) > smallb_thresh)%Z),
+    exists (q : float),
+      (step2_real_div_longu (Vlong a) (Vlong b)) = Vfloat q /\
+        (Rabs ((B2R _ _ q) - (IZR (Int64.unsigned a)) / (IZR (Int64.unsigned b))) <= (32767/65536))%R /\
+      is_finite _ _ q = true.
 Proof.
   intros.
-  unfold twostep_div_longu.
-  set (b' := Int64.unsigned b) in *.
+  unfold step2_real_div_longu.
+  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_longu.
+  unfold Float.mul, Float.of_longu.
+  set (re := (@eq_refl Datatypes.comparison Lt)) in *.
+  pose proof (Int64.unsigned_range b) as b_RANGE.
+  pose proof (Int64.unsigned_range a) as a_RANGE.
+  change Int64.modulus with 18446744073709551616%Z in *.
   set (a' := Int64.unsigned a) in *.
-Admitted.
+  set (b' := Int64.unsigned b) in *.
+  assert (IZR (1 + smallb_thresh) <= IZR b' <= 18446744073709551615)%R as b_RANGE'.
+  { split; apply IZR_le; lia.
+  }
+  assert(0 <= IZR a' <= 18446744073709551615)%R as a_RANGE'.
+  { split; apply IZR_le; lia.
+  }
+  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_div_longu_bigb_correct :
+  forall a b
+           (b_BIG : ((Int64.unsigned b) > smallb_thresh)%Z)
+           (b_NOT_TOO_BIG : ((Int64.unsigned b) <= Int64.max_signed)%Z),
+      step2_div_longu (Vlong a) (Vlong b) = Vlong (Int64.repr (Int64.unsigned a / Int64.unsigned b))%Z.
+Proof.
+  intros.
+  pose proof (Int64.unsigned_range b) as b_RANGE.
+  pose proof (Int64.unsigned_range a) as a_RANGE.
+  change Int64.modulus with 18446744073709551616%Z in *.
+  set (a' := (Int64.unsigned 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(0 <= IZR a' <= 18446744073709551615)%R as a_RANGE'.
+  { split; apply IZR_le; lia.
+  }
+  unfold smallb_thresh in *; cbn in b_RANGE'.
   assert (0 < b')%Z as b_NOT0 by lia.
-  assert (b' <= Int64.max_signed)%Z as b_NOTBIG.
-  { change Int64.max_signed with (9223372036854775807)%Z.
-    unfold smallb_thresh in b_RANGE.
-    lia.
+
+  destruct (step2_real_div_longu_bigb_correct a b b_BIG) as (q & C1R & C1E & C1F).
+  fold a' b' in C1E.
+
+  assert ((0 <=? ZnearestE (B2R 53 1024 q))=true)%Z as q_LOW.
+  { apply Zle_imp_le_bool.
+    set (q' :=  B2R 53 1024 q) in *.
+    assert (-32767 / 65536 <= q')%R as LOWROUND.
+    { replace q' with (IZR a' / IZR b' + (q' - IZR a' / IZR b'))%R by ring.
+      gappa.
+    }
+    destruct (Rcase_abs q').
+    { replace (ZnearestE q') with 0%Z. lia.
+      symmetry.
+      apply Znearest_imp.
+      apply Rabs_lt.
+      split; lra.
+    }
+    apply Znearest_lub.
+    lra.
   }
+  assert ((ZnearestE (B2R 53 1024 q) <=? Int64.max_unsigned)=true)%Z as q_HIGH.
+  { apply Zle_imp_le_bool.
+    change Int64.max_unsigned with (18446744073709551615)%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.
+    gappa.
+  }
+  
+  unfold step2_div_longu, step2_div_longu'.
+  rewrite C1R.
   cbn.
-  rewrite (step2_div_long_bigb_correct (Int64.sub a (Int64.mul b q1)) b r1_SMALL b_NOT0 b_NOTBIG).
-  unfold Int64.add, Int64.sub, Int64.mul, Int64.divu.
-  fold q1' b' a'.
-  rewrite <- unsigned_repr_sub.
-  rewrite <- unsigned_repr_add.
-  rewrite Int64.signed_repr ; cycle 1.
-  {
-    change Int64.min_signed with (-9223372036854775808)%Z.
-    change Int64.max_signed with (9223372036854775807)%Z.
-    unfold smallb_approx_range in *.
-    lia.
+  unfold Float.to_longu_ne.
+  rewrite (ZofB_ne_range_correct _ _ q _ _).
+  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'.
+  rewrite Int64.signed_zero.
+  set (q'' :=  (ZnearestE q')) in *.
+  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.
   }
-  rewrite Z.add_comm.
-  rewrite <- Z.div_add by lia.
-  replace (a' - b' * q1' + q1' * b')%Z with a' by ring.
-  rewrite Int64.eq_false ; cycle 1.
-  { intro Z. unfold b' in b_NOT0. rewrite Z in b_NOT0.
-    rewrite Int64.unsigned_zero in b_NOT0.
-    lia.
+  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.
+  assert ((Int64.signed (Int64.sub (Int64.mul (Int64.repr q'') b) a))
+          = q''*b'-a')%Z as SIMPL.
+  { admit. }
+  rewrite SIMPL.
+  destruct (Z_lt_dec a' (b' * q'')) as [LT | GE].
+  { replace (b' * q'' >? a')%Z with true in QUOTIENT by lia.
+    rewrite zlt_true by lia.
+    replace q'' with (1 + (a' / b'))%Z by lia.
+    unfold Int64.add.
+    apply Int64.eqm_samerepr.
+    apply Int64.eqm_trans with (y := ((1 + a' / b') + (-1))%Z).
+    { apply Int64.eqm_add.
+      apply Int64.eqm_sym.
+      apply Int64.eqm_unsigned_repr.
+      rewrite Int64.unsigned_mone.
+      replace (-1)%Z with (0 - 1)%Z by ring.
+      apply Int64.eqm_add.
+      exists 1%Z.
+      lia.
+      apply Int64.eqm_refl.
+    }
+    replace (1 + a' / b' + -1)%Z with (a'/b')%Z by ring.
+    apply Int64.eqm_refl.
   }
-  reflexivity.
-Qed.
-*)
+  rewrite zlt_false by lia.
+  replace (b' * q'' >? a')%Z with false in QUOTIENT by lia.
+  congruence.
+Admitted.