From 0fb548cbf78c19430e60827f7c253fc42aa25e05 Mon Sep 17 00:00:00 2001 From: David Monniaux Date: Tue, 25 Jan 2022 10:54:19 +0100 Subject: some progress? --- kvx/FPDivision64.v | 311 +++++++++++++++++++++++++++++++---------------------- 1 file changed, 185 insertions(+), 126 deletions(-) (limited to 'kvx') 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. -- cgit