progress in proofs

1 files changed, 30 insertions, 224 deletions

diff --git a/kvx/FPDivision64.v b/kvx/FPDivision64.v
index f4d727be..e7ba56e3 100644
--- a/kvx/FPDivision64.v
+++ b/kvx/FPDivision64.v
@@ -283,27 +283,27 @@ Proof.
   gappa.
 Qed.
 
-Definition rough_approx_inv_longu b :=
+Definition step1_real_inv_longu b :=
   let invb_s := ExtValues.invfs (Val.maketotal (Val.singleoflongu b)) in
   Val.floatofsingle invb_s.
 
-Definition rough_approx_inv_thresh := (3/33554432)%R.
+Definition step1_real_inv_thresh := (3/33554432)%R.
 (* 8.94069671630859e-8 *)
 
 Local Notation "'rd'" := (round radix2 (FLT_exp (-1074) 53) ZnearestE).
 Local Notation "'rs'" := (round radix2 (FLT_exp (-149) 24) ZnearestE).
 
-Theorem rough_approx_inv_longu_correct :
+Theorem step1_real_inv_longu_correct :
   forall b,
     (0 < Int64.unsigned b)%Z ->
     exists (f : float),
-      (rough_approx_inv_longu (Vlong b)) = Vfloat f /\
+      (step1_real_inv_longu (Vlong b)) = Vfloat f /\
         (B2R _ _ f) = (rd (rs (1 / rs (IZR (Int64.unsigned b))))) /\
         is_finite _ _ f = true /\
         Bsign _ _ f = false.
 Proof.
   intros b NONZ.
-  unfold rough_approx_inv_longu.
+  unfold step1_real_inv_longu.
   cbn.
   econstructor.
   split.
@@ -402,18 +402,18 @@ Proof.
   auto.
 Qed.
 
-Theorem rough_approx_inv_longu_correct1 :
+Theorem step1_real_inv_longu_correct1 :
   forall b,
     (Int64.unsigned b = 1%Z) ->
     exists f,
-    (rough_approx_inv_longu (Vlong b)) = Vfloat f /\
+    (step1_real_inv_longu (Vlong b)) = Vfloat f /\
         (B2R _ _ f) = 1%R /\
         is_finite _ _ f = true /\
         Bsign _ _ f = false.
 Proof.
   intros b EQ1.
   assert (0 < Int64.unsigned b)%Z as b_RANGE by lia.
-  destruct (rough_approx_inv_longu_correct b b_RANGE) as (f & C1E & C1R & C1F & C1S).
+  destruct (step1_real_inv_longu_correct b b_RANGE) as (f & C1E & C1R & C1F & C1S).
   rewrite EQ1 in C1R.
   exists f.
   repeat split; try assumption.
@@ -457,20 +457,20 @@ Proof.
   lia.
 Qed.
 
+Definition step1_real_div_longu a b :=
+  Val.mulf (Val.maketotal (Val.floatoflongu a)) (step1_real_inv_longu b).
 
-Definition rough_approx_div_longu a b :=
-  Val.mulf (Val.maketotal (Val.floatoflongu a)) (rough_approx_inv_longu b).
-Definition smallb_threshold := 4398046511104%Z.
-Definition rough_approx_div_longu_smallb_delta := 3917010173952%Z.
+Definition step1_div_longu a b :=
+  Val.maketotal (Val.longuoffloat (Val.mulf (Val.maketotal (Val.floatoflongu a)) (step1_real_inv_longu b))).
 
 Definition step1_real_quotient (a b : R) :=
              rd ((rd (a)) * (rd (rs (1 / rs (b))))).
   
-Theorem rough_approx_div_longu_correct:
+Theorem step1_real_div_longu_correct:
   forall a b,
     (1 < Int64.unsigned b)%Z ->
     exists (q : float),
