1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
|
(**
This file is part of the Flocq formalization of floating-point
arithmetic in Coq: http://flocq.gforge.inria.fr/
Copyright (C) 2010-2018 Sylvie Boldo
#<br />#
Copyright (C) 2010-2018 Guillaume Melquiond
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 3 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
COPYING file for more details.
*)
(** * Helper function and theorem for computing the rounded quotient of two floating-point numbers. *)
From Coq Require Import Lia.
Require Import Raux Defs Generic_fmt Float_prop Digits Bracket.
Set Implicit Arguments.
Set Strongly Strict Implicit.
Section Fcalc_div.
Variable beta : radix.
Notation bpow e := (bpow beta e).
Variable fexp : Z -> Z.
(** Computes a mantissa of precision p, the corresponding exponent,
and the position with respect to the real quotient of the
input floating-point numbers.
The algorithm performs the following steps:
- Shift dividend mantissa so that it has at least p2 + p digits.
- Perform the Euclidean division.
- Compute the position according to the division remainder.
Complexity is fine as long as p1 <= 2p and p2 <= p.
*)
Lemma mag_div_F2R :
forall m1 e1 m2 e2,
(0 < m1)%Z -> (0 < m2)%Z ->
let e := ((Zdigits beta m1 + e1) - (Zdigits beta m2 + e2))%Z in
(e <= mag beta (F2R (Float beta m1 e1) / F2R (Float beta m2 e2)) <= e + 1)%Z.
Proof.
intros m1 e1 m2 e2 Hm1 Hm2.
rewrite <- (mag_F2R_Zdigits beta m1 e1) by now apply Zgt_not_eq.
rewrite <- (mag_F2R_Zdigits beta m2 e2) by now apply Zgt_not_eq.
apply mag_div.
now apply F2R_neq_0, Zgt_not_eq.
now apply F2R_neq_0, Zgt_not_eq.
Qed.
Definition Fdiv_core m1 e1 m2 e2 e :=
let (m1', m2') :=
if Zle_bool e (e1 - e2)%Z
then (m1 * Zpower beta (e1 - e2 - e), m2)%Z
else (m1, m2 * Zpower beta (e - (e1 - e2)))%Z in
let '(q, r) := Z.div_eucl m1' m2' in
(q, new_location m2' r loc_Exact).
Theorem Fdiv_core_correct :
forall m1 e1 m2 e2 e,
(0 < m1)%Z -> (0 < m2)%Z ->
let '(m, l) := Fdiv_core m1 e1 m2 e2 e in
inbetween_float beta m e (F2R (Float beta m1 e1) / F2R (Float beta m2 e2)) l.
Proof.
intros m1 e1 m2 e2 e Hm1 Hm2.
unfold Fdiv_core.
match goal with |- context [if ?b then ?b1 else ?b2] => set (m12 := if b then b1 else b2) end.
case_eq m12 ; intros m1' m2' Hm.
assert ((F2R (Float beta m1 e1) / F2R (Float beta m2 e2) = IZR m1' / IZR m2' * bpow e)%R /\ (0 < m2')%Z) as [Hf Hm2'].
{ unfold F2R, Zminus ; simpl.
destruct (Zle_bool e (e1 - e2)) eqn:He' ; injection Hm ; intros ; subst.
- split ; try easy.
apply Zle_bool_imp_le in He'.
rewrite mult_IZR, IZR_Zpower by lia.
unfold Zminus ; rewrite 2!bpow_plus, 2!bpow_opp.
field.
repeat split ; try apply Rgt_not_eq, bpow_gt_0.
now apply IZR_neq, Zgt_not_eq.
- apply Z.leb_gt in He'.
split ; cycle 1.
{ apply Z.mul_pos_pos with (1 := Hm2).
apply Zpower_gt_0 ; lia. }
rewrite mult_IZR, IZR_Zpower by lia.
unfold Zminus ; rewrite bpow_plus, bpow_opp, bpow_plus, bpow_opp.
field.
repeat split ; try apply Rgt_not_eq, bpow_gt_0.
now apply IZR_neq, Zgt_not_eq. }
clearbody m12 ; clear Hm Hm1 Hm2.
generalize (Z_div_mod m1' m2' (Z.lt_gt _ _ Hm2')).
destruct (Z.div_eucl m1' m2') as (q, r).
intros (Hq, Hr).
rewrite Hf.
unfold inbetween_float, F2R. simpl.
rewrite Hq, 2!plus_IZR, mult_IZR.
apply inbetween_mult_compat.
apply bpow_gt_0.
destruct (Z_lt_le_dec 1 m2') as [Hm2''|Hm2''].
- replace 1%R with (IZR m2' * /IZR m2')%R.
apply new_location_correct ; try easy.
apply Rinv_0_lt_compat.
now apply IZR_lt.
constructor.
field.
now apply IZR_neq, Zgt_not_eq.
field.
now apply IZR_neq, Zgt_not_eq.
- assert (r = 0 /\ m2' = 1)%Z as [-> ->] by (clear -Hr Hm2'' ; lia).
unfold Rdiv.
rewrite Rmult_1_l, Rplus_0_r, Rinv_1, Rmult_1_r.
now constructor.
Qed.
Definition Fdiv (x y : float beta) :=
let (m1, e1) := x in
let (m2, e2) := y in
let e' := ((Zdigits beta m1 + e1) - (Zdigits beta m2 + e2))%Z in
let e := Z.min (Z.min (fexp e') (fexp (e' + 1))) (e1 - e2) in
let '(m, l) := Fdiv_core m1 e1 m2 e2 e in
(m, e, l).
Theorem Fdiv_correct :
forall x y,
(0 < F2R x)%R -> (0 < F2R y)%R ->
let '(m, e, l) := Fdiv x y in
(e <= cexp beta fexp (F2R x / F2R y))%Z /\
inbetween_float beta m e (F2R x / F2R y) l.
Proof.
intros [m1 e1] [m2 e2] Hm1 Hm2.
apply gt_0_F2R in Hm1.
apply gt_0_F2R in Hm2.
unfold Fdiv.
generalize (mag_div_F2R m1 e1 m2 e2 Hm1 Hm2).
set (e := Zminus _ _).
set (e' := Z.min (Z.min (fexp e) (fexp (e + 1))) (e1 - e2)).
intros [H1 H2].
generalize (Fdiv_core_correct m1 e1 m2 e2 e' Hm1 Hm2).
destruct Fdiv_core as [m' l].
apply conj.
apply Z.le_trans with (1 := Z.le_min_l _ _).
unfold cexp.
destruct (Zle_lt_or_eq _ _ H1) as [H|H].
- replace (fexp (mag _ _)) with (fexp (e + 1)).
apply Z.le_min_r.
clear -H1 H2 H ; apply f_equal ; lia.
- replace (fexp (mag _ _)) with (fexp e).
apply Z.le_min_l.
clear -H1 H2 H ; apply f_equal ; lia.
Qed.
End Fcalc_div.
|
