Fappli_Axpy.v 6.88 KB
Newer Older
BOLDO Sylvie's avatar
BOLDO Sylvie committed
1
(**
2 3 4
This file is part of the Flocq formalization of floating-point
arithmetic in Coq: http://flocq.gforge.inria.fr/

5
Copyright (C) 2010-2013 Sylvie Boldo
BOLDO Sylvie's avatar
BOLDO Sylvie committed
6
#<br />#
7
Copyright (C) 2010-2013 Guillaume Melquiond
8 9 10 11 12 13 14 15 16 17 18 19

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.
*)

BOLDO Sylvie's avatar
BOLDO Sylvie committed
20
(** * Necessary conditions to compute a*x+y faithfully *)
BOLDO Sylvie's avatar
BOLDO Sylvie committed
21 22 23 24 25 26 27 28 29 30
Require Import Fcore.

Section Axpy.

Variable beta : radix.
Notation bpow e := (bpow beta e).

Variable prec : Z.
Variable Hp : Zlt 0 prec.

BOLDO Sylvie's avatar
BOLDO Sylvie committed
31
(* FLX first - probably correct in FLTand maybe FTZ? *)
BOLDO Sylvie's avatar
BOLDO Sylvie committed
32 33

Notation format := (generic_format beta (FLX_exp prec)).
34
Notation cexp := (canonic_exp beta (FLX_exp prec)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
35 36
Notation ulp := (ulp beta (FLX_exp prec)).

BOLDO Sylvie's avatar
BOLDO Sylvie committed
37 38 39
Theorem pred_gt_0: forall f,
  format f -> (0 < f)%R -> (0 < pred beta (FLX_exp prec) f)%R.
intros f Hf Zf.
40
unfold pred, Fcore_ulp.ulp, canonic_exp, FLX_exp.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
41 42 43 44 45 46 47 48 49 50 51
destruct (ln_beta beta f) as (ef,Hef).
simpl.
assert (Zf2: (f <>0)%R) by now apply Rgt_not_eq.
specialize (Hef Zf2).
rewrite Rabs_right in Hef.
2: now apply Rgt_ge.
case (Zle_lt_or_eq 1 prec).
omega.
intros Hp1.
case Req_bool_spec; intros H; apply Rlt_Rminus;
    apply Rlt_le_trans with (2:=proj1 Hef);
52
    apply bpow_lt; omega.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
53 54 55 56
(* special case for p = 1 *)
intros Hp1.
case Req_bool_spec; intros H; apply Rlt_Rminus.
apply Rlt_le_trans with (2:=proj1 Hef).
57
apply bpow_lt; omega.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
58 59 60 61 62 63 64
rewrite <- Hp1.
case (proj1 Hef); trivial.
intros H'.
now contradict H.
Qed.


Guillaume Melquiond's avatar
Guillaume Melquiond committed
65 66
Definition MinOrMax x f :=
   ((f = round beta (FLX_exp prec) Zfloor x)
67
     \/ (f = round beta (FLX_exp prec) Zceil x)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
68 69 70 71 72 73

Theorem MinOrMax_opp: forall x f,
   MinOrMax x f <->  MinOrMax (-x) (-f).
assert (forall x f, MinOrMax x f -> MinOrMax (- x) (- f)).
unfold MinOrMax; intros x f [H|H].
right.
74
now rewrite H, round_UP_opp.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
75
left.
76
now rewrite H, round_DN_opp.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
77 78 79 80 81 82 83
intros x f; split; intros H1.
now apply H.
rewrite <- (Ropp_involutive x), <- (Ropp_involutive f).
now apply H.
Qed.


84
Theorem implies_DN_lt_ulp:
BOLDO Sylvie's avatar
BOLDO Sylvie committed
85
  forall x f, format f ->
BOLDO Sylvie's avatar
BOLDO Sylvie committed
86
  (0 < f <= x)%R ->
Guillaume Melquiond's avatar
Guillaume Melquiond committed
87
  (Rabs (f-x) < ulp f)%R ->
88
  (f = round beta (FLX_exp prec) Zfloor x)%R.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
89
intros x f Hf Hxf1 Hxf2.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
90 91
apply sym_eq.
replace x with (f+-(f-x))%R by ring.
92
apply round_DN_succ; trivial.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
93 94 95 96 97 98
apply Hxf1.
replace (- (f - x))%R with (Rabs (f-x)).
split; trivial; apply Rabs_pos.
rewrite Rabs_left1; trivial.
now apply Rle_minus.
Qed.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
99

BOLDO Sylvie's avatar
BOLDO Sylvie committed
100
Theorem MinOrMax_ulp_pred:
BOLDO Sylvie's avatar
BOLDO Sylvie committed
101
  forall x f, format f ->
102
  (0 < f)%R ->
Guillaume Melquiond's avatar
Guillaume Melquiond committed
103
  (Rabs (f-x) < ulp (pred beta (FLX_exp prec) f))%R ->
104
  MinOrMax x f.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
105
intros x f Ff Zf Hf.
106 107 108 109 110 111 112 113 114 115
case (Rlt_or_le x f); intros Hxf2.
(* *)
(* assert (ulp (f-ulp f) = ulp f)%R.

admit.

right; apply sym_eq.
replace f with ((f-ulp f) + (ulp (f-ulp f)))%R.
2: rewrite H; ring.
replace x with ((f-ulp f)+-(f-ulp f-x))%R by ring.
116
apply round_UP_succ; trivial.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
117

118 119 120 121 122 123
apply Hxf1.
replace (- (f - x))%R with (Rabs (f-x)).
split; trivial; apply Rabs_pos.
rewrite Rabs_left1; trivial.
now apply Rle_minus.
*)
BOLDO Sylvie's avatar
BOLDO Sylvie committed
124

125 126 127
admit.
(* *)
left.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
128 129 130 131 132 133 134 135 136 137
apply  implies_DN_lt_ulp; trivial.
now split.
apply Rlt_le_trans with (1:=Hf).
apply ulp_monotone.
clear; intros m n H.
unfold FLX_exp.
omega.
now apply pred_gt_0.
left.
apply pred_lt.
138 139 140 141
Qed.


Theorem implies_MinOrMax_bpow:
BOLDO Sylvie's avatar
BOLDO Sylvie committed
142
  forall x f, format f ->
143
  f = bpow (ln_beta beta f) ->
Guillaume Melquiond's avatar
Guillaume Melquiond committed
144
  (Rabs (f-x) < /2* ulp f)%R ->
BOLDO Sylvie's avatar
BOLDO Sylvie committed
145
  MinOrMax x f.
146 147 148
intros x f Hf1 Hf2 Hxf.


BOLDO Sylvie's avatar
BOLDO Sylvie committed
149 150 151 152 153
Admitted.




BOLDO Sylvie's avatar
BOLDO Sylvie committed
154

155
Variable choice : Z -> bool.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
156 157 158 159 160 161 162

Variable a1 x1 y1 a x y:R.

Hypothesis Ha: format a.
Hypothesis Hx: format x.
Hypothesis Hy: format y.

163 164
Notation t := (round beta (FLX_exp prec) (Znearest choice) (a*x)).
Notation u := (round beta (FLX_exp prec) (Znearest choice) (t+y)).
BOLDO Sylvie's avatar
BOLDO Sylvie committed
165

BOLDO Sylvie's avatar
BOLDO Sylvie committed
166 167
(*
Axpy_aux1 : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 < u
Guillaume Melquiond's avatar
Guillaume Melquiond committed
168
   => abs(t-a*x)          <=  Fulp(b)(Fpred(b)(u))/4
BOLDO Sylvie's avatar
BOLDO Sylvie committed
169 170 171 172 173 174
   => abs(y1-y+a1*x1-a*x) <   Fulp(b)(Fpred(b)(u))/4
   => MinOrMax?(y1+a1*x1,u)


Axpy_aux1_aux1 : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 < u
  => Fnormal?(b)(t) => radix*abs(t) <= Fpred(b)(u)
Guillaume Melquiond's avatar
Guillaume Melquiond committed
175
  => abs(t-a*x)  <=  Fulp(b)(Fpred(b)(u))/4
BOLDO Sylvie's avatar
BOLDO Sylvie committed
176 177 178

Axpy_aux1_aux2 : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 < u
  => Fsubnormal?(b)(t) => 1-dExp(b) <= Fexp(Fpred(b)(u))
Guillaume Melquiond's avatar
Guillaume Melquiond committed
179
  => abs(t-a*x) <=  Fulp(b)(Fpred(b)(u))/4
BOLDO Sylvie's avatar
BOLDO Sylvie committed
180 181 182 183 184 185 186 187

Axpy_aux2 : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 < u
  => Fsubnormal?(b)(t) => u=t+y
  => abs(y1-y+a1*x1-a*x) < Fulp(b)(Fpred(b)(u))/4
  => MinOrMax?(y1+a1*x1,u)


Axpy_aux3 : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 < u
Guillaume Melquiond's avatar
Guillaume Melquiond committed
188
  => Fsubnormal?(b)(t)
BOLDO Sylvie's avatar
BOLDO Sylvie committed
189 190 191 192 193 194 195 196 197 198
  => -dExp(b) = Fexp(Fpred(b)(u)) =>  1-dExp(b) <= Fexp(u)
  => abs(y1-y+a1*x1-a*x) < Fulp(b)(Fpred(b)(u))/4
  => MinOrMax?(y1+a1*x1,u)


AxpyPos : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 < u
  => (Fnormal?(b)(t) => radix*abs(t) <= Fpred(b)(u))
  => abs(y1-y+a1*x1-a*x) < Fulp(b)(Fpred(b)(u))/4
  => MinOrMax?(y1+a1*x1,u)

Guillaume Melquiond's avatar
Guillaume Melquiond committed
199 200
Axpy_opt_aux1_aux1 : lemma Fnormal?(b)(t) => Fnormal?(b)(u) => 0 < u
    => Prec(b) >= 3 =>
BOLDO Sylvie's avatar
BOLDO Sylvie committed
201 202 203 204
   (1+radix*(1+1/(2*abs(Fnum(u))))*(1+1/abs(Fnum(Fpred(b)(u)))))/(1-1/(2*abs(Fnum(t))))
      <= 1+radix+radix^(4-Prec(b))

Axpy_opt_aux1 : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 < u
Guillaume Melquiond's avatar
Guillaume Melquiond committed
205
  =>  Prec(b) >= 3 => Fnormal?(b)(t)
BOLDO Sylvie's avatar
BOLDO Sylvie committed
206 207 208 209
  => (radix+1+radix^(4-Prec(b)))*abs(a*x) <= abs(y)
  =>  radix*abs(t) <= Fpred(b)(u)

Axpy_opt_aux2 : lemma  Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 < u
Guillaume Melquiond's avatar
Guillaume Melquiond committed
210
  =>  Prec(b) >= 6 =>  Fnormal?(b)(t)
BOLDO Sylvie's avatar
BOLDO Sylvie committed
211 212 213 214
  => (radix+1+radix^(4-Prec(b)))*abs(a*x) <= abs(y)
  =>  abs(y)*radix^(1-Prec(b))/(radix+1) < Fulp(b)(Fpred(b)(u))

Axpy_opt_aux3 : lemma  Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 < u
Guillaume Melquiond's avatar
Guillaume Melquiond committed
215
  =>  Prec(b) >= 3  =>  Fsubnormal?(b)(t)
BOLDO Sylvie's avatar
BOLDO Sylvie committed
216 217 218 219 220
  => (radix+1+radix^(4-Prec(b)))*abs(a*x) <= abs(y)
  =>  abs(y)*radix^(1-Prec(b))/(radix+radix/2) <= Fulp(b)(Fpred(b)(u))


Axpy_optPos : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 < u
Guillaume Melquiond's avatar
Guillaume Melquiond committed
221
  => Prec(b) >= 6
BOLDO Sylvie's avatar
BOLDO Sylvie committed
222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237
  => (radix+1+radix^(4-Prec(b)))*abs(a*x) <= abs(y)
  => abs(y1-y+a1*x1-a*x) < abs(y)*radix^(1-Prec(b))/(6*radix)
  => MinOrMax?(y1+a1*x1,u)

Axpy_optZero : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u) => 0 = u
  => (radix+1+radix^(4-Prec(b)))*abs(a*x) <= abs(y)
  => abs(y1-y+a1*x1-a*x) < abs(y)*radix^(1-Prec(b))/(6*radix)
  => MinOrMax?(y1+a1*x1,u)

Axpy_opt : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u)
  => Prec(b) >= 6
  => (radix+1+radix^(4-Prec(b)))*abs(a*x) <= abs(y)
  => abs(y1-y+a1*x1-a*x) < abs(y)*radix^(1-Prec(b))/(6*radix)
  => MinOrMax?(y1+a1*x1,u)

Axpy_simpl : lemma Closest?(b)(a*x,t) => Closest?(b)(t+y,u)
Guillaume Melquiond's avatar
Guillaume Melquiond committed
238
  => Prec(b) >= 24 => radix=2
BOLDO Sylvie's avatar
BOLDO Sylvie committed
239 240 241 242
  => (3+1/100000)*abs(a*x) <= abs(y)
  => abs(y1-y+a1*x1-a*x) < abs(y)*2^(1-Prec(b))/12
  => MinOrMax?(y1+a1*x1,u)
*)
BOLDO Sylvie's avatar
BOLDO Sylvie committed
243 244 245 246 247 248 249 250

Theorem Axpy_opt :
  (6 <= prec)%Z ->
 ((bpow 1 +1 + bpow (4-prec))* Rabs (a*x) <= Rabs y)%R ->
   (Rabs (y1 - y + a1 * x1 - a * x) <=
      bpow (1-prec) / (6*bpow 1) * Rabs y)%R ->
         (MinOrMax (a1 * x1 + y1) u).
Admitted.
BOLDO Sylvie's avatar
BOLDO Sylvie committed
251 252 253


End Axpy.