find_quotient

1 files changed, 64 insertions, 0 deletions

diff --git a/kvx/FPDivision64.v b/kvx/FPDivision64.v
index 8f0d2fa1..b5052917 100644
--- a/kvx/FPDivision64.v
+++ b/kvx/FPDivision64.v
@@ -968,6 +968,7 @@ Proof.
 Qed.
  *)
 
+(*
 Lemma find_quotient:
   forall (a b : Z)
          (b_POS : (b > 0)%Z)
@@ -1040,6 +1041,69 @@ Proof.
   
   split; lia.
 Qed.
+ *)
+
+Lemma find_quotient:
+  forall (a b : Z)
+         (b_POS : (b > 0)%Z)
+         (invb : R)
+         (qr : R)
+         (GAP : (Rabs (IZR a / IZR b - qr) < /2)%R),
+    (a / b)%Z =
+      let q := ZnearestE qr in
+      if (b*q >? a)%Z
+      then (q-1)%Z
+      else q.
+Proof.
+  intros.
+  set (q := ZnearestE qr).
+  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)) qr) as ROUND.
+  fold q in ROUND.
+  set (q' := IZR q) in *.
+  
+  pose proof (Rabs_triang (a' / b' - qr)
+                          (qr - q'))%R as TRIANGLE.
+  replace ((a' / b' - qr) +  (qr - q'))%R with
+    (a' / b' - q')%R in TRIANGLE by ring.
+  rewrite <- Rabs_Ropp in ROUND.
+  replace (- (q' - qr))%R with (qr - q')%R in ROUND by ring.
+  assert (Z.abs (a - b*q) < b)%Z as DELTA.
+  { apply lt_IZR.
+    rewrite <- Rabs_Zabs.
+    rewrite minus_IZR.
+    rewrite mult_IZR.
+    fold a' q' b'.
+    apply Rmult_lt_reg_r with (r := (/b')%R).
+    { apply Rinv_0_lt_compat. lra. }
+    rewrite Rinv_r by lra.
+    replace (/ b')%R with (/ Rabs(b'))%R ; cycle 1.
+    { f_equal.
+      apply Rabs_pos_eq. lra. }
+    rewrite <- Rabs_Rinv by lra.
+    rewrite <- Rabs_mult.
+    replace ((a' - b' * q') * / b')%R with (a'/b' - q')%R by (field ; lra).
+    lra.
+  }
+   
+  pose proof (Zgt_cases (b * q) a)%Z as CASE.
+  destruct (_ >? _)%Z.
+  { apply Zdiv_unique with (b := (a - (q-1)*b)%Z).
+    ring.
+    split; lia.
+  }
+
+  apply Zdiv_unique with (b := (a - q*b)%Z).
+  ring.
+  split; lia.
+Qed.
 
 Definition step2_real_div_long a b :=
   Val.mulf (Val.maketotal (Val.floatoflong a)) (step1_real_inv_longu b).