-      (rough_approx_div_longu (Vlong a) (Vlong b)) = Vfloat q /\
+      (step1_real_div_longu (Vlong a) (Vlong b)) = Vfloat q /\
         (B2R _ _ q) = step1_real_quotient (IZR (Int64.unsigned a))
                                           (IZR (Int64.unsigned b)) /\
         is_finite _ _ q = true /\
@@ -478,8 +478,8 @@ Theorem rough_approx_div_longu_correct:
 Proof.
   intros a b b_NON01.
   assert (0 < Int64.unsigned b)%Z as b_NON0 by lia.
-  destruct (rough_approx_inv_longu_correct b b_NON0) as (f & C1E & C1R & C1F & C1S).
-  unfold rough_approx_div_longu.
+  destruct (step1_real_inv_longu_correct b b_NON0) as (f & C1E & C1R & C1F & C1S).
+  unfold step1_real_div_longu.
   rewrite C1E.
   cbn.
   set (b' := Int64.unsigned b) in *.
@@ -536,224 +536,30 @@ Proof.
   auto.
 Qed.
 
-Definition smallb_thresh := 4398046511104%Z.
+Definition smallb_thresh :=       4398046511104%Z.
 
-Definition smallb_approx_range := 1649280000000%R.
+Definition smallb_approx_range := 2200000000000%R.
 Lemma step1_smallb :
   forall a b
-    (a_RANGE : (0 <= a <= 18446744073709551615)%R)
+    (a_RANGE : (1 <= a <= 18446744073709551615)%R)
     (b_RANGE : (1 <= b <= IZR smallb_thresh)%R),
-    (Rabs((step1_real_quotient a b) - a/b) <= smallb_approx_range)%R.
+    (Rabs((step1_real_quotient a b) * b - a) <= smallb_approx_range)%R.
 Proof.
   intros.
   unfold smallb_thresh in b_RANGE.
   unfold smallb_approx_range.
-  replace (step1_real_quotient a b - a/b)%R with
-    ((rd ((rd (a)) * (rd (rs (1 / rs (b)))))
-             - ((rd (a)) * (rd (rs (1 / rs (b)))))) +
-	   (rd(a)) * (rd (rs (1 / rs (b))) - 1/b) +
-       (rd(a) - a)/b)%R; cycle 1.
-  { unfold step1_real_quotient.
+  unfold step1_real_quotient.
+  set (q := ((rd (a)) * (rd (rs (1 / rs (b)))))%R) in *.
+  replace ((rd q) *b - a)%R with
+     ((rd(q)-q)/q * rd(a) * (1 + (rd (rs (1 / rs (b))) - 1/b)/(1/b)) +
+  (rd (a)) * ((rd (rs (1 / rs (b))) - 1 / b) / (1/b)) +
+        (rd(a) - a))%R; cycle 1.
+  { unfold q.
     field.
-    lra.
-  }
-  gappa.
-Qed.
-
-(*
-Definition rough_approx_div_thresh := 1683627180032%R.
-
-Theorem rough_approx_div_longu_correct :
-  forall a b,
-    (2 <= Int64.unsigned b)%Z ->
-    exists (q : float),
-      (rough_approx_div_longu (Vlong a) (Vlong b)) = Vfloat q /\
-        (Rabs (B2R _ _ q - (IZR (Int64.unsigned a)) / (IZR (Int64.unsigned b)))
-         <= rough_approx_div_thresh)%R /\
-        is_finite _ _ q = true /\
-        Bsign _ _ q = false.
-Proof.
-  intros a b b_NON01.
-  unfold rough_approx_div_longu.
-  assert (0 < Int64.unsigned b)%Z as b_NONZ by lia.
-  destruct (rough_approx_inv_longu_correct b b_NONZ) as (f & C1R & C1F & C1D & C1S).
-  rewrite C1R.
-  cbn.
-
-  Local Transparent Float.to_longu Float.mul.
-  unfold Float.to_longu_ne, Float.mul, Float.of_longu.
-  set (re := (@eq_refl Datatypes.comparison Lt)).
-
-  pose proof (Int64.unsigned_range a) as a_RANGE.
-  change Int64.modulus with 18446744073709551616%Z in a_RANGE.
-  set (a' := Int64.unsigned a) in *.
-  assert (0 <= IZR a' <= 18446744073709551615)%R as IZR_a_RANGE.
-  { split ; apply IZR_le; lia.
-  }
-
-  pose proof (Int64.unsigned_range b) as b_RANGE.
-  change Int64.modulus with 18446744073709551616%Z in b_RANGE.
-  set (b' := Int64.unsigned b) in *.
-  assert (2 <= IZR b' <= 18446744073709551615)%R as IZR_b_RANGE.
-  { split ; apply IZR_le; lia.
-  }
-
-  set (a_d := (round radix2 (FLT_exp (-1074) 53) ZnearestE (IZR a'))).
-  assert (0 <= a_d <= 18446744073709551616)%R as a_d_RANGE.
-  { unfold a_d.
-    gappa.
-  }
-
-  set (f_r := (B2R _ _ f)) in *.
-  assert (0 <= f_r <= 4194305/8388608)%R as f_RANGE.
-  { split.
-    { apply Bsign_false_nonneg. trivial .}
-    unfold rough_approx_inv_thresh in C1D.
-    replace f_r with ((f_r - 1 / IZR b') + 1/ IZR b')%R by (field; lra).
-    set (delta := (f_r - 1 / IZR b')%R) in *.
-    gappa.
-  }
-  
-  pose proof (BofZ_correct 53 1024 re re a') as C2.
-  rewrite Rlt_bool_true in C2; cycle 1.
-  { clear C2.
-    apply Rabs_relax with (b := bpow radix2 64).
-    { apply bpow_lt. lia. }
-    cbn.
-    gappa.
-  }
-  destruct C2 as (C2R & C2F & C2S).
-
-  pose proof (Bmult_correct  53 1024 re re Float.binop_nan mode_NE
-                             (BofZ 53 1024 re re a') f) as C3.
-  rewrite Rlt_bool_true in C3; cycle 1.
-  { clear C3.
-    rewrite C2R.
-    fold f_r.
-    
-    apply Rabs_relax with (b := bpow radix2 65).
-    { apply bpow_lt ; lia. }
-    cbn.
+    split. lra.
+    split. gappa.
     gappa.
   }
-  rewrite C1F in C3.
-  rewrite C2F in C3.
-  cbn in C3.
-  destruct C3 as (C3R & C3F & C3S).
-
-  econstructor. split. reflexivity.
-  rewrite C3R.
-  rewrite C2R.
-  rewrite C3S; cycle 1.
-  { apply is_finite_not_is_nan.
-    assumption. }
-  rewrite C2S.
-  rewrite C1S.
-  rewrite C3F.
-  rewrite Zlt_bool_false by lia.
-  split. 2: auto; fail.
-
-  cbn.
-  unfold rough_approx_div_thresh.
-  fold f_r.
-  set (rd := round radix2 (FLT_exp (-1074) 53) ZnearestE).
-  replace (rd (rd (IZR a') * f_r) - IZR a' / IZR b')%R with
-   ((rd (rd (IZR a') * f_r) - (rd (IZR a') * f_r))
-    + ((rd (IZR a') - IZR a') * f_r)
-    + (IZR a') * (f_r - 1 / IZR b'))%R by (field ; lra).
-
-  assert (Rabs (rd (IZR a') - IZR a') <= 1024)%R as DELTA1.
-  { unfold rd. gappa. }
-  set (delta1 := (rd (IZR a') - IZR a')%R) in *.
-  
-  assert (Rabs (rd (rd (IZR a') * f_r) - rd (IZR a') * f_r) <= 1024)%R as DELTA2.
-  { unfold rd. gappa. }
-  set (delta2 := (rd (rd (IZR a') * f_r) - rd (IZR a') * f_r)%R) in *.
-  unfold rough_approx_inv_thresh in C1D.
-  set (delta3 := (f_r - 1 / IZR b')%R) in *.
+  unfold q.
   gappa.
 Qed.
- *)
-
-Definition rough_approx_euclid_div_longu a b :=
-  Val.maketotal (Val.longuoffloat_ne (rough_approx_div_longu a b)).
-
-
-Theorem rough_approx_div_longu_smallb :
-  forall a b,
-    (2 <= Int64.unsigned b <= smallb_threshold)%Z ->
-    exists (qr : int64),
-      (rough_approx_euclid_div_longu (Vlong a) (Vlong b)) = Vlong qr /\
-        (Z.abs (Int64.unsigned a - Int64.unsigned b* Int64.unsigned qr) <= rough_approx_div_longu_smallb_delta)%Z /\
-        (0 <= Int64.unsigned qr <= 11529215046068469760)%Z.
-Proof.
-  intros a b b_RANGE.
-  unfold rough_approx_euclid_div_longu.
-  destruct (rough_approx_div_longu_correct a b) as (q & C1R & C1D & C1F & C1S).
-  lia.
-  rewrite C1R.
-  cbn.
-  set (qr := B2R 53 1024 q) in *.
-  pose proof (Int64.unsigned_range a) as a_RANGE.
-  change Int64.modulus with 18446744073709551616%Z in a_RANGE.
-  set (a' := Int64.unsigned a) in *.
-  assert (0 <= IZR a' <= 18446744073709551615)%R as IZR_a_RANGE.
-  { split ; apply IZR_le; lia.
-  }
-
-  set (b' := Int64.unsigned b) in *.
-  assert (2 <= IZR b' <= IZR smallb_threshold)%R as IZR_b_RANGE.
-  { split ; apply IZR_le; lia.
-  }
-
-  assert (0 <= qr <= 11529215046068469760)%R as q_RANGE.
-  { split.
-    { unfold qr.
-      apply Bsign_false_nonneg.
-      assumption.
-    }
-    unfold rough_approx_div_thresh in C1D.
-    replace qr with ((qr - IZR a' / IZR b') + (IZR a' / IZR b'))%R by (field ; lra).
-    set (delta := (qr - IZR a' / IZR b')%R) in *.
-    unfold smallb_threshold in *.
-    gappa.
-  }
-  unfold Float.to_longu_ne.
-  pose proof (ZofB_ne_range_correct _ _ q 0 Int64.max_unsigned) as C2.
-  rewrite C1F in C2.
-  cbn in C2.
-  fold qr in C2.
-  assert (0 <= ZnearestE qr <= 11529215046068469760)%Z as round_RANGE.
-  { split.
-    { apply Znearest_IZR_le.
-      lra.
-    }
-    apply Znearest_le_IZR.
-    lra.
-  }
-  rewrite Zle_bool_true in C2 by lia.
-  rewrite Zle_bool_true in C2 by lia.
-  change Int64.max_unsigned with 18446744073709551615%Z.
-  cbn in C2.
-  rewrite C2.
-  cbn.
-  econstructor. split. reflexivity.
-  rewrite Int64.unsigned_repr; cycle 1.
-  { change Int64.max_unsigned with 18446744073709551615%Z.
-    lia.
-  }
-  split.
-  2: assumption.
-  unfold rough_approx_div_thresh, rough_approx_div_longu_smallb_delta in *.
-  apply le_IZR.
-  rewrite <- Rabs_Zabs.
-  rewrite minus_IZR.
-  rewrite mult_IZR.
-  replace  (IZR a' - IZR b' * IZR (ZnearestE qr))%R with
-    (- IZR b' * (IZR (ZnearestE qr) - qr)
-     - IZR b' * (qr - IZR a' / IZR b'))%R by (field ; lra).
-  pose proof (Znearest_imp2  (fun x : Z => negb (Z.even x)) qr) as ROUND.
-  set (delta1 := (IZR (ZnearestE qr) - qr)%R) in *.
-  set (delta2 := (qr - IZR a' / IZR b')%R) in *.
-  unfold smallb_threshold in *.
-