some progress

1 files changed, 75 insertions, 0 deletions

diff --git a/kvx/FPDivision64.v b/kvx/FPDivision64.v
index f9557335..767e7ae0 100644
--- a/kvx/FPDivision64.v
+++ b/kvx/FPDivision64.v
@@ -952,3 +952,78 @@ Proof.
   set (delta2 := (q_r * IZR b' - IZR a')%R) in *.
   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.
+  { admit. }
+  apply Zdiv_unique with (b := (a - q*b)%Z).
+  ring.
+  apply range_up_le.
+  rewrite minus_IZR.
+  rewrite mult_IZR.
+  fold a' b' q'.
+  
+  
+  Search (Rabs _ * Rabs _)%R.
+Definition step2_real_div_long a b :=
+  Val.mulf (Val.maketotal (Val.floatoflong a)) (step1_real_inv_longu b).
+
+Definition step2_div_long a b :=
+  Val.maketotal (Val.longuoffloat_ne (step2_real_div_long a)).