diff --git a/src/Calc/Fcalc_div.v b/src/Calc/Fcalc_div.v new file mode 100644 index 0000000000000000000000000000000000000000..1645ee25cfbeed4d07568fa89fe6744d9bce8fe8 --- /dev/null +++ b/src/Calc/Fcalc_div.v @@ -0,0 +1,139 @@ +Require Import Fcore. +Require Import Fcalc_bracket. +Require Import Fcalc_digits. + +Section Fcalc_div. + +Variable beta : radix. +Notation bpow e := (bpow beta e). + +Variable fexp : Z -> Z. +Hypothesis prop_exp : valid_exp fexp. +Notation format := (generic_format beta fexp). + +(* + * 1. Shift dividend mantissa so that it has at least p2 + p digits. + * 2. Perform the euclidean division. + * 3. Compute position with remainder. + * 4. Round to p digits. + * + * Complexity is fine as long as p1 <= 2p and p2 <= p. + *) +Definition Fdiv_aux prec m1 e1 m2 e2 := + let d1 := digits beta m1 in + let d2 := digits beta m2 in + let e := (e1 - e2)%Z in + let (m, e') := + match (d2 + prec - d1)%Z with + | Zpos p => (m1 * Zpower_pos (radix_val beta) p, e + Zneg p)%Z + | _ => (m1, e) + end in + let '(q, r) := Zdiv_eucl m m2 in + (q, e', new_location m2 r loc_Exact). + +Theorem Fdiv_aux_correct : + forall prec m1 e1 m2 e2, + (0 < prec)%Z -> + (0 < m1)%Z -> (1 < m2)%Z -> + let '(m, e, l) := Fdiv_aux prec m1 e1 m2 e2 in + (prec <= digits beta m)%Z /\ + inbetween_float beta m e (F2R (Float beta m1 e1) / F2R (Float beta m2 e2)) l. +Proof. +intros prec m1 e1 m2 e2 Hprec Hm1 Hm2. +unfold Fdiv_aux. +set (d1 := digits beta m1). +set (d2 := digits beta m2). +case_eq + (match (d2 + prec - d1)%Z with + | Zpos p => ((m1 * Zpower_pos (radix_val beta) p)%Z, (e1 - e2 + Zneg p)%Z) + | _ => (m1, (e1 - e2)%Z) + end). +intros m' e' Hme. +(* . the shifted mantissa m' has enough digits *) +assert (Hs: F2R (Float beta m' (e' + e2)) = F2R (Float beta m1 e1) /\ (0 < m')%Z /\ (d2 + prec <= digits beta m')%Z). +replace (d2 + prec)%Z with (d2 + prec - d1 + d1)%Z by ring. +destruct (d2 + prec - d1)%Z as [|p|p] ; + unfold d1 ; + inversion Hme. +ring_simplify (e1 - e2 + e2)%Z. +repeat split. +now rewrite <- H0. +apply Zle_refl. +replace (e1 - e2 + Zneg p + e2)%Z with (e1 - Zpos p)%Z by (unfold Zminus ; simpl ; ring). +fold (Zpower (radix_val beta) (Zpos p)). +split. +pattern (Zpos p) at 1 ; replace (Zpos p) with (e1 - (e1 - Zpos p))%Z by ring. +apply sym_eq. +apply F2R_change_exp. +assert (0 < Zpos p)%Z by easy. +omega. +split. +apply Zmult_lt_0_compat. +exact Hm1. +apply Zpower_gt_0. +generalize (radix_prop beta). +omega. +easy. +rewrite digits_shift. +rewrite Zplus_comm. +apply Zle_refl. +apply sym_not_eq. +now apply Zlt_not_eq. +easy. +split. +now ring_simplify (e1 - e2 + e2)%Z. +assert (Zneg p < 0)%Z by easy. +omega. +(* . *) +destruct Hs as (Hs1, (Hs2, Hs3)). +rewrite <- Hs1. +assert (Hm2': (0 < m2)%Z) by omega. +generalize (Z_div_mod m' m2 (Zlt_gt _ _ Hm2')). +destruct (Zdiv_eucl m' m2) as (q, r). +intros (Hq, Hr). +split. +(* . the result mantissa q has enough digits *) +cut (digits beta m' <= d2 + digits beta q)%Z. omega. +unfold d2. +rewrite Hq. +assert (Hq': (0 < q)%Z). +apply Zmult_lt_reg_r with m2. +exact Hm2'. +assert (m2 < m')%Z. +apply digits_lt with beta. +exact Hs2. +unfold d2 in Hs3. +omega. +cut (q * m2 = m' - r)%Z. omega. +rewrite Hq. +ring. +apply Zle_trans with (digits beta (m2 + q + m2 * q)). +apply digits_le. +now rewrite <- Hq. +omega. +apply digits_mult_strong. +omega. +exact Hq'. +(* . the location is correctly computed *) +unfold inbetween_float, F2R. simpl. +rewrite bpow_add, plus_Z2R. +rewrite Hq, plus_Z2R, mult_Z2R. +replace ((Z2R m2 * Z2R q + Z2R r) * (bpow e' * bpow e2) / (Z2R m2 * bpow e2))%R + with ((Z2R q + Z2R r / Z2R m2) * bpow e')%R. +apply inbetween_mult_compat. +apply bpow_gt_0. +replace (Z2R 1) with (Z2R m2 * /Z2R m2)%R. +apply new_location_correct ; try easy. +apply Rinv_0_lt_compat. +now apply (Z2R_lt 0). +now constructor. +apply Rinv_r. +apply Rgt_not_eq. +now apply (Z2R_lt 0). +field. +split ; apply Rgt_not_eq. +apply bpow_gt_0. +now apply (Z2R_lt 0). +Qed. + +End Fcalc_div. \ No newline at end of file diff --git a/src/Makefile.am b/src/Makefile.am index 94a5af608519ded7878dbc0460b0ddc18c810262..2046256d58e35d1c057c4098385d6bacf5c377d9 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -13,6 +13,7 @@ FILES = \ Core/Fcore.v \ Calc/Fcalc_bracket.v \ Calc/Fcalc_digits.v \ + Calc/Fcalc_div.v \ Calc/Fcalc_ops.v \ Prop/Fprop_nearest.v