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
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
|
From Flocq Require Import Core Digits Operations Round Bracket Sterbenz
Binary Round_odd Bits.
Require Archi.
Require Import Coqlib.
Require Import Compopts.
Require Import AST.
Require Import Integers.
Require Import Floats.
Require Import Op.
Require Import CminorSel.
Require Import OpHelpers.
Require Import ExtFloats.
Require Import DecBoolOps.
Require Import Chunks.
Require Import Builtins.
Require Import Values Globalenvs.
Require Compopts.
Require Import Psatz.
Require Import IEEE754_extra.
From Gappa Require Import Gappa_tactic.
Local Open Scope cminorsel_scope.
Local Open Scope string_scope.
Local Open Scope error_monad_scope.
Definition approx_inv b :=
let invb_s := Eop Oinvfs ((Eop Osingleofintu ((Eletvar 0%nat):::Enil)):::Enil) in
let invb_d := Eop Ofloatofsingle (invb_s ::: Enil) in
let b_d := Eop Ofloatoflongu ((Eop Ocast32unsigned ((Eletvar 1%nat):::Enil)):::Enil) in
let invb_d_var := Eletvar (0%nat) in
let one := Eop (Ofloatconst ExtFloat.one) Enil in
let alpha := Eop Ofmsubf (one ::: invb_d_var ::: b_d ::: Enil) in
let x := Eop Ofmaddf (invb_d_var ::: alpha ::: invb_d_var ::: Enil) in
Elet b (Elet invb_d x).
Definition approx_inv_thresh := (0.000000000000001)%R.
(*
Lemma BofZ_exact_zero:
forall (prec emax : Z) (prec_gt_0_ : Prec_gt_0 prec)
(Hmax : (prec < emax)%Z),
B2R (BofZ prec emax 0%Z (Hmax := Hmax)) = 0%R /\
is_finite (BofZ prec emax 0%Z (Hmax := Hmax)) = true /\
Bsign prec emax (BofZ prec emax 0%Z (Hmax := Hmax)) = false.
Proof.
intros.
apply BofZ_exact.
pose proof (Z.pow_nonneg 2 prec).
lia.
Qed.
*)
Theorem approx_inv_correct :
forall (ge : genv) (sp: val) cmenv memenv
(le : letenv) (expr_b : expr) (b : int)
(EVAL_b : eval_expr ge sp cmenv memenv le expr_b (Vint b))
(b_nz : ((Int.unsigned b) > 0)%Z),
exists f : float,
eval_expr ge sp cmenv memenv le (approx_inv expr_b) (Vfloat f) /\
is_finite _ _ f = true /\ (Rabs((B2R _ _ f) - (1 / IZR (Int.unsigned b))) <= approx_inv_thresh)%R.
Proof.
intros. unfold approx_inv.
econstructor. constructor.
{ repeat econstructor.
{ eassumption. }
{ reflexivity. } }
set (invb_d := (Float.of_single (ExtFloat32.inv (Float32.of_intu b)))).
set (b_d := (Float.of_longu (Int64.repr (Int.unsigned b)))).
cbn.
pose proof (Int.unsigned_range b) as RANGE.
change Int.modulus with 4294967296%Z in RANGE.
set (b' := Int.unsigned b) in *.
assert (0 <= IZR b' <= 4294967295)%R as RANGE'.
{ split.
{ apply IZR_le. lia. }
apply IZR_le. lia.
}
pose proof (BofZ_correct 24 128 eq_refl eq_refl b') as C1.
rewrite Rlt_bool_true in C1; cycle 1.
{ clear C1. cbn.
admit.
}
rewrite Zlt_bool_false in C1 by lia.
destruct C1 as (C1E & C1F & C1S).
Local Transparent Float32.of_intu Float32.of_int Float32.div.
unfold ExtFloat32.inv, ExtFloat32.one, Float32.of_intu, Float32.of_int, Float32.div in invb_d.
fold b' in invb_d.
Check BofZ.
change (Int.signed (Int.repr 1%Z)) with 1%Z in invb_d.
pose proof (Bdiv_correct 24 128 eq_refl eq_refl Float32.binop_nan mode_NE
(BofZ 24 128 eq_refl eq_refl 1)
(BofZ 24 128 eq_refl eq_refl b')) as C2.
rewrite Rlt_bool_true in C2; cycle 1.
{ clear C2. admit. }
assert (B2R 24 128 (BofZ 24 128 eq_refl eq_refl b') <> 0%R) as NONZ.
{ admit. }
destruct (C2 NONZ) as (C2E & C2F & C2S).
clear C2 NONZ.
Local Transparent Float.of_single.
unfold Float.of_single in invb_d.
pose proof (Bconv_correct 24 128 53 1024 (@eq_refl Datatypes.comparison Lt)
(@eq_refl Datatypes.comparison Lt) Float.of_single_nan mode_NE
(Bdiv 24 128 (@eq_refl Datatypes.comparison Lt)
(@eq_refl Datatypes.comparison Lt) Float32.binop_nan mode_NE
(BofZ 24 128 (@eq_refl Datatypes.comparison Lt)
(@eq_refl Datatypes.comparison Lt) 1)
(BofZ 24 128 (@eq_refl Datatypes.comparison Lt)
(@eq_refl Datatypes.comparison Lt) b'))) as C3.
fold invb_d in C3.
rewrite Rlt_bool_true in C3; cycle 1.
{ clear C3. admit. }
change (is_finite 24 128 (BofZ 24 128 eq_refl eq_refl 1)) with true in C2F.
destruct (C3 C2F) as (C3E & C3F & C3S).
unfold Float.fma.
assert (is_finite _ _ (Float.neg invb_d) = true) as invb_d_F.
{ Local Transparent Float.neg.
unfold Float.neg.
rewrite is_finite_Bopp.
assumption.
}
Set Printing Implicit.
pose proof (Bfma_correct 53 1024 (@eq_refl Datatypes.comparison Lt)
(@eq_refl Datatypes.comparison Lt) Float.fma_nan mode_NE
(Float.neg invb_d) b_d ExtFloat.one) as C4.
Admitted.
Definition approx_div a b :=
let a_var := Eletvar (1%nat) in
let b_var := Eletvar (0%nat) in
let a_d := Eop Ofloatoflongu ((Eop Ocast32unsigned (a_var ::: Enil)) ::: Enil) in
let qr := Eop Olonguoffloat_ne ((Eop Omulf (a_d:::(approx_inv b_var):::Enil)):::Enil) in
let qr_var := Eletvar 0%nat in
let rem := Eop Omsubl ((Eop Ocast32unsigned ((Eletvar (2%nat)):::Enil)):::
qr_var :::
(Eop Ocast32unsigned ((Eletvar (1%nat)):::Enil)):::Enil) in
let qr_m1 := Eop (Oaddlimm (Int64.repr (-1)%Z)) (qr_var:::Enil) in
let cases := Eop (Osel (Ccompl0 Clt) Tlong)
(qr_m1 ::: qr_var ::: rem ::: Enil) in (* (Elet qr cases) *)
Eop Olowlong ((Elet a (Elet (lift b) (Elet qr cases))) ::: Enil).
Open Scope Z.
Definition div_approx_reals (a b : Z) (x : R) :=
let q:=ZnearestE x in
let r:=a-q*b in
if r <? 0
then q-1
else q.
Lemma floor_ball1:
forall x : R, forall y : Z,
(Rabs (x - IZR y) < 1)%R -> Zfloor x = (y-1)%Z \/ Zfloor x = y.
Proof.
intros x y BALL.
apply Rabs_lt_inv in BALL.
case (Rcompare_spec x (IZR y)); intro CMP.
- left. apply Zfloor_imp.
ring_simplify (y-1+1).
rewrite minus_IZR. lra.
- subst.
rewrite Zfloor_IZR. right. reflexivity.
- right. apply Zfloor_imp.
rewrite plus_IZR. lra.
Qed.
Theorem div_approx_reals_correct:
forall a b : Z, forall x : R,
b > 0 ->
(Rabs (x - IZR a/ IZR b) < 1/2)%R ->
div_approx_reals a b x = (a/b)%Z.
Proof.
intros a b x bPOS GAP.
assert (0 < IZR b)%R by (apply IZR_lt ; lia).
unfold div_approx_reals.
pose proof (Znearest_imp2 (fun x => negb (Z.even x)) x) as NEAR.
assert (Rabs (IZR (ZnearestE x) - IZR a/ IZR b) < 1)%R as BALL.
{ pose proof (Rabs_triang (IZR (ZnearestE x) - x)
(x - IZR a/ IZR b)) as TRI.
ring_simplify (IZR (ZnearestE x) - x + (x - IZR a / IZR b))%R in TRI.
lra.
}
clear GAP NEAR.
rewrite Rabs_minus_sym in BALL.
pose proof (floor_ball1 _ _ BALL) as FLOOR.
clear BALL.
rewrite Zfloor_div in FLOOR by lia.
pose proof (Z_div_mod_eq_full a b) as DIV_MOD.
assert (0 < b) as bPOS' by lia.
pose proof (Z.mod_pos_bound a b bPOS') as MOD_BOUNDS.
case Z.ltb_spec; intro; destruct FLOOR; lia.
Qed.
Opaque approx_inv.
Theorem approx_div_correct :
forall (ge : genv) (sp: val) cmenv memenv
(le : letenv) (expr_a expr_b : expr) (a b : int)
(EVAL_a : eval_expr ge sp cmenv memenv le expr_a (Vint a))
(EVAL_b : eval_expr ge sp cmenv memenv le expr_b (Vint b))
(b_nz : (Int.unsigned b > 0)%Z),
eval_expr ge sp cmenv memenv le (approx_div expr_a expr_b)
(Vint (Int.divu a b)).
Proof.
intros.
assert (eval_expr ge sp cmenv memenv (Vint b :: Vint a :: le)
(Eletvar 0) (Vint b)) as EVAL_b'.
{ constructor. reflexivity. }
destruct (approx_inv_correct ge sp cmenv memenv (Vint b :: Vint a :: le) (Eletvar 0) b EVAL_b' b_nz) as (inv_b & inv_b_eval & inv_b_finite & inv_b_correct).
unfold approx_div.
repeat econstructor.
{ exact EVAL_a. }
{ apply eval_lift. exact EVAL_b. }
exact inv_b_eval.
cbn. f_equal.
rewrite <- Float.of_intu_of_longu.
unfold Float.to_longu_ne.
rewrite ZofB_ne_range_correct by lia.
set (prod := (Float.mul (Float.of_intu a) inv_b)).
set (a' := Int.unsigned a) in *.
set (b' := Int.unsigned b) in *.
set (prod' := (B2R _ _ prod)).
set (prod_r := ZnearestE prod') in *.
replace (_ && _ && _) with true by admit.
cbn.
f_equal.
unfold Int.divu.
rewrite <- div_approx_reals_correct with (x := B2R _ _ prod).
2: lia.
{
unfold div_approx_reals.
fold a' b' prod' prod_r.
unfold Int64.mul.
rewrite Int64.unsigned_repr by admit.
rewrite Int64.unsigned_repr by admit.
unfold Int64.sub.
rewrite Int64.unsigned_repr by admit.
rewrite Int64.unsigned_repr by admit.
case Z.ltb_spec; intro CMP.
- replace (Int64.lt (Int64.repr (a' - prod_r * b')) Int64.zero) with true by admit.
cbn.
f_equal.
admit.
- replace (Int64.lt (Int64.repr (a' - prod_r * b')) Int64.zero) with false by admit.
cbn.
f_equal.
admit.
}
fold prod'.
admit.
Admitted.
|
