stdlib/ufloat.mlw

   1 (** {1 Unbounded floating-point numbers}
   2
   3
   4 See also {h <a href="https://inria.hal.science/hal-04343157">this report</a>}
   5
   6 *)
   7
   8 module UFloat
   9   use real.RealInfix
  10   use real.FromInt
  11   use real.Abs
  12   use ieee_float.RoundingMode
  13
  14   type t
  15
  16   val function uround mode real : t
  17   val function to_real t : real
  18   val function of_int int : t
  19   axiom to_real_of_int : forall x [of_int x]. to_real (of_int x) = from_int x
  20
  21   val function uzero : t
  22   axiom uzero_spec : to_real uzero = 0.0
  23
  24   val function uone : t
  25   axiom uone_spec : to_real uone = 1.0
  26
  27   val function utwo : t
  28   axiom utwo_spec : to_real utwo = 2.0
  29
  30   constant eps:real
  31   constant eta:real
  32   axiom eps_bounds : 0. <. eps <. 1.
  33   axiom eta_bounds : 0. <. eta <. 1.
  34
  35   (* To avoid "inline_trivial" to break the forward_propagation strategy *)
  36   meta "inline:no" function eps
  37   meta "inline:no" function eta
  38
  39   let function uadd (x y:t) : t
  40     (* TODO: Do we want the two first assertions in our context ?
  41       We only use them to prove the addition lemma *)
  42     ensures { abs (to_real result -. (to_real x +. to_real y)) <=. abs (to_real x) }
  43     ensures { abs (to_real result -. (to_real x +. to_real y)) <=. abs (to_real y) }
  44     ensures {
  45       abs (to_real result -. (to_real x +. to_real y))
  46       <=. abs (to_real x +. to_real y) *. eps
  47     }
  48   = uround RNE (to_real x +. to_real y)
  49
  50   let function usub (x y:t) : t
  51     (* TODO: Do we want the two first assertions in our context ?
  52       We only use them to prove the addition lemma *)
  53     ensures { abs (to_real result -. (to_real x -. to_real y)) <=. abs (to_real x) }
  54     ensures { abs (to_real result -. (to_real x -. to_real y)) <=. abs (to_real y) }
  55     ensures {
  56       abs (to_real result -. (to_real x -. to_real y))
  57       <=. abs (to_real x -. to_real y) *. eps
  58     }
  59   = uround RNE (to_real x -. to_real y)
  60
  61   let function umul (x y:t) : t
  62     ensures {
  63       abs (to_real result -. (to_real x *. to_real y))
  64         <=. abs (to_real x *. to_real y) *. eps +. eta
  65     }
  66   = uround RNE (to_real x *. to_real y)
  67
  68   let function udiv (x y:t) : t
  69     requires { to_real y <> 0. }
  70     ensures {
  71       abs (to_real result -. (to_real x /. to_real y))
  72         <=. abs (to_real x /. to_real y) *. eps +. eta
  73     }
  74   = uround RNE (to_real x /. to_real y)
  75
  76   let function uminus (x:t) : t
  77     ensures { to_real result = -. (to_real x) }
  78   = uround RNE (-. (to_real x))
  79
  80   predicate is_exact (uop : t -> t -> t) (x y :t)
  81
  82   (* Exact division but can still underflow, giving eta as error *)
  83   let function udiv_exact (x y:t) : t
  84     requires { to_real y <> 0. }
  85     requires { is_exact udiv x y }
  86     ensures { abs (to_real result -. (to_real x /. to_real y)) <=. eta }
  87   = uround RNE (to_real x /. to_real y)
  88
  89   (** Infix operators *)
  90   let function ( ++. ) (x:t) (y:t) : t = uadd x y
  91   let function ( --. ) (x:t) (y:t) : t = usub x y
  92   let function ( **. ) (x:t) (y:t) : t = umul x y
  93   (* Why3 doesn't support abbreviations so we need to add the requires *)
  94   let function ( //. ) (x:t) (y:t) : t
  95     requires { to_real y <> 0. }
  96   = udiv x y
  97   let function ( --._ ) (x:t) : t = uminus x
  98   let function ( ///. ) (x:t) (y:t) : t
  99     requires { to_real y <> 0. }
 100     requires { is_exact udiv x y }
 101   = udiv_exact x y
 102
 103   (* Some constants *)
 104   constant u0:t
 105   axiom to_real_u0 : to_real u0 = 0.0
 106   constant u1:t
 107   axiom to_real_u1 : to_real u1 = 1.0
 108   constant u2:t
 109   axiom to_real_u2 : to_real u2 = 2.0
 110   constant u4:t
 111   axiom to_real_u4 : to_real u4 = 4.0
 112   constant u8:t
 113   axiom to_real_u8 : to_real u8 = 8.0
 114   constant u16:t
 115   axiom to_real_u16 : to_real u16 = 16.0
 116   constant u32:t
 117   axiom to_real_u32 : to_real u32 = 32.0
 118   constant u64:t
 119   axiom to_real_u64 : to_real u64 = 64.0
 120   constant u128:t
 121   axiom to_real_u128 : to_real u128 = 128.0
 122   constant u256:t
 123   axiom to_real_u256 : to_real u256 = 256.0
 124   constant u512:t
 125   axiom to_real_u512 : to_real u512 = 512.0
 126   constant u1024:t
 127   axiom to_real_u1024 : to_real u1024 = 1024.0
 128   constant u2048:t
 129   axiom to_real_u2048 : to_real u2048 = 2048.0
 130   constant u4096:t
 131   axiom to_real_u4096 : to_real u4096 = 4096.0
 132   constant u8192:t
 133   axiom to_real_u8192 : to_real u8192 = 8192.0
 134   constant u16384:t
 135   axiom to_real_u16384 : to_real u16384 = 16384.0
 136   constant u32768:t
 137   axiom to_real_u32768 : to_real u32768 = 32768.0
 138   constant u65536:t
 139   axiom to_real_u65536 : to_real u65536 = 65536.0
 140
 141   predicate is_positive_power_of_2 (x:t) =
 142     x = u1 \/ x = u2 || x = u4 || x = u8 || x = u16 || x = u32 || x = u64
 143     || x = u128 \/ x = u256 || x = u4096 || x = u8192 || x = u16384 || x = u32768
 144     || x = u65536
 145
 146   axiom div_by_positive_power_of_2_is_exact :
 147     forall x y. is_positive_power_of_2 y -> is_exact udiv x y
 148 end
 149
 150 (** {2 Single-precision unbounded floats} *)
 151 module USingle
 152   use real.RealInfix
 153
 154   type usingle
 155
 156   constant eps:real = 0x1p-24 /. (1. +. 0x1p-24)
 157   constant eta:real = 0x1p-150
 158
 159   clone export UFloat with
 160     type t = usingle,
 161     constant eps = eps,
 162     constant eta = eta,
 163     axiom.
 164 end
 165
 166
 167 (** {3 Double-precision unbounded floats} *)
 168 module UDouble
 169   use real.RealInfix
 170   type udouble
 171
 172   constant eps:real = 0x1p-53 /. (1. +. 0x1p-53)
 173   constant eta:real = 0x1p-1075
 174
 175   clone export UFloat with
 176     type t = udouble,
 177     constant eps = eps,
 178     constant eta = eta,
 179     axiom.
 180 end
 181
 182 (* Helper lemmas to help the proof of propagation lemmas *)
 183 module HelperLemmas
 184   use real.RealInfix
 185   use real.Abs
 186
 187   let ghost div_order_compat (x y z:real)
 188     requires { x <=. y }
 189     requires { 0. <. z }
 190     ensures { x /. z <=. y /. z }
 191     = ()
 192
 193   let ghost div_order_compat2 (x y z:real)
 194     requires { x <=. y }
 195     requires { 0. >. z }
 196     ensures { y /. z <=. x /. z }
 197     = ()
 198
 199   let ghost mult_err (x x_exact x_factor x_rel x_abs y:real)
 200     requires { 0. <=. x_rel }
 201     requires { 0. <=. x_abs }
 202     requires { abs x_exact <=. x_factor }
 203     requires { abs (x -. x_exact) <=. x_rel *. x_factor +. x_abs }
 204     ensures { abs (x *. y -. x_exact *. y) <=. x_rel *. abs (x_factor *. y) +. x_abs *. abs y }
 205   =
 206   assert {
 207     y >=. 0. ->
 208     abs (x *. y -. x_exact *. y) <=. abs (x_rel *. x_factor *. y) +. x_abs *. abs y
 209     by
 210       (x_exact -. x_rel *. x_factor -. x_abs) *. y <=. x *. y <=. (x_exact +. x_rel *. x_factor +. x_abs) *. y
 211   };
 212   assert {
 213     y <. 0. ->
 214     abs (x *. y -. x_exact *. y) <=. abs (x_rel *. x_factor *. y) +. x_abs *. abs y
 215     by
 216       (x_exact +. x_rel *. x_factor +. x_abs) *. y <=. x *. y <=. (x_exact -. x_rel *. x_factor -. x_abs) *. y
 217   }
 218
 219   let ghost mult_err_combine (x x_exact x_factor x_rel x_abs y exact_y y_factor y_rel y_abs:real)
 220     requires { 0. <=. x_rel }
 221     requires { 0. <=. y_rel }
 222     requires { 0. <=. x_abs }
 223     requires { 0. <=. y_abs }
 224     requires { abs x_exact <=. x_factor }
 225     requires { abs exact_y <=. y_factor }
 226     requires { abs (x -. x_exact) <=. x_rel *. x_factor +. x_abs }
 227     requires { abs (y -. exact_y) <=. y_rel *. y_factor +. y_abs }
 228     ensures {
 229       abs (x *. y -. x_exact *. exact_y)
 230       <=. (x_rel +. y_rel +. x_rel *. y_rel) *. (x_factor *. y_factor)
 231         +. (y_abs +. y_abs *. x_rel) *. x_factor
 232         +. (x_abs +. x_abs *. y_rel) *. y_factor
 233         +. x_abs *. y_abs
 234     }
 235   =
 236   mult_err x x_exact x_factor x_rel x_abs y;
 237   mult_err y exact_y y_factor y_rel y_abs x_exact;
 238   mult_err y exact_y y_factor y_rel y_abs x_factor;
 239   assert {
 240     abs (x *. y -. x_exact *. exact_y) <=. (y_rel *. x_factor *. y_factor) +. (y_abs *. x_factor) +. (x_rel *. abs (x_factor *. y)) +. x_abs *. abs y
 241   };
 242   assert {
 243     abs (x *. y -. x_exact *. exact_y) <=. (y_rel *. x_factor *. y_factor) +. (x_rel *. (x_factor *. y_factor *. (1. +. y_rel) +. x_factor *. y_abs)) +. y_abs *. x_factor +. x_abs *. abs y
 244     by
 245       abs (x_factor *. y) <=. x_factor *. y_factor *. (1. +. y_rel) +. x_factor *. y_abs
 246   };
 247   assert {
 248     x_abs *. abs y <=. x_abs *. (y_factor *. (1. +. y_rel) +. y_abs)
 249   }
 250
 251   use real.ExpLog
 252
 253   let ghost exp_approx_err (x x_approx x_factor a b :real)
 254     requires { abs (x_approx -. x) <=. x_factor *. a +. b }
 255     requires { x <=. x_factor }
 256     ensures {
 257       abs (exp(x_approx) -. exp(x)) <=. exp(x) *. (exp(a *. x_factor +. b) -. 1.)
 258     }
 259   =
 260   assert {
 261     exp(x_approx) <=. exp(x) +. exp(x) *. (exp(a *. x_factor +. b) -. 1.)
 262     by
 263       exp (x_approx) <=. exp(x) *. exp (a *. x_factor +. b)
 264   };
 265   assert {
 266     exp(x_approx) >=. exp(x) -. exp(x) *. (exp(a *. x_factor +. b) -. 1.)
 267     by
 268       exp (x_approx) >=. exp(x) *. exp (-. a *. x_factor -. b)
 269     so
 270       exp(x_approx) -. exp(x) >=. exp(x) *. (exp (-. a *. x_factor -. b) -. 1.)
 271     so
 272       exp(a *. x_factor +. b) +. exp(-.a *. x_factor -. b) >=. 2.
 273     so
 274       -. exp(a *. x_factor +. b) +. 1. <=. exp(-.a *. x_factor -. b) -. 1.
 275     so
 276       exp(x) *. ((-. exp(a *. x_factor +. b)) +. 1.) <=. exp(x) *. (exp(-. a *. x_factor -. b) -. 1.)
 277     so
 278       -. exp(x) *. (exp(a *. x_factor +. b) -. 1.) <=. exp(x) *. (exp(-. a *. x_factor -. b) -. 1.)
 279   }
 280
 281   let lemma log_1_minus_x (x:real)
 282     requires { 0. <=. abs x <. 1. }
 283     ensures { log (1. +. x) <=. -. log (1. -. x) }
 284   =
 285     assert { 1. +. x <=. 1. /. (1. -. x) };
 286     assert { forall x y z. 0. <=. x -> 0. <. y -> 0. <=. z -> x *. y <=. z -> x <=. z /. y };
 287     assert { exp (-.log (1. -. x)) = 1. /. (1. -. x) }
 288
 289   let lemma log2_1_minus_x (x:real)
 290     requires { 0. <=. abs x <. 1. }
 291     ensures { log2 (1. +. x) <=. -. log2 (1. -. x) }
 292   =
 293   div_order_compat (log (1. +. x)) (-. log (1. -. x)) (log 2.);
 294   log_1_minus_x x
 295
 296   let lemma log10_1_minus_x (x:real)
 297     requires { 0. <=. abs x <. 1. }
 298     ensures { log10 (1. +. x) <=. -. log10 (1. -. x) }
 299   =
 300   div_order_compat (log (1. +. x)) (-. log (1. -. x)) (log 10.);
 301   log_1_minus_x x
 302
 303   let ghost log_approx_err (x x_approx x_factor a b :real)
 304     requires { abs (x_approx -. x) <=. x_factor *. a +. b }
 305     requires { 0. <. (x -. a *. x_factor -. b) }
 306     requires { 0. <. x <=. x_factor }
 307     ensures {
 308       abs (log x_approx -. log x) <=. -. log(1. -. ((a *. x_factor +. b) /. x))
 309     }
 310   =
 311     assert { a *. x_factor +. b  = x *. ((a *. x_factor +. b) /. x) };
 312     assert {
 313       log (x *. (1. -. (a *. x_factor +. b) /. x))
 314       <=. log x_approx
 315       <=. log (x *. (1. +. (a *. x_factor +. b) /.x))
 316       by
 317         0. <.  (x -. (a *. x_factor +. b)) <=. x_approx
 318     };
 319     log_1_minus_x ((a *. x_factor +. b) /. x)
 320
 321   let ghost log2_approx_err (x x_approx x_factor a b :real)
 322     requires { abs (x_approx -. x) <=. x_factor *. a +. b }
 323     requires { 0. <. (x -. a *. x_factor -. b) }
 324     requires { 0. <. x <=. x_factor }
 325     ensures {
 326       abs (log2 x_approx -. log2 x) <=. -. log2(1. -. ((a *. x_factor +. b) /. x))
 327     }
 328   =
 329     assert { a *. x_factor +. b  = x *. ((a *. x_factor +. b) /. x) };
 330     assert {
 331       log2 (x *. (1. -. (a *. x_factor +. b) /. x))
 332       <=. log2 x_approx
 333       <=. log2 (x *. (1. +. (a *. x_factor +. b) /.x))
 334       by
 335         0. <.  (x -. (a *. x_factor +. b)) <=. x_approx
 336     };
 337     log2_1_minus_x ((a *. x_factor +. b) /. x)
 338
 339   let ghost log10_approx_err (x x_approx x_factor a b :real)
 340     requires { abs (x_approx -. x) <=. x_factor *. a +. b }
 341     requires { 0. <. (x -. a *. x_factor -. b) }
 342     requires { 0. <. x <=. x_factor }
 343     ensures {
 344       abs (log10 x_approx -. log10 x) <=. -. log10(1. -. ((a *. x_factor +. b) /. x))
 345     }
 346   =
 347     assert { a *. x_factor +. b  = x *. ((a *. x_factor +. b) /. x) };
 348     assert {
 349       log10 (x *. (1. -. (a *. x_factor +. b) /. x))
 350       <=. log10 x_approx
 351       <=. log10 (x *. (1. +. (a *. x_factor +. b) /.x))
 352       by
 353         0. <.  (x -. (a *. x_factor +. b)) <=. x_approx
 354     };
 355     log10_1_minus_x ((a *. x_factor +. b) /. x)
 356
 357   use real.Trigonometry
 358
 359   lemma sin_of_approx : forall x y. abs (sin x -. sin y) <=. abs (x -. y)
 360   lemma cos_of_approx : forall x y. abs (cos x -. cos y) <=. abs (x -. y)
 361
 362   use real.Sum
 363   use int.Int
 364   use real.FromInt
 365
 366   let rec ghost sum_approx_err (fi_rel fi_abs:real) (f f_exact f_factor : int -> real) (a b:int)
 367     requires { a <= b }
 368     requires { forall i. a <= i < b -> abs (f i -. f_exact i) <=. f_factor i *. fi_rel +. fi_abs }
 369     variant { b - a }
 370     ensures { abs (sum f a b -. sum f_exact a b) <=. fi_rel *. sum f_factor a b +. fi_abs *. from_int (b-a) }
 371   =
 372   if (a < b) then
 373     begin
 374       sum_approx_err fi_rel fi_abs f f_exact f_factor a (b - 1)
 375     end
 376
 377 end
 378
 379 (** {4 Single propagation lemmas} *)
 380 module USingleLemmas
 381   use real.RealInfix
 382   use real.FromInt
 383   use real.Abs
 384   use USingle
 385
 386   let lemma uadd_single_error_propagation (x_f y_f r: usingle) (x x_factor x_rel x_abs y y_factor y_rel y_abs : real)
 387     requires {
 388       abs (to_real x_f -. x) <=. x_rel *. x_factor +. x_abs
 389     }
 390     requires {
 391       abs (to_real y_f -. y) <=. y_rel *. y_factor +. y_abs
 392     }
 393     requires { abs x <=. x_factor }
 394     requires { abs y <=. y_factor }
 395     (* TODO: Use (0 <=. x_rel \/ (x_factor = 0 /\ x_abs = 0)), same for y. *)
 396     requires { 0. <=. x_rel }
 397     requires { 0. <=. y_rel }
 398     requires { 0. <=. x_abs }
 399     requires { 0. <=. y_abs }
 400     requires { r = (x_f ++. y_f) }
 401     ensures {
 402       abs (to_real r -. (x +. y)) <=.
 403       (x_rel +. y_rel +. eps) *. (x_factor +. y_factor)
 404           +. ((1. +. eps +. y_rel) *. x_abs +. (1. +. eps +. x_rel) *. y_abs)
 405     }
 406   =
 407   let ghost delta = abs (to_real (x_f ++. y_f) -. (to_real x_f +. to_real y_f)) in
 408   assert {
 409     0. <=. x_rel /\ 0. <=. y_rel ->
 410     delta <=.
 411       (eps +. y_rel) *. x_factor +. (eps +. x_rel) *. y_factor
 412       +. (x_rel +. eps) *. y_abs +. (y_rel +. eps) *. x_abs
 413     by
 414       (delta <=. x_factor *. x_rel +. x_abs +. x_factor
 415       so
 416         x_factor +. x_abs <=. eps *. (y_factor +. y_abs) ->
 417         delta <=. (eps +. x_rel) *. y_factor
 418         +. (eps +. y_rel) *. x_factor
 419         +. (y_rel +. eps) *. x_abs +. (x_rel +. eps) *. y_abs
 420       by
 421         delta <=. eps *. (y_factor +. y_abs) *. x_rel
 422               +. (eps *. (y_factor +. y_abs)))
 423       /\
 424       (delta <=. y_factor *. y_rel +. y_abs +. y_factor
 425       so
 426       abs y_factor +. y_abs <=. eps *. (x_factor +. x_abs) ->
 427       delta <=. (eps +. y_rel) *. x_factor
 428         +. (eps +. x_rel) *. y_factor
 429         +. (x_rel +. eps) *. y_abs +. (y_rel +. eps) *. x_abs
 430       by
 431         delta <=. eps *. (x_factor +. x_abs) *. y_rel
 432               +. (eps *. (x_factor +. x_abs)))
 433       /\
 434       (
 435        (eps *. (x_factor +. x_abs) <. abs y_factor +. y_abs /\
 436        eps *. (y_factor +. y_abs) <. abs x_factor +. x_abs) ->
 437        (delta <=.
 438        (eps +. y_rel) *. x_factor +. (eps +. x_rel) *. y_factor
 439       +. (x_rel +. eps) *. y_abs +. (y_rel +. eps) *. x_abs
 440       by
 441         abs (to_real x_f +. to_real y_f) <=.
 442         abs (to_real x_f -. x) +. x_factor +. (abs (to_real y_f -. y) +. y_factor)
 443       so
 444         x_factor *. x_rel <=. (y_factor +. y_abs) /. eps *. x_rel /\
 445         y_factor *. y_rel <=. (x_factor +. x_abs) /. eps *. y_rel))
 446   }
 447
 448   let lemma usub_single_error_propagation (x_f y_f r : usingle) (x x_factor x_rel x_abs y y_factor y_rel y_abs : real)
 449     requires {
 450       abs (to_real x_f -. x) <=. x_rel *. x_factor +. x_abs
 451     }
 452     requires {
 453       abs (to_real y_f -. y) <=. y_rel *. y_factor +. y_abs
 454     }
 455     requires { abs x <=. x_factor }
 456     requires { abs y <=. y_factor }
 457     requires { 0. <=. x_abs }
 458     requires { 0. <=. y_abs }
 459     requires { 0. <=. x_rel }
 460     requires { 0. <=. y_rel }
 461     requires { r = x_f --. y_f }
 462     ensures {
 463       abs (to_real r -. (x -. y))
 464       <=. (x_rel +. y_rel +. eps) *. (x_factor +. y_factor)
 465           +. ((1. +. eps +. y_rel) *. x_abs +. (1. +. eps +. x_rel) *. y_abs)
 466     }
 467   = uadd_single_error_propagation x_f (--. y_f) r x x_factor x_rel x_abs (-. y) y_factor y_rel y_abs
 468
 469   use HelperLemmas
 470
 471   let lemma umul_single_error_propagation (x_f y_f r : usingle) (x x_factor x_rel x_abs y y_factor y_rel y_abs : real)
 472     requires {
 473       abs (to_real x_f -. x) <=. x_rel *. x_factor +. x_abs
 474     }
 475     requires {
 476       abs (to_real y_f -. y) <=. y_rel *. y_factor +. y_abs
 477     }
 478     requires { abs x <=. x_factor }
 479     requires { abs y <=. y_factor }
 480     requires { 0. <=. x_rel }
 481     requires { 0. <=. y_rel }
 482     requires { 0. <=. x_abs }
 483     requires { 0. <=. y_abs }
 484     requires { r = x_f **. y_f }
 485     ensures {
 486       abs (to_real r -. (x *. y)) <=.
 487         (eps +. (x_rel +. y_rel +. x_rel *. y_rel) *. (1. +. eps)) *. (x_factor *. y_factor)
 488         +. (((y_abs +. y_abs *. x_rel) *. x_factor
 489         +. (x_abs +. x_abs *. y_rel) *. y_factor
 490         +. x_abs *. y_abs) *. (1. +. eps) +. eta)
 491     }
 492   =
 493   assert {
 494     to_real x_f *. to_real y_f -. abs (to_real x_f *. to_real y_f) *. eps -. eta
 495     <=. to_real (x_f **. y_f)
 496     <=. to_real x_f *. to_real y_f +. abs (to_real x_f *. to_real y_f) *. eps +. eta
 497   };
 498     assert { abs (x *. y) <=. x_factor *. y_factor by
 499        abs x *. abs y <=. x_factor *. abs y = abs y *. x_factor <=. y_factor *. x_factor };
 500   mult_err_combine (to_real x_f) x x_factor x_rel x_abs (to_real y_f) y y_factor y_rel y_abs
 501
 502   use real.ExpLog
 503
 504   let lemma log_single_error_propagation (logx_f x_f : usingle)
 505         (x_exact x_factor log_rel log_abs x_rel x_abs : real)
 506     requires { abs (to_real x_f -. x_exact) <=. x_rel *. x_factor +. x_abs }
 507     requires {
 508       abs (to_real logx_f -. log(to_real x_f))
 509        <=. log_rel *. abs (log (to_real x_f)) +. log_abs
 510     }
 511     requires { 0. <. x_exact <=. x_factor }
 512     requires { 0. <. (x_exact -. x_rel *. x_factor -. x_abs) }
 513     requires { 0. <=. log_rel }
 514     ensures {
 515       abs (to_real logx_f -. log (x_exact))
 516         <=. log_rel *. abs (log (x_exact)) +.
 517           (-. log (1. -. ((x_rel *. x_factor +. x_abs) /. x_exact)) *. (1. +. log_rel)
 518           +. log_abs)
 519     }
 520   =
 521   log_approx_err x_exact (to_real x_f) x_factor x_rel x_abs;
 522   assert {
 523    abs (log (to_real x_f)) *. log_rel
 524     <=. (abs (log (x_exact)) -. log (1.0 -. (((x_rel *. x_factor) +. x_abs) /. x_exact))) *. log_rel
 525   }
 526
 527   let lemma log2_single_error_propagation (log2x_f x_f : usingle)
 528         (x_exact x_factor log_rel log_abs x_rel x_abs : real)
 529     requires { abs (to_real x_f -. x_exact) <=. x_rel *. x_factor +. x_abs }
 530     requires {
 531       abs (to_real log2x_f -. log2(to_real x_f))
 532        <=. log_rel *. abs (log2 (to_real x_f)) +. log_abs
 533     }
 534     requires { 0. <. x_exact <=. x_factor }
 535     requires { 0. <. (x_exact -. x_rel *. x_factor -. x_abs) }
 536     requires { 0. <=. log_rel }
 537     ensures {
 538       abs (to_real log2x_f -. log2 (x_exact))
 539         <=. log_rel *. abs (log2 (x_exact)) +.
 540           (-. log2 (1. -. ((x_rel *. x_factor +. x_abs) /. x_exact)) *. (1. +. log_rel)
 541           +. log_abs)
 542     }
 543   =
 544   log2_approx_err x_exact (to_real x_f) x_factor x_rel x_abs;
 545   assert {
 546    abs (log2 (to_real x_f)) *. log_rel
 547     <=. (abs (log2 (x_exact)) -. log2 (1.0 -. (((x_rel *. x_factor) +. x_abs) /. x_exact))) *. log_rel
 548   }
 549
 550   let lemma log10_single_error_propagation (log10x_f x_f : usingle)
 551         (x_exact x_factor log_rel log_abs x_rel x_abs : real)
 552     requires { abs (to_real x_f -. x_exact) <=. x_rel *. x_factor +. x_abs }
 553     requires {
 554       abs (to_real log10x_f -. log10(to_real x_f))
 555        <=. log_rel *. abs (log10 (to_real x_f)) +. log_abs
 556     }
 557     requires { 0. <. x_exact <=. x_factor }
 558     requires { 0. <. (x_exact -. x_rel *. x_factor -. x_abs) }
 559     requires { 0. <=. log_rel }
 560     ensures {
 561       abs (to_real log10x_f -. log10 (x_exact))
 562         <=. log_rel *. abs (log10 (x_exact)) +.
 563           (-. log10 (1. -. ((x_rel *. x_factor +. x_abs) /. x_exact)) *. (1. +. log_rel)
 564           +. log_abs)
 565     }
 566   =
 567   log10_approx_err x_exact (to_real x_f) x_factor x_rel x_abs;
 568   assert {
 569    abs (log10 (to_real x_f)) *. log_rel
 570     <=. (abs (log10 (x_exact)) -. log10 (1.0 -. (((x_rel *. x_factor) +. x_abs) /. x_exact))) *. log_rel
 571   }
 572
 573   let lemma exp_single_error_propagation (expx_f x_f : usingle)
 574         (x_exact x_factor exp_rel exp_abs x_rel x_abs : real)
 575     requires { abs (to_real x_f -. x_exact) <=. x_rel *. x_factor +. x_abs }
 576     requires {
 577       abs (to_real expx_f -. exp(to_real x_f))
 578         <=. exp_rel *. exp (to_real x_f) +. exp_abs
 579     }
 580     requires { x_exact <=. x_factor }
 581     requires { 0. <=. exp_rel <=. 1. }
 582     ensures {
 583       abs (to_real expx_f -. exp (x_exact))
 584       <=. (exp_rel +. (exp(x_rel *. x_factor +. x_abs) -. 1.) *. (1. +. exp_rel)) *. exp(x_exact)
 585         +. exp_abs
 586     }
 587   =
 588     exp_approx_err x_exact (to_real x_f) x_factor x_rel x_abs;
 589     assert {
 590       exp x_exact *. (1. -. exp_rel) -.
 591       exp x_exact *. (exp (x_rel *. x_factor +. x_abs) -. 1.) *. (1. -. exp_rel)
 592       -. exp_abs
 593       <=. to_real expx_f
 594       by
 595         (exp x_exact -. exp x_exact *. (exp (x_rel *. x_factor +. x_abs) -. 1.))
 596          *. (1. -. exp_rel) -. exp_abs
 597         <=. exp (to_real x_f) *. (1. -. exp_rel) -. exp_abs
 598         <=. to_real expx_f
 599     };
 600     assert {
 601       to_real expx_f <=. (exp(x_exact) +. exp(x_exact)*.(exp(x_rel *. x_factor +. x_abs) -. 1.))*. (1. +. exp_rel) +. exp_abs
 602       by
 603         to_real expx_f <=. exp(to_real x_f) *. (1. +. exp_rel) +. exp_abs
 604     };
 605
 606
 607   use real.Trigonometry
 608
 609   let lemma sin_single_error_propagation (sinx_f x_f : usingle)
 610         (x_exact x_factor sin_rel sin_abs x_rel x_abs : real)
 611     requires { abs (to_real x_f -. x_exact) <=. x_rel *. x_factor +. x_abs }
 612     requires {
 613       abs (to_real sinx_f -. sin(to_real x_f))
 614         <=. sin_rel *. abs (sin (to_real x_f)) +. sin_abs
 615     }
 616     requires { x_exact <=. x_factor }
 617     requires { 0. <=. sin_rel }
 618     ensures {
 619       abs (to_real sinx_f -. sin (x_exact))
 620       <=. sin_rel *. abs(sin(x_exact))
 621           +. (((x_rel *. x_factor +. x_abs) *. (1. +. sin_rel)) +. sin_abs)
 622     }
 623   =
 624   assert {
 625   abs (sin (to_real x_f)) *. sin_rel
 626   <=. (abs (sin x_exact) +. (x_rel *. x_factor +. x_abs)) *. sin_rel
 627   }
 628
 629   let lemma cos_single_error_propagation (cosx_f x_f : usingle)
 630         (x_exact x_factor cos_rel cos_abs x_rel x_abs : real)
 631     requires { abs (to_real x_f -. x_exact) <=. x_rel *. x_factor +. x_abs }
 632     requires {
 633       abs (to_real cosx_f -. cos(to_real x_f))
 634         <=. cos_rel *. abs (cos (to_real x_f)) +. cos_abs
 635     }
 636     requires { x_exact <=. x_factor }
 637     requires { 0. <=. cos_rel }
 638     ensures {
 639       abs (to_real cosx_f -. cos (x_exact))
 640       <=. cos_rel *. abs(cos(x_exact))
 641           +. (((x_rel *. x_factor +. x_abs) *. (1. +. cos_rel)) +. cos_abs)
 642     }
 643   =
 644   assert {
 645   abs (cos (to_real x_f)) *. cos_rel
 646   <=. (abs (cos x_exact) +. (x_rel *. x_factor +. x_abs)) *. cos_rel
 647   }
 648
 649   use real.Sum
 650   use int.Int
 651   use real.FromInt
 652
 653   function real_fun (f:int -> usingle) : int -> real = fun i -> to_real (f i)
 654
 655   let lemma sum_single_error_propagation (x : usingle)
 656                 (f : int -> usingle) (f_exact f_factor f_factor' : int -> real) (n:int)
 657                 (sum_rel sum_abs f_rel f_abs : real)
 658     requires {
 659       forall i. 0 <= i < n ->
 660         abs ((real_fun f) i -. f_exact i) <=. f_rel *. f_factor i +. f_abs
 661     }
 662     requires {
 663       forall i. 0 <= i < n ->
 664       f_factor i -. f_rel *. f_factor i -. f_abs <=. f_factor' i <=. f_factor i +. f_rel *. f_factor i +. f_abs
 665     }
 666     requires {
 667       abs (to_real x -. (sum (real_fun f) 0 n))
 668         <=. sum_rel *. (sum f_factor' 0 n) +. sum_abs
 669     }
 670     requires { 0. <=. sum_rel }
 671     requires { 0 <= n }
 672     ensures {
 673       abs (to_real x -. sum f_exact 0 n)
 674       <=. (f_rel +. (sum_rel *. (1. +. f_rel))) *. sum f_factor 0 n +.
 675         ((f_abs *. from_int n *.(1. +. sum_rel)) +. sum_abs)
 676     }
 677   =
 678   sum_approx_err f_rel f_abs (real_fun f) f_exact f_factor 0 n;
 679   sum_approx_err f_rel f_abs f_factor' f_factor f_factor 0 n;
 680   assert {
 681     sum_rel *. sum f_factor' 0 n <=.
 682     sum_rel *. (sum f_factor 0 n +. ((f_rel *. sum f_factor 0 n) +. (f_abs *. from_int n)))
 683   }
 684
 685   (* We don't have an error on y_f because in practice we won't have an exact division with an approximate divisor *)
 686   let lemma udiv_exact_single_error_propagation (x_f y_f r: usingle) (x x_factor x_rel x_abs : real)
 687     requires {
 688       abs (to_real x_f -. x) <=. x_rel *. x_factor +. x_abs
 689     }
 690     requires { abs x <=. x_factor }
 691     requires { 0. <=. x_rel }
 692     requires { 0. <=. x_abs }
 693     requires { 0. <> to_real y_f }
 694     requires { is_exact udiv x_f y_f }
 695     requires { r = x_f ///. y_f }
 696     ensures {
 697       abs (to_real r -. (x /. (to_real y_f))) <=.
 698         x_rel *. (x_factor /. abs (to_real y_f)) +. ((x_abs /. abs (to_real y_f)) +. eta)
 699     }
 700   =
 701   let lemma y_f_pos ()
 702     requires { 0. <. to_real y_f }
 703     ensures {
 704       abs (to_real r -. (x /. (to_real y_f))) <=.
 705         x_rel *. (x_factor /. to_real y_f) +. ((x_abs /. to_real y_f) +. eta)
 706     }
 707     =
 708     div_order_compat (to_real x_f) (x +. x_rel *. x_factor +. x_abs) (to_real y_f);
 709     div_order_compat (x -. x_rel *. x_factor -. x_abs) (to_real x_f) (to_real y_f) in
 710   let lemma y_f_neg ()
 711     requires { to_real y_f <. 0. }
 712     ensures {
 713       abs (to_real r -. (x /. (to_real y_f))) <=.
 714         x_rel *. (x_factor /. abs (to_real y_f)) +. ((x_abs /. abs (to_real y_f)) +. eta)
 715     }
 716     =
 717     div_order_compat2 (to_real x_f) (x +. x_rel *. x_factor +. x_abs) (to_real y_f);
 718     (* TODO: Prove this somehow *)
 719     assert {
 720       forall x y. y <> 0.0 -> x /. y <=. abs x /. abs y
 721       by abs (x /. y) = abs (x *. inv y) = abs x *. abs (inv y) = abs x *. inv (abs y) = abs x /. abs y
 722     };
 723     assert {
 724       (x -. x_rel *. x_factor -. x_abs) /. to_real y_f
 725       <=. x /. (to_real y_f) +. ((x_rel *. x_factor) +. x_abs) /. abs (to_real y_f)
 726       by
 727         (-. x_rel *. x_factor -. x_abs) /. to_real y_f
 728         <=. (x_rel *. x_factor +. x_abs) /. abs (to_real y_f)
 729     };
 730     div_order_compat2 (x -. x_rel *. x_factor -. x_abs) (to_real x_f) (to_real y_f);
 731   in ()
 732 end
 733
 734 (** {5 Double propagation lemmas} *)
 735 module UDoubleLemmas
 736   use real.RealInfix
 737   use real.FromInt
 738   use real.Abs
 739   use UDouble
 740
 741   let lemma uadd_double_error_propagation (x_f y_f r : udouble) (x x_factor x_rel x_abs y y_factor y_rel y_abs : real)
 742     requires {
 743       abs (to_real x_f -. x) <=. x_rel *. x_factor +. x_abs
 744     }
 745     requires {
 746       abs (to_real y_f -. y) <=. y_rel *. y_factor +. y_abs
 747     }
 748     requires { abs x <=. x_factor }
 749     requires { abs y <=. y_factor }
 750     (* TODO: Use (0 <=. x_rel \/ (x_factor = 0 /\ x_abs = 0)), same for y. *)
 751     requires { 0. <=. x_rel }
 752     requires { 0. <=. y_rel }
 753     requires { 0. <=. x_abs }
 754     requires { 0. <=. y_abs }
 755     requires { r = x_f ++. y_f }
 756     ensures {
 757       abs (to_real r -. (x +. y)) <=.
 758       (x_rel +. y_rel +. eps) *. (x_factor +. y_factor)
 759           +. ((1. +. eps +. y_rel) *. x_abs +. (1. +. eps +. x_rel) *. y_abs)
 760     }
 761   =
 762   let ghost delta = abs (to_real (x_f ++. y_f) -. (to_real x_f +. to_real y_f)) in
 763   assert {
 764     0. <=. x_rel /\ 0. <=. y_rel ->
 765     delta <=.
 766       (eps +. y_rel) *. x_factor +. (eps +. x_rel) *. y_factor
 767       +. (x_rel +. eps) *. y_abs +. (y_rel +. eps) *. x_abs
 768     by
 769       (delta <=. x_factor *. x_rel +. x_abs +. x_factor
 770       so
 771         x_factor +. x_abs <=. eps *. (y_factor +. y_abs) ->
 772         delta <=. (eps +. x_rel) *. y_factor
 773         +. (eps +. y_rel) *. x_factor
 774         +. (y_rel +. eps) *. x_abs +. (x_rel +. eps) *. y_abs
 775       by
 776         delta <=. eps *. (y_factor +. y_abs) *. x_rel
 777               +. (eps *. (y_factor +. y_abs)))
 778       /\
 779       (delta <=. y_factor *. y_rel +. y_abs +. y_factor
 780       so
 781       abs y_factor +. y_abs <=. eps *. (x_factor +. x_abs) ->
 782       delta <=. (eps +. y_rel) *. x_factor
 783         +. (eps +. x_rel) *. y_factor
 784         +. (x_rel +. eps) *. y_abs +. (y_rel +. eps) *. x_abs
 785       by
 786         delta <=. eps *. (x_factor +. x_abs) *. y_rel
 787               +. (eps *. (x_factor +. x_abs)))
 788       /\
 789       (
 790        (eps *. (x_factor +. x_abs) <. abs y_factor +. y_abs /\
 791        eps *. (y_factor +. y_abs) <. abs x_factor +. x_abs) ->
 792        (delta <=.
 793        (eps +. y_rel) *. x_factor +. (eps +. x_rel) *. y_factor
 794       +. (x_rel +. eps) *. y_abs +. (y_rel +. eps) *. x_abs
 795       by
 796         abs (to_real x_f +. to_real y_f) <=.
 797         abs (to_real x_f -. x) +. x_factor +. (abs (to_real y_f -. y) +. y_factor)
 798       so
 799         x_factor *. x_rel <=. (y_factor +. y_abs) /. eps *. x_rel /\
 800         y_factor *. y_rel <=. (x_factor +. x_abs) /. eps *. y_rel))
 801   }
 802
 803   let lemma usub_double_error_propagation (x_f y_f r : udouble) (x x_factor x_rel x_abs y y_factor y_rel y_abs : real)
 804     requires {
 805       abs (to_real x_f -. x) <=. x_rel *. x_factor +. x_abs
 806     }
 807     requires {
 808       abs (to_real y_f -. y) <=. y_rel *. y_factor +. y_abs
 809     }
 810     requires { abs x <=. x_factor }
 811     requires { abs y <=. y_factor }
 812     requires { 0. <=. x_abs }
 813     requires { 0. <=. y_abs }
 814     requires { 0. <=. x_rel }
 815     requires { 0. <=. y_rel }
 816     requires { r = x_f --. y_f }
 817     ensures {
 818       abs (to_real r -. (x -. y))
 819       <=. (x_rel +. y_rel +. eps) *. (x_factor +. y_factor)
 820           +. ((1. +. eps +. y_rel) *. x_abs +. (1. +. eps +. x_rel) *. y_abs)
 821     }
 822   = uadd_double_error_propagation x_f (--. y_f) r x x_factor x_rel x_abs (-. y) y_factor y_rel y_abs
 823
 824   use HelperLemmas
 825
 826   let lemma umul_double_error_propagation (x_f y_f r : udouble) (x x_factor x_rel x_abs y y_factor y_rel y_abs : real)
 827     requires {
 828       abs (to_real x_f -. x) <=. x_rel *. x_factor +. x_abs
 829     }
 830     requires {
 831       abs (to_real y_f -. y) <=. y_rel *. y_factor +. y_abs
 832     }
 833     requires { abs x <=. x_factor }
 834     requires { abs y <=. y_factor }
 835     requires { 0. <=. x_rel }
 836     requires { 0. <=. y_rel }
 837     requires { 0. <=. x_abs }
 838     requires { 0. <=. y_abs }
 839     requires { r = x_f **. y_f }
 840     ensures {
 841       abs (to_real r -. (x *. y)) <=.
 842         (eps +. (x_rel +. y_rel +. x_rel *. y_rel) *. (1. +. eps)) *. (x_factor *. y_factor)
 843         +. (((y_abs +. y_abs *. x_rel) *. x_factor
 844         +. (x_abs +. x_abs *. y_rel) *. y_factor
 845         +. x_abs *. y_abs) *. (1. +. eps) +. eta)
 846     }
 847   =
 848   assert {
 849     to_real x_f *. to_real y_f -. abs (to_real x_f *. to_real y_f) *. eps -. eta
 850     <=. to_real (x_f **. y_f)
 851     <=. to_real x_f *. to_real y_f +. abs (to_real x_f *. to_real y_f) *. eps +. eta
 852   };
 853     assert { abs (x *. y) <=. x_factor *. y_factor by
 854        abs x *. abs y <=. x_factor *. abs y = abs y *. x_factor <=. y_factor *. x_factor };
 855   mult_err_combine (to_real x_f) x x_factor x_rel x_abs (to_real y_f) y y_factor y_rel y_abs
 856
 857   use real.ExpLog
 858
 859   let lemma log_double_error_propagation (logx_f x_f : udouble)
 860         (x_exact x_factor log_rel log_abs x_rel x_abs : real)
 861     requires { abs (to_real x_f -. x_exact) <=. x_rel *. x_factor +. x_abs }
 862     requires {
 863       abs (to_real logx_f -. log(to_real x_f))
 864        <=. log_rel *. abs (log (to_real x_f)) +. log_abs
 865     }
 866     requires { 0. <. x_exact <=. x_factor }
 867     requires { 0. <. (x_exact -. x_rel *. x_factor -. x_abs) }
 868     requires { 0. <=. log_rel }
 869     ensures {
 870       abs (to_real logx_f -. log (x_exact))
 871         <=. log_rel *. abs (log (x_exact)) +.
 872           (-. log (1. -. ((x_rel *. x_factor +. x_abs) /. x_exact)) *. (1. +. log_rel)
 873           +. log_abs)
 874     }
 875   =
 876   log_approx_err x_exact (to_real x_f) x_factor x_rel x_abs;
 877   assert {
 878    abs (log (to_real x_f)) *. log_rel
 879     <=. (abs (log (x_exact)) -. log (1.0 -. (((x_rel *. x_factor) +. x_abs) /. x_exact))) *. log_rel
 880   }
 881
 882   let lemma log2_double_error_propagation (log2x_f x_f : udouble)
 883         (x_exact x_factor log_rel log_abs x_rel x_abs : real)
 884     requires { abs (to_real x_f -. x_exact) <=. x_rel *. x_factor +. x_abs }
 885     requires {
 886       abs (to_real log2x_f -. log2(to_real x_f))
 887        <=. log_rel *. abs (log2 (to_real x_f)) +. log_abs
 888     }
 889     requires { 0. <. x_exact <=. x_factor }
 890     requires { 0. <. (x_exact -. x_rel *. x_factor -. x_abs) }
 891     requires { 0. <=. log_rel }
 892     ensures {
 893       abs (to_real log2x_f -. log2 (x_exact))
 894         <=. log_rel *. abs (log2 (x_exact)) +.
 895           (-. log2 (1. -. ((x_rel *. x_factor +. x_abs) /. x_exact)) *. (1. +. log_rel)
 896           +. log_abs)
 897     }
 898   =
 899   log2_approx_err x_exact (to_real x_f) x_factor x_rel x_abs;
 900   assert {
 901    abs (log2 (to_real x_f)) *. log_rel
 902     <=. (abs (log2 (x_exact)) -. log2 (1.0 -. (((x_rel *. x_factor) +. x_abs) /. x_exact))) *. log_rel
 903   }
 904
 905   let lemma log10_double_error_propagation (log10x_f x_f : udouble)
 906         (x_exact x_factor log_rel log_abs x_rel x_abs : real)
 907     requires { abs (to_real x_f -. x_exact) <=. x_rel *. x_factor +. x_abs }
 908     requires {
 909       abs (to_real log10x_f -. log10(to_real x_f))
 910        <=. log_rel *. abs (log10 (to_real x_f)) +. log_abs
 911     }
 912     requires { 0. <. x_exact <=. x_factor }
 913     requires { 0. <. (x_exact -. x_rel *. x_factor -. x_abs) }
 914     requires { 0. <=. log_rel }
 915     ensures {
 916       abs (to_real log10x_f -. log10 (x_exact))
 917         <=. log_rel *. abs (log10 (x_exact)) +.
 918           (-. log10 (1. -. ((x_rel *. x_factor +. x_abs) /. x_exact)) *. (1. +. log_rel)
 919           +. log_abs)
 920     }
 921   =
 922   log10_approx_err x_exact (to_real x_f) x_factor x_rel x_abs;
 923   assert {
 924    abs (log10 (to_real x_f)) *. log_rel
 925     <=. (abs (log10 (x_exact)) -. log10 (1.0 -. (((x_rel *. x_factor) +. x_abs) /. x_exact))) *. log_rel
 926   }
 927
 928   let lemma exp_double_error_propagation (expx_f x_f : udouble)
 929         (x_exact x_factor exp_rel exp_abs x_rel x_abs : real)
 930     requires { abs (to_real x_f -. x_exact) <=. x_rel *. x_factor +. x_abs }
 931     requires {
 932       abs (to_real expx_f -. exp(to_real x_f))
 933         <=. exp_rel *. exp (to_real x_f) +. exp_abs
 934     }
 935     requires { x_exact <=. x_factor }
 936     requires { 0. <=. exp_rel <=. 1. }
 937     ensures {
 938       abs (to_real expx_f -. exp (x_exact))
 939       <=. (exp_rel +. (exp(x_rel *. x_factor +. x_abs) -. 1.) *. (1. +. exp_rel)) *. exp(x_exact)
 940         +. exp_abs
 941     }
 942   =
 943     exp_approx_err x_exact (to_real x_f) x_factor x_rel x_abs;
 944     assert {
 945       exp x_exact *. (1. -. exp_rel) -.
 946       exp x_exact *. (exp (x_rel *. x_factor +. x_abs) -. 1.) *. (1. -. exp_rel)
 947       -. exp_abs
 948       <=. to_real expx_f
 949       by
 950         (exp x_exact -. exp x_exact *. (exp (x_rel *. x_factor +. x_abs) -. 1.))
 951          *. (1. -. exp_rel) -. exp_abs
 952         <=. exp (to_real x_f) *. (1. -. exp_rel) -. exp_abs
 953         <=. to_real expx_f
 954     };
 955     assert {
 956       to_real expx_f <=. (exp(x_exact) +. exp(x_exact)*.(exp(x_rel *. x_factor +. x_abs) -. 1.))*. (1. +. exp_rel) +. exp_abs
 957       by
 958         to_real expx_f <=. exp(to_real x_f) *. (1. +. exp_rel) +. exp_abs
 959     };
 960
 961
 962   use real.Trigonometry
 963
 964   let lemma sin_double_error_propagation (sinx_f x_f : udouble)
 965         (x_exact x_factor sin_rel sin_abs x_rel x_abs : real)
 966     requires { abs (to_real x_f -. x_exact) <=. x_rel *. x_factor +. x_abs }
 967     requires {
 968       abs (to_real sinx_f -. sin(to_real x_f))
 969         <=. sin_rel *. abs (sin (to_real x_f)) +. sin_abs
 970     }
 971     requires { x_exact <=. x_factor }
 972     requires { 0. <=. sin_rel }
 973     ensures {
 974       abs (to_real sinx_f -. sin (x_exact))
 975       <=. sin_rel *. abs(sin(x_exact))
 976           +. (((x_rel *. x_factor +. x_abs) *. (1. +. sin_rel)) +. sin_abs)
 977     }
 978   =
 979   assert {
 980   abs (sin (to_real x_f)) *. sin_rel
 981   <=. (abs (sin x_exact) +. (x_rel *. x_factor +. x_abs)) *. sin_rel
 982   }
 983
 984   let lemma cos_double_error_propagation (cosx_f x_f : udouble)
 985         (x_exact x_factor cos_rel cos_abs x_rel x_abs : real)
 986     requires { abs (to_real x_f -. x_exact) <=. x_rel *. x_factor +. x_abs }
 987     requires {
 988       abs (to_real cosx_f -. cos(to_real x_f))
 989         <=. cos_rel *. abs (cos (to_real x_f)) +. cos_abs
 990     }
 991     requires { x_exact <=. x_factor }
 992     requires { 0. <=. cos_rel }
 993     ensures {
 994       abs (to_real cosx_f -. cos (x_exact))
 995       <=. cos_rel *. abs(cos(x_exact))
 996           +. (((x_rel *. x_factor +. x_abs) *. (1. +. cos_rel)) +. cos_abs)
 997     }
 998   =
 999   assert {
1000   abs (cos (to_real x_f)) *. cos_rel
1001   <=. (abs (cos x_exact) +. (x_rel *. x_factor +. x_abs)) *. cos_rel
1002   }
1003
1004   use real.Sum
1005   use int.Int
1006   use real.FromInt
1007
1008   function real_fun (f:int -> udouble) : int -> real = fun i -> to_real (f i)
1009
1010   let lemma sum_double_error_propagation (x : udouble)
1011                 (f : int -> udouble) (f_exact f_factor f_factor' : int -> real) (n:int)
1012                 (sum_rel sum_abs f_rel f_abs : real)
1013     requires {
1014       forall i. 0 <= i < n ->
1015         abs ((real_fun f) i -. f_exact i) <=. f_rel *. f_factor i +. f_abs
1016     }
1017     requires {
1018       forall i. 0 <= i < n ->
1019       f_factor i -. f_rel *. f_factor i -. f_abs <=. f_factor' i <=. f_factor i +. f_rel *. f_factor i +. f_abs
1020     }
1021     requires {
1022       abs (to_real x -. (sum (real_fun f) 0 n))
1023         <=. sum_rel *. (sum f_factor' 0 n) +. sum_abs
1024     }
1025     requires { 0. <=. sum_rel }
1026     requires { 0 <= n }
1027     ensures {
1028       abs (to_real x -. sum f_exact 0 n)
1029       <=. (f_rel +. (sum_rel *. (1. +. f_rel))) *. sum f_factor 0 n +.
1030         ((f_abs *. from_int n *.(1. +. sum_rel)) +. sum_abs)
1031     }
1032   =
1033   sum_approx_err f_rel f_abs (real_fun f) f_exact f_factor 0 n;
1034   sum_approx_err f_rel f_abs f_factor' f_factor f_factor 0 n;
1035   assert {
1036     sum_rel *. sum f_factor' 0 n <=.
1037     sum_rel *. (sum f_factor 0 n +. ((f_rel *. sum f_factor 0 n) +. (f_abs *. from_int n)))
1038   }
1039
1040   (* We don't have an error on y_f because in practice we won't have an exact division with an approximate divisor *)
1041   let lemma udiv_exact_single_error_propagation (x_f y_f r: udouble) (x x_factor x_rel x_abs : real)
1042     requires {
1043       abs (to_real x_f -. x) <=. x_rel *. x_factor +. x_abs
1044     }
1045     requires { abs x <=. x_factor }
1046     requires { 0. <=. x_rel }
1047     requires { 0. <=. x_abs }
1048     requires { 0. <> to_real y_f }
1049     requires { is_exact udiv x_f y_f }
1050     requires { r = x_f ///. y_f }
1051     ensures {
1052       abs (to_real r -. (x /. (to_real y_f))) <=.
1053         x_rel *. (x_factor /. abs (to_real y_f)) +. ((x_abs /. abs (to_real y_f)) +. eta)
1054     }
1055   =
1056   let lemma y_f_pos ()
1057     requires { 0. <. to_real y_f }
1058     ensures {
1059       abs (to_real r -. (x /. (to_real y_f))) <=.
1060         x_rel *. (x_factor /. to_real y_f) +. ((x_abs /. to_real y_f) +. eta)
1061     }
1062     =
1063     div_order_compat (to_real x_f) (x +. x_rel *. x_factor +. x_abs) (to_real y_f);
1064     div_order_compat (x -. x_rel *. x_factor -. x_abs) (to_real x_f) (to_real y_f) in
1065   let lemma y_f_neg ()
1066     requires { to_real y_f <. 0. }
1067     ensures {
1068       abs (to_real r -. (x /. (to_real y_f))) <=.
1069         x_rel *. (x_factor /. abs (to_real y_f)) +. ((x_abs /. abs (to_real y_f)) +. eta)
1070     }
1071     =
1072     div_order_compat2 (to_real x_f) (x +. x_rel *. x_factor +. x_abs) (to_real y_f);
1073     (* TODO: Prove this somehow *)
1074     assert {
1075       forall x y. x /. y <=. abs x /. abs y
1076     };
1077     assert {
1078       (x -. x_rel *. x_factor -. x_abs) /. to_real y_f
1079       <=. x /. (to_real y_f) +. ((x_rel *. x_factor) +. x_abs) /. abs (to_real y_f)
1080       by
1081         (-. x_rel *. x_factor -. x_abs) /. to_real y_f
1082         <=. (x_rel *. x_factor +. x_abs) /. abs (to_real y_f)
1083     };
1084     div_order_compat2 (x -. x_rel *. x_factor -. x_abs) (to_real x_f) (to_real y_f);
1085   in ()
1086 end