ieee_float.mlw 32.5 KB
Newer Older
1 2
(** {1 Formalization of Floating-Point Arithmetic}

3
  Full float theory (with infinities and NaN).
4 5 6 7 8

  A note on intended semantics: we use the same idea as the SMTLIB floating
  point theory, that defers any inconsistencies to the "parent" document.
  Hence, in doubt, the correct axiomatisation is one that implements
  [BTRW14] "An Automatable Formal Semantics for IEEE-754 Floating-Point
9
  Arithmetic", which in turn defers any inconsistencies to IEEE-754.
10 11 12 13 14 15 16 17 18

  This theory is split into two parts: the first part talks about IEEE
  operations and this is what you should use as a user, the second part is
  internal only and is an axiomatisation for the provers that do not
  natively support floating point. You should not use any symbols you find
  there in your verification conditions as solvers with native floating
  point support will leave them uninterpreted.
*)

19 20
(** {2 Rounding Modes} *)

21
module RoundingMode
22
  type mode = RNE | RNA | RTP | RTN | RTZ
23 24 25 26 27 28
  (** {h <ul>
     <li>RNE : Round Nearest ties to Even
     <li>RNA : Round Nearest ties to Away
     <li>RTP : Round Towards Positive
     <li>RTN : Round Towards Negative
     <li>RTZ : Round Towards Zero
Guillaume Melquiond's avatar
Guillaume Melquiond committed
29
     </ul>} *)
30 31

  predicate to_nearest (m:mode) = m = RNE \/ m = RNA
32 33
end

34
module GenericFloat
35

36 37
  use int.Int
  use bv.Pow2int
38 39 40
  use real.Abs as Abs
  use real.FromInt as FromInt
  use real.Truncate as Truncate
41
  use real.RealInfix
42 43
  use export RoundingMode

44
  (** {2 Part I - Public Interface}   *)
45 46

  constant eb : int
47 48
  (** the number of bits in the exponent.  *)

49
  constant sb : int
50
  (** the number of bits in the significand, including the hidden bit. *)
51 52 53 54

  axiom eb_gt_1: 1 < eb
  axiom sb_gt_1: 1 < sb

55
  (** {3 Sorts} *)
56 57

  type t
MARCHE Claude's avatar
MARCHE Claude committed
58 59
  (** abstract type to denote floating-point numbers, including the special values
      for infinities and NaNs *)
60

61
  (** {3 Constructors and Constants} *)
62

63
  val constant zeroF     : t   (** +0.0 *)
64 65 66
  (* exp_bias = 2^(sb - 1) - 1 *)
  (* max_finite_exp = 2^sb - 2 - exp_bias = exp_bias *)
  (* max_significand = (2^eb + 2^eb - 1) * 2^(1-eb) *)
67
  (* max_value = (2^eb + 2^eb - 1) * 2^(1-eb) * 2 ^ max_finite_exp *)
68
  (* [m;x] = ( 1 + m * 2^(1-eb) ) * 2^( x - exp_bias ) *)
69 70
  (* max_value = [2^(eb-1); 2^sb - 2] *)

71
  (** {3 Operators} *)
72

73 74 75 76
  val function add mode t t : t
  val function sub mode t t : t
  val function mul mode t t : t
  val function div mode t t : t
MARCHE Claude's avatar
MARCHE Claude committed
77
    (** The four basic operations, rounded in the given mode *)
78

79 80 81 82
  val function abs         t   : t   (** Absolute value *)
  val function neg         t   : t   (** Opposite *)
  val function fma  mode t t t : t   (** Fused multiply-add: x * y + z *)
  val function sqrt mode   t   : t   (** Square root *)
83

84 85 86 87 88
  let function (.-_) (x:t)   : t = neg x
  let function (.+)  (x y:t) : t = add RNE x y
  let function (.-)  (x y:t) : t = sub RNE x y
  let function (.*)  (x y:t) : t = mul RNE x y
  let function (./)  (x y:t) : t = div RNE x y
MARCHE Claude's avatar
MARCHE Claude committed
89
    (** Notations for operations in the default mode RNE *)
90

91
  val function roundToIntegral mode t : t
MARCHE Claude's avatar
MARCHE Claude committed
92
    (** Rounding to an integer *)
93 94 95

  function min       t t : t
  function max       t t : t
MARCHE Claude's avatar
MARCHE Claude committed
96 97 98 99
  (** Minimum and Maximum
      Note that we have to follow IEEE-754 and SMTLIB here. Two things to
     note in particular:
     1) min(-0, 0) is either 0 or -0, there is a choice
100

MARCHE Claude's avatar
MARCHE Claude committed
101
     2) if either argument is NaN then the other argument is returned
102

103 104 105 106
     Due to the unclear status of min and max in IEEE norm, we
     intentionally not provide these function as "val"s to be used in
     programs

MARCHE Claude's avatar
MARCHE Claude committed
107
   *)
108

MARCHE Claude's avatar
MARCHE Claude committed
109
  (** {3 Comparisons} *)
110

111 112 113 114 115
  val predicate le t t
  val predicate lt t t
  let predicate ge (x:t) (y:t) = le y x
  let predicate gt (x:t) (y:t) = lt y x
  val predicate eq t t
MARCHE Claude's avatar
MARCHE Claude committed
116
  (** equality on floats, different from = since not (eq NaN NaN) *)
117

118 119 120 121 122
  let predicate (.<=) (x:t) (y:t) = le x y
  let predicate (.<)  (x:t) (y:t) = lt x y
  let predicate (.>=) (x:t) (y:t) = ge x y
  let predicate (.>)  (x:t) (y:t) = gt x y
  let predicate (.=)  (x:t) (y:t) = eq x y
MARCHE Claude's avatar
MARCHE Claude committed
123 124
  (** Notations *)

125

126
  (** {3 Classification of numbers} *)
127 128 129

  predicate is_normal    t
  predicate is_subnormal t
130
  predicate is_zero      t
131 132 133 134 135
  predicate is_infinite  t
  predicate is_nan       t
  predicate is_positive  t
  predicate is_negative  t

136
  (** helper predicate for zeros, normals and subnormals. not defined
137 138 139 140 141 142 143 144 145 146
     so that the axiomatisation below can use it without talking about
     subnormals *)
  predicate is_finite    t

  predicate is_plus_infinity  (x:t) = is_infinite x /\ is_positive x
  predicate is_minus_infinity (x:t) = is_infinite x /\ is_negative x
  predicate is_plus_zero      (x:t) = is_zero x /\ is_positive x
  predicate is_minus_zero     (x:t) = is_zero x /\ is_negative x
  predicate is_not_nan        (x:t) = is_finite x \/ is_infinite x

147 148 149 150
  axiom is_not_nan: forall x:t. is_not_nan x <-> not (is_nan x)

  axiom is_not_finite: forall x:t.
    not (is_finite x) <-> (is_infinite x \/ is_nan x)
151

MARCHE Claude's avatar
MARCHE Claude committed
152 153


154
  (** {3 Conversions from other sorts} *)
155 156 157 158 159 160 161 162 163

  (* from bitvec binary interchange                           *)
  (*   partly done with from_binary (for literals only)       *)
  (* from another fp                 - see FloatConverter     *)
  (* from real                                                *)
  (*   partly done with (!) (for literals only)               *)
  (* from unsigned integer bitvector - see Float_BV_Converter *)
  (* from signed integer bitvector                            *)

164
  (** {3 Conversions to other sorts} *)
165 166 167 168 169 170

  (* to unsigned integer bitvector   - see Float_BV_Converter *)
  (* to signed integer bitvector                              *)
  (* to real *)
  function to_real   t    : real

171
  (** {2 Part II - Private Axiomatisation} *)
172

173
  (** {3 Constructors and Constants} *)
174

175
  axiom zeroF_is_positive : is_positive zeroF
176
  axiom zeroF_is_zero     : is_zero zeroF
177

178 179
  axiom zero_to_real : forall x [is_zero x].
    is_zero x <-> is_finite x /\ to_real x = 0.0
180

181
  (** {3 Conversions from other sorts} *)
182 183 184 185

  (* with mathematical int *)
  (* note that these conversions do not feature in SMTLIB *)

186 187
  (* intended semantics for of_int are the same as (_ to_fp eb sb) with a   *)
  (* suitably sized bitvector, large enough to hold x                       *)
188 189 190
  (* note values >= than the below should result in infinities              *)
  (*    float32 : 0x1ffffff * 2^103                                         *)
  (*    float64 : 0x3fffffffffffff * 2^970                                  *)
191
  (* also note that this function never yields a subnormal, or a NaN, or -0 *)
192 193
  function of_int (m:mode) (x:int) : t

194
  (** {3 Conversions to other sorts} *)
195

196 197
  (* Intended semantics for to_int are the same as (_ fp.to_sbv) with a    *)
  (* suitably sized bitvector. Safe minimum sizes are given below:         *)
198 199 200 201 202 203 204 205 206 207 208 209
  (*    float32 : 129                                                      *)
  (*    float64 : 1025                                                     *)
  (* In particular this function should be uninterpreted for infinities    *)
  (* and NaN. Take care that no conclusion can be made on the result based *)
  (* on the size of the bitvector chosen in those cases, i.e. this should  *)
  (* not hold:                                                             *)
  (*    to_int +INF < 2 ** 2048     // nope                                *)
  (*    to_int +INF > 0             // nope                                *)
  function to_int (m:mode) (x:t) : int

  axiom zero_of_int : forall m. zeroF = of_int m 0

210
  (** {3 Arithmetic} *)
211

212 213
  (* The intended meaning for round is the rounding for floating point as  *)
  (* described on p17 of IEEE-754. For results where the corresponding     *)
214 215 216 217 218 219 220 221
  (* floating point result would be infinity or NaN this function should   *)
  (* be uninterpreted.                                                     *)
  (*                                                                       *)
  (* Note that this means round (+INF) > 0 is not true.                    *)
  (* Note also this means round (2*INF) > round (INF) is not true either.  *)
  function round mode real : real

  constant max_real : real (* defined when cloning *)
222 223
  constant max_int  : int

224
  constant emax : int = pow2 (eb - 1)
225

226
  axiom max_int : max_int = pow2 emax - pow2 (emax - sb)
227
  axiom max_real_int: max_real = FromInt.from_int max_int
228 229 230

  predicate in_range (x:real) = -. max_real <=. x <=. max_real

231
  predicate in_int_range (i:int) = - max_int <= i <= max_int
232

233
  axiom is_finite: forall x:t. is_finite x -> in_range (to_real x)
234 235 236 237 238 239

  (* used as a condition to propagate is_finite *)
  predicate no_overflow (m:mode) (x:real) = in_range (round m x)

  (* Axioms on round *)

240
  axiom Bounded_real_no_overflow [@W:non_conservative_extension:N] :
241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261
    forall m:mode, x:real. in_range x -> no_overflow m x

  axiom Round_monotonic :
    forall m:mode, x y:real. x <=. y -> round m x <=. round m y

  axiom Round_idempotent :
    forall m1 m2:mode, x:real. round m1 (round m2 x) = round m2 x

  axiom Round_to_real :
    forall m:mode, x:t. is_finite x -> round m (to_real x) = to_real x

  (** rounding up and down *)
  axiom Round_down_le:
    forall x:real. round RTN x <=. x
  axiom Round_up_ge:
    forall x:real. round RTP x >=. x
  axiom Round_down_neg:
    forall x:real. round RTN (-.x) = -. round RTP x
  axiom Round_up_neg:
    forall x:real. round RTP (-.x) = -. round RTN x

262 263 264
  (* The biggest representable integer whose predecessor (i.e. -1) is  representable *)
  constant pow2sb : int (* defined when cloning *)
  axiom pow2sb: pow2sb = pow2 sb
265

266 267
  (* range in which every integer is representable *)
  predicate in_safe_int_range (i: int) = - pow2sb <= i <= pow2sb
268 269 270 271 272

  (* round and integers *)

  axiom Exact_rounding_for_integers:
    forall m:mode, i:int.
273
      in_safe_int_range i ->
274 275 276 277 278 279 280 281 282 283 284 285
        round m (FromInt.from_int i) = FromInt.from_int i

  (** {3 Comparisons} *)

  (** Comparison predicates *)

  predicate same_sign (x y : t) =
    (is_positive x /\ is_positive y) \/ (is_negative x /\ is_negative y)
  predicate diff_sign (x y : t) =
    (is_positive x /\ is_negative y) \/ (is_negative x /\ is_positive y)

  axiom feq_eq: forall x y.
286
    is_finite x -> is_finite y -> not (is_zero x) -> x .= y -> x = y
287 288 289 290 291 292

  axiom eq_feq: forall x y.
    is_finite x -> is_finite y -> x = y -> x .= y

  axiom eq_refl: forall x. is_finite x -> x .= x

293
  axiom eq_sym :
MARCHE Claude's avatar
MARCHE Claude committed
294
    forall x y. x .= y -> y .= x
295

296
  axiom eq_trans :
MARCHE Claude's avatar
MARCHE Claude committed
297
    forall x y z. x .= y -> y .= z -> x .= z
298 299 300 301 302 303 304 305

  axiom eq_zero: zeroF .= (.- zeroF)

  axiom eq_to_real_finite: forall x y.
    is_finite x /\ is_finite y -> (x .= y <-> to_real x = to_real y)

  axiom eq_special: forall x y. x .= y ->
       (is_not_nan x  /\ is_not_nan y
306 307
    /\ ((is_finite x /\ is_finite y)
        \/ (is_infinite x /\ is_infinite y /\ same_sign x y)))
308

309 310
  axiom lt_finite: forall x y [lt x y].
    is_finite x /\ is_finite y -> (lt x y <-> to_real x <. to_real y)
311

312 313
  axiom le_finite: forall x y [le x y].
    is_finite x /\ is_finite y -> (le x y <-> to_real x <=. to_real y)
314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329

  lemma le_lt_trans:
    forall x y z:t. x .<= y /\ y .< z -> x .< z

  lemma lt_le_trans:
    forall x y z:t. x .< y /\ y .<= z -> x .< z

  lemma le_ge_asym:
    forall x y:t. x .<= y /\ x .>= y -> x .= y

  lemma not_lt_ge: forall x y:t.
    not (x .< y) /\ is_not_nan x /\ is_not_nan y -> x .>= y

  lemma not_gt_le: forall x y:t.
    not (x .> y) /\ is_not_nan x /\ is_not_nan y -> x .<= y

330
  axiom le_special: forall x y [le x y]. le x y ->
331 332 333 334
      ((is_finite x /\ is_finite y)
   \/ ((is_minus_infinity x /\ is_not_nan y)
   \/  (is_not_nan x /\ is_plus_infinity y)))

335
  axiom lt_special: forall x y [lt x y]. lt x y ->
336 337 338 339
      ((is_finite x /\ is_finite y)
   \/ ((is_minus_infinity x /\ is_not_nan y /\ not (is_minus_infinity y))
   \/  (is_not_nan x /\ not (is_plus_infinity x) /\ is_plus_infinity y)))

340 341
  axiom lt_lt_finite: forall x y z. lt x y -> lt y z -> is_finite y

342 343 344 345 346 347 348 349 350 351 352
  (*  lemmas on sign *)
  axiom positive_to_real: forall x[is_positive x|to_real x >=. 0.0].
    is_finite x -> is_positive x -> to_real x >=. 0.0
  axiom to_real_positive: forall x[is_positive x].
    is_finite x -> to_real x >. 0.0 -> is_positive x

  axiom negative_to_real: forall x [is_negative x|to_real x <=. 0.0].
    is_finite x -> is_negative x -> to_real x <=. 0.0
  axiom to_real_negative: forall x [is_negative x].
    is_finite x -> to_real x <. 0.0 -> is_negative x

353 354 355 356
  axiom negative_xor_positive: forall x.
    not (is_positive x /\ is_negative x)
  axiom negative_or_positive: forall x.
    is_not_nan x -> is_positive x \/ is_negative x
357 358 359 360 361 362 363 364 365 366 367 368 369 370

  lemma diff_sign_trans:
    forall x y z:t. (diff_sign x y /\ diff_sign y z) -> same_sign x z

  lemma diff_sign_product:
    forall x y:t.
      (is_finite x /\ is_finite y /\ to_real x *. to_real y <. 0.0) ->
        diff_sign x y

  lemma same_sign_product:
    forall x y:t.
      (is_finite x /\ is_finite y /\ same_sign x y) ->
        to_real x *. to_real y >=. 0.0

371 372 373
  predicate product_sign (z x y: t) =
    (same_sign x y -> is_positive z) /\ (diff_sign x y -> is_negative z)

374
  (** {3 Overflow} *)
375

376 377 378
  (* This predicate is used to tell what is the result of a rounding
     in case of overflow in the axioms specifying add/sub/mul and fma
     *)
379 380 381 382 383 384 385 386 387 388 389
  predicate overflow_value (m:mode) (x:t) =
    match m with
    | RTN -> if is_positive x then is_finite x /\ to_real x = max_real
                              else is_infinite x
    | RTP -> if is_positive x then is_infinite x
                              else is_finite x /\ to_real x = -. max_real
    | RTZ -> if is_positive x then is_finite x /\ to_real x = max_real
                              else is_finite x /\ to_real x = -. max_real
    | (RNA | RNE) -> is_infinite x
    end

390 391
  (* This predicate is used to tell what is the sign of zero in the
     axioms specifying add and sub *)
392 393 394 395 396 397 398 399 400 401
  predicate sign_zero_result (m:mode) (x:t) =
    is_zero x ->
    match m with
    | RTN -> is_negative x
    | _   -> is_positive x
    end

  (** {3 binary operations} *)

  axiom add_finite: forall m:mode, x y:t [add m x y].
402
    is_finite x -> is_finite y -> no_overflow m (to_real x +. to_real y) ->
403
    is_finite (add m x y) /\
404
    to_real (add m x y) = round m (to_real x +. to_real y)
405

406 407
  lemma add_finite_rev: forall m:mode, x y:t [add m x y].
    is_finite (add m x y) ->
408
    is_finite x /\ is_finite y
409

410 411 412 413 414 415
  lemma add_finite_rev_n: forall m:mode, x y:t [add m x y].
    to_nearest m ->
    is_finite (add m x y) ->
    no_overflow m (to_real x +. to_real y) /\
    to_real (add m x y) = round m (to_real x +. to_real y)

416 417
  axiom sub_finite: forall m:mode, x y:t [sub m x y].
    is_finite x -> is_finite y -> no_overflow m (to_real x -. to_real y) ->
418
    is_finite (sub m x y) /\
419
    to_real (sub m x y) = round m (to_real x -. to_real y)
420

421 422
  lemma sub_finite_rev: forall m:mode, x y:t [sub m x y].
    is_finite (sub m x y) ->
423
    is_finite x /\ is_finite y
424

425 426 427 428 429 430
  lemma sub_finite_rev_n: forall m:mode, x y:t [sub m x y].
    to_nearest m ->
    is_finite (sub m x y) ->
    no_overflow m (to_real x -. to_real y) /\
    to_real (sub m x y) = round m (to_real x -. to_real y)

431 432
  axiom mul_finite: forall m:mode, x y:t [mul m x y].
    is_finite x -> is_finite y -> no_overflow m (to_real x *. to_real y) ->
433 434 435
    is_finite (mul m x y) /\
    to_real (mul m x y) = round m (to_real x *. to_real y)

436 437
  lemma mul_finite_rev: forall m:mode, x y:t [mul m x y].
    is_finite (mul m x y) ->
438
    is_finite x /\ is_finite y
439

440 441 442 443 444 445
  lemma mul_finite_rev_n: forall m:mode, x y:t [mul m x y].
    to_nearest m ->
    is_finite (mul m x y) ->
    no_overflow m (to_real x *. to_real y) /\
    to_real (mul m x y) = round m (to_real x *. to_real y)

446 447
  axiom div_finite: forall m:mode, x y:t [div m x y].
    is_finite x -> is_finite y ->
448
    not is_zero y -> no_overflow m (to_real x /. to_real y) ->
449 450 451
    is_finite (div m x y) /\
    to_real (div m x y) = round m (to_real x /. to_real y)

452 453
  lemma div_finite_rev: forall m:mode, x y:t [div m x y].
    is_finite (div m x y) ->
454
    (is_finite x /\ is_finite y /\ not is_zero y) \/
455 456
    (is_finite x /\ is_infinite y /\ to_real (div m x y) = 0.)

457 458 459 460 461 462
  lemma div_finite_rev_n: forall m:mode, x y:t [div m x y].
    to_nearest m ->
    is_finite (div m x y) -> is_finite y ->
    no_overflow m (to_real x /. to_real y) /\
    to_real (div m x y) = round m (to_real x /. to_real y)

463 464 465
  axiom neg_finite: forall x:t [neg x].
    is_finite x ->
    is_finite (neg x) /\
466 467
    to_real (neg x) = -. to_real x

468 469 470 471 472
  lemma neg_finite_rev: forall x:t [neg x].
    is_finite (neg x) ->
    is_finite x /\
    to_real (neg x) = -. to_real x

473 474
  axiom abs_finite: forall x:t [abs x].
    is_finite x ->
475 476 477 478
    is_finite (abs x) /\
    to_real (abs x) = Abs.abs (to_real x) /\
    is_positive (abs x)

479 480 481 482 483
  lemma abs_finite_rev: forall x:t [abs x].
    is_finite (abs x) ->
    is_finite x /\
    to_real (abs x) = Abs.abs (to_real x)

484 485
  axiom abs_universal : forall x:t [abs x]. not (is_negative (abs x))

486 487 488 489 490
  axiom fma_finite: forall m:mode, x y z:t [fma m x y z].
    is_finite x -> is_finite y -> is_finite z ->
    no_overflow m (to_real x *. to_real y +. to_real z) ->
    is_finite (fma m x y z) /\
    to_real (fma m x y z) = round m (to_real x *. to_real y +. to_real z)
491

492 493
  lemma fma_finite_rev: forall m:mode, x y z:t [fma m x y z].
    is_finite (fma m x y z) ->
494
    is_finite x /\ is_finite y /\ is_finite z
495

496 497 498 499 500 501
  lemma fma_finite_rev_n: forall m:mode, x y z:t [fma m x y z].
    to_nearest m ->
    is_finite (fma m x y z) ->
    no_overflow m (to_real x *. to_real y +. to_real z) /\
    to_real (fma m x y z) = round m (to_real x *. to_real y +. to_real z)

502
  use real.Square as S
503

504
  axiom sqrt_finite: forall m:mode, x:t [sqrt m x].
505
    is_finite x -> to_real x >=. 0. ->
506
    is_finite (sqrt m x) /\
507
    to_real (sqrt m x) = round m (S.sqrt (to_real x))
508

509 510 511
  lemma sqrt_finite_rev: forall m:mode, x:t [sqrt m x].
    is_finite (sqrt m x) ->
    is_finite x /\ to_real x >=. 0. /\
512
    to_real (sqrt m x) = round m (S.sqrt (to_real x))
513

514 515 516
  predicate same_sign_real (x:t) (r:real) =
    (is_positive x /\ r >. 0.0) \/ (is_negative x /\ r <. 0.0)

517
  axiom add_special: forall m:mode, x y:t [add m x y].
518 519 520 521 522 523 524 525 526 527 528 529 530 531
    let r = add m x y in
    (is_nan x \/ is_nan y -> is_nan r)
    /\
    (is_finite x /\ is_infinite y -> is_infinite r /\ same_sign r y)
    /\
    (is_infinite x /\ is_finite y -> is_infinite r /\ same_sign r x)
    /\
    (is_infinite x /\ is_infinite y /\ same_sign x y
      -> is_infinite r /\ same_sign r x)
    /\
    (is_infinite x /\ is_infinite y /\ diff_sign x y -> is_nan r)
    /\
    (is_finite x /\ is_finite y /\ not no_overflow m (to_real x +. to_real y)
      -> same_sign_real r (to_real x +. to_real y) /\ overflow_value m r)
532 533 534
    /\
    (is_finite x /\ is_finite y
      -> if same_sign x y then same_sign r x else sign_zero_result m r)
535

536
  axiom sub_special: forall m:mode, x y:t [sub m x y].
537 538 539 540 541 542 543 544 545 546 547 548 549
    let r = sub m x y in
    (is_nan x \/ is_nan y -> is_nan r)
    /\
    (is_finite x /\ is_infinite y -> is_infinite r /\ diff_sign r y)
    /\
    (is_infinite x /\ is_finite y -> is_infinite r /\ same_sign r x)
    /\
    (is_infinite x /\ is_infinite y /\ same_sign x y -> is_nan r)
    /\
    (is_infinite x /\ is_infinite y /\ diff_sign x y
      -> is_infinite r /\ same_sign r x)
    /\
    (is_finite x /\ is_finite y /\ not no_overflow m (to_real x -. to_real y)
550
      -> same_sign_real r (to_real x -. to_real y) /\ overflow_value m r)
551 552 553
    /\
    (is_finite x /\ is_finite y
      -> if diff_sign x y then same_sign r x else sign_zero_result m r)
554

555
  axiom mul_special: forall m:mode, x y:t [mul m x y].
556 557 558 559 560 561 562 563 564 565 566
    let r = mul m x y in
       (is_nan x \/ is_nan y -> is_nan r)
    /\ (is_zero x /\ is_infinite y -> is_nan r)
    /\ (is_finite x /\ is_infinite y /\ not (is_zero x)
         -> is_infinite r)
    /\ (is_infinite x /\ is_zero y -> is_nan r)
    /\ (is_infinite x /\ is_finite y /\ not (is_zero  y)
         -> is_infinite r)
    /\ (is_infinite x /\ is_infinite y -> is_infinite r)
    /\ (is_finite x /\ is_finite y /\ not no_overflow m (to_real x *. to_real y)
         -> overflow_value m r)
567
    /\ (not is_nan r -> product_sign r x y)
568

569
  axiom div_special: forall m:mode, x y:t [div m x y].
570 571 572 573 574 575
    let r = div m x y in
       (is_nan x \/ is_nan y -> is_nan r)
    /\ (is_finite x /\ is_infinite y -> is_zero r)
    /\ (is_infinite x /\ is_finite y -> is_infinite r)
    /\ (is_infinite x /\ is_infinite y -> is_nan r)
    /\ (is_finite x /\ is_finite y /\ not (is_zero y) /\
576 577
        not no_overflow m (to_real x /. to_real y)
         -> overflow_value m r)
578
    /\ (is_finite x /\ is_zero y /\ not (is_zero x)
579
         -> is_infinite r)
580
    /\ (is_zero x /\ is_zero y -> is_nan r)
581 582 583 584 585 586
    /\ (not is_nan r -> product_sign r x y)

  axiom neg_special: forall x:t [neg x].
       (is_nan x -> is_nan (neg x))
    /\ (is_infinite x -> is_infinite (neg x))
    /\ (not is_nan x -> diff_sign x (neg x))
587

588 589 590 591
  axiom abs_special: forall x:t [abs x].
       (is_nan x -> is_nan (abs x))
    /\ (is_infinite x -> is_infinite (abs x))
    /\ (not is_nan x -> is_positive (abs x))
592

593 594
  axiom fma_special: forall m:mode, x y z:t [fma m x y z].
    let r = fma m x y z in
595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618
       (is_nan x \/ is_nan y \/ is_nan z -> is_nan r)
    /\ (is_zero x /\ is_infinite y -> is_nan r)
    /\ (is_infinite x /\ is_zero y -> is_nan r)
    /\ (is_finite x /\ not (is_zero x) /\ is_infinite y /\ is_finite z
        -> is_infinite r /\ product_sign r x y)
    /\ (is_finite x /\ not (is_zero x) /\ is_infinite y /\ is_infinite z
        -> (if product_sign z x y then is_infinite r /\ same_sign r z
            else is_nan r))
    /\ (is_infinite x /\ is_finite y /\ not (is_zero y) /\ is_finite z
        -> is_infinite r /\ product_sign r x y)
    /\ (is_infinite x /\ is_finite y /\ not (is_zero y) /\ is_infinite z
        -> (if product_sign z x y then is_infinite r /\ same_sign r z
            else is_nan r))
    /\ (is_infinite x /\ is_infinite y /\ is_finite z
        -> is_infinite r /\ product_sign r x y)
    /\ (is_finite x /\ is_finite y /\ is_infinite z
        -> is_infinite r /\ same_sign r z)
    /\ (is_infinite x /\ is_infinite y /\ is_infinite z
        -> (if product_sign z x y then is_infinite r /\ same_sign r z
            else is_nan r))
    /\ (is_finite x /\ is_finite y /\ is_finite z /\
        not no_overflow m (to_real x *. to_real y +. to_real z)
        -> same_sign_real r (to_real x *. to_real y +. to_real z)
        /\ overflow_value m r)
619 620 621 622 623 624 625 626 627
    /\ (is_finite x /\ is_finite y /\ is_finite z
        -> if product_sign z x y then same_sign r z
           else (to_real x *. to_real y +. to_real z = 0.0 ->
                 if m = RTN then is_negative r else is_positive r))

  axiom sqrt_special: forall m:mode, x:t [sqrt m x].
    let r = sqrt m x in
       (is_nan x -> is_nan r)
    /\ (is_plus_infinity x -> is_plus_infinity r)
628 629
    /\ (is_minus_infinity x -> is_nan r)
    /\ (is_finite x /\ to_real x <. 0.0 -> is_nan r)
630 631
    /\ (is_zero x -> same_sign r x)
    /\ (is_finite x /\ to_real x >. 0.0 -> is_positive r)
632

633 634 635
  (* exact arithmetic with integers *)

  axiom of_int_add_exact: forall m n, i j.
636 637
    in_safe_int_range i -> in_safe_int_range j ->
    in_safe_int_range (i + j) -> eq (of_int m (i + j)) (add n (of_int m i) (of_int m j))
638 639

  axiom of_int_sub_exact: forall m n, i j.
640 641
    in_safe_int_range i -> in_safe_int_range j ->
    in_safe_int_range (i - j) -> eq (of_int m (i - j)) (sub n (of_int m i) (of_int m j))
642 643

  axiom of_int_mul_exact: forall m n, i j.
644 645
    in_safe_int_range i -> in_safe_int_range j ->
    in_safe_int_range (i * j) -> eq (of_int m (i * j)) (mul n (of_int m i) (of_int m j))
646

647
  (* min and max *)
648 649 650 651 652 653 654 655

  lemma Min_r : forall x y:t. y .<= x -> (min x y) .= y
  lemma Min_l : forall x y:t. x .<= y -> (min x y) .= x
  lemma Max_r : forall x y:t. y .<= x -> (max x y) .= x
  lemma Max_l : forall x y:t. x .<= y -> (max x y) .= y

  (* _____________ *)

656
  use real.Truncate as Truncate
657

658
  (* This predicate specify if a float is finite and is an integer *)
659 660
  predicate is_int (x:t)

661 662 663
  (** characterizing integers *)

  (* by construction *)
664 665
  axiom zeroF_is_int: is_int zeroF

666
  axiom of_int_is_int: forall m, x.
667
    in_int_range x -> is_int (of_int m x)
668

669 670
  axiom big_float_is_int: forall m i.
    is_finite i ->
671
    i .<= neg (of_int m pow2sb) \/ (of_int m pow2sb) .<= i ->
672
      is_int i
673 674 675 676

  axiom roundToIntegral_is_int: forall m:mode, x:t. is_finite x ->
    is_int (roundToIntegral m x)

677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712
  (* by propagation *)
  axiom eq_is_int: forall x y. eq x y -> is_int x -> is_int y

  axiom add_int: forall x y m. is_int x -> is_int y ->
    is_finite (add m x y) -> is_int (add m x y)

  axiom sub_int: forall x y m. is_int x -> is_int y ->
    is_finite (sub m x y) -> is_int (sub m x y)

  axiom mul_int: forall x y m. is_int x -> is_int y ->
    is_finite (mul m x y) -> is_int (mul m x y)

  axiom fma_int: forall x y z m. is_int x -> is_int y -> is_int z ->
    is_finite (fma m x y z) -> is_int (fma m x y z)

  axiom neg_int: forall x. is_int x -> is_int (neg x)

  axiom abs_int: forall x. is_int x -> is_int (abs x)

  (** basic properties of float integers *)

  axiom is_int_of_int: forall x m m'.
    is_int x -> eq x (of_int m' (to_int m x))

  axiom is_int_to_int: forall m x.
    is_int x -> in_int_range (to_int m x)

  axiom is_int_is_finite: forall x.
    is_int x -> is_finite x

  axiom int_to_real: forall m x.
    is_int x -> to_real x = FromInt.from_int (to_int m x)

(*  axiom int_mode: forall m1 m2 x.
    is_int x -> to_int m1 x = to_int m2 x  etc ...*)

713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731
  (** rounding ints *)

  axiom truncate_int: forall m:mode, i:t. is_int i ->
    roundToIntegral m i .= i

  (** truncate *)

  axiom truncate_neg: forall x:t.
    is_finite x -> is_negative x -> roundToIntegral RTZ x = roundToIntegral RTP x

  axiom truncate_pos: forall x:t.
    is_finite x -> is_positive x -> roundToIntegral RTZ x = roundToIntegral RTN x

  (** ceil *)

  axiom ceil_le: forall x:t. is_finite x -> x .<= (roundToIntegral RTP x)

  axiom ceil_lest: forall x y:t. x .<= y /\ is_int y -> (roundToIntegral RTP x) .<= y

732 733
  axiom ceil_to_real: forall x:t.
    is_finite x ->
734 735 736
      to_real (roundToIntegral RTP x) = FromInt.from_int (Truncate.ceil (to_real x))

  axiom ceil_to_int: forall m:mode, x:t.
737
    is_finite x ->
738 739 740 741 742 743 744 745
      to_int m (roundToIntegral RTP x) = Truncate.ceil (to_real x)

  (** floor *)

  axiom floor_le: forall x:t. is_finite x -> (roundToIntegral RTN x) .<= x

  axiom floor_lest: forall x y:t. y .<= x /\ is_int y -> y .<= (roundToIntegral RTN x)

746 747
  axiom floor_to_real: forall x:t.
    is_finite x ->
748 749 750
      to_real (roundToIntegral RTN x) = FromInt.from_int (Truncate.floor (to_real x))

  axiom floor_to_int: forall m:mode, x:t.
751
    is_finite x ->
752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776
      to_int m (roundToIntegral RTN x) = Truncate.floor (to_real x)

  (* Rna *)

  axiom RNA_down:
    forall x:t. (x .- (roundToIntegral RTN x)) .< ((roundToIntegral RTP x) .- x) ->
      roundToIntegral RNA x = roundToIntegral RTN x

  axiom RNA_up:
    forall x:t. (x .- (roundToIntegral RTN x)) .> ((roundToIntegral RTP x) .- x) ->
      roundToIntegral RNA x = roundToIntegral RTP x

  axiom RNA_down_tie:
    forall x:t. (x .- (roundToIntegral RTN x)) .= ((roundToIntegral RTP x) .- x) ->
      is_negative x -> roundToIntegral RNA x = roundToIntegral RTN x

  axiom RNA_up_tie:
    forall x:t. ((roundToIntegral RTP x) .- x) .= (x .- (roundToIntegral RTN x)) ->
    is_positive x -> roundToIntegral RNA x = roundToIntegral RTP x

  (* to_int *)
  axiom to_int_roundToIntegral: forall m:mode, x:t.
    to_int m x = to_int m (roundToIntegral m x)

  axiom to_int_monotonic: forall m:mode, x y:t.
777
    is_finite x -> is_finite y -> le x y -> to_int m x <= to_int m y
778 779

  axiom to_int_of_int: forall m:mode, i:int.
780
    in_safe_int_range i ->
781 782 783 784 785
      to_int m (of_int m i) = i

  axiom eq_to_int: forall m, x y. is_finite x -> x .= y ->
    to_int m x = to_int m y

786 787 788
  axiom neg_to_int: forall m x.
    is_int x -> to_int m (neg x) = - (to_int m x)

789 790 791 792
  axiom roundToIntegral_is_finite  : forall m:mode, x:t. is_finite x ->
    is_finite (roundToIntegral m x)
end

793 794
(** {2 Conversions to/from bitvectors} *)

795
module Float_BV_Converter
796 797 798 799
  use bv.BV8 as BV8
  use bv.BV16 as BV16
  use bv.BV32 as BV32
  use bv.BV64 as BV64
800
  use RoundingMode
801 802 803 804

  (* with unsigned int as bitvector *)
  type t                        (* float *)

805 806 807 808
  val function of_ubv8  mode BV8.t  : t
  val function of_ubv16 mode BV16.t : t
  val function of_ubv32 mode BV32.t : t
  val function of_ubv64 mode BV64.t : t
809

810 811 812 813
  val function to_ubv8  mode t : BV8.t
  val function to_ubv16 mode t : BV16.t
  val function to_ubv32 mode t : BV32.t
  val function to_ubv64 mode t : BV64.t
814

815
  use real.RealInfix
816
  use real.FromInt as FromInt
817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840

  predicate is_finite t
  predicate le t t
  function to_real         t : real
  function round   mode real : real

  (** of unsigned bv axioms  *)
  (* only true for big enough floats...  *)

  axiom of_ubv8_is_finite : forall m, x. is_finite (of_ubv8  m x)
  axiom of_ubv16_is_finite: forall m, x. is_finite (of_ubv16 m x)
  axiom of_ubv32_is_finite: forall m, x. is_finite (of_ubv32 m x)
  axiom of_ubv64_is_finite: forall m, x. is_finite (of_ubv64 m x)

  axiom of_ubv8_monotonic :
    forall m, x y. BV8.ule  x y -> le (of_ubv8  m x) (of_ubv8  m y)
  axiom of_ubv16_monotonic:
    forall m, x y. BV16.ule x y -> le (of_ubv16 m x) (of_ubv16 m y)
  axiom of_ubv32_monotonic:
    forall m, x y. BV32.ule x y -> le (of_ubv32 m x) (of_ubv32 m y)
  axiom of_ubv64_monotonic:
    forall m, x y. BV64.ule x y -> le (of_ubv64 m x) (of_ubv64 m y)

  axiom of_ubv8_to_real : forall m, x.
841
    to_real (of_ubv8  m x) = FromInt.from_int (BV8.t'int x)
842
  axiom of_ubv16_to_real: forall m, x.
843
    to_real (of_ubv16 m x) = FromInt.from_int (BV16.t'int x)
844 845
  (* of_ubv32_to_real is defined at cloning *)
  axiom of_ubv64_to_real: forall m, x.
846
    to_real (of_ubv64 m x) = round m (FromInt.from_int (BV64.t'int x))
847 848
end

849 850
(** {2 Standard simple precision floats (32 bits)} *)

851
module Float32
852
  use int.Int as Int
853
  use real.Real
854

855
  type t = < float 8 24 >
856

857
  constant pow2sb : int = 16777216
858 859 860 861
  constant max_real : real = 0x1.FFFFFEp127

  clone export GenericFloat with
    type t = t,
862 863
    constant eb = t'eb,
    constant sb = t'sb,
864
    constant max_real = max_real,
865
    constant pow2sb = pow2sb,
866 867
    function to_real = t'real,
    predicate is_finite = t'isFinite,
868 869 870
    goal eb_gt_1,
    goal sb_gt_1,
    goal max_int,
871 872
    goal pow2sb,
    axiom . (* TODO: "lemma"? "goal"? *)
873

874 875
  lemma round_bound_ne :
    forall x:real [round RNE x].
876 877
      no_overflow RNE x ->
        x - 0x1p-24 * Abs.abs(x) - 0x1p-150 <= round RNE x <= x + 0x1p-24 * Abs.abs(x) + 0x1p-150
878 879

  lemma round_bound :
880
    forall m:mode, x:real [round m x].
881
      no_overflow m x ->
882
      x - 0x1p-23 * Abs.abs(x) - 0x1p-149 <= round m x <= x + 0x1p-23 * Abs.abs(x) + 0x1p-149
883 884
end

885 886
(** {2 Standard double precision floats (64 bits)} *)

887
module Float64
888
  use int.Int as Int
889
  use real.Real
890

891
  type t = < float 11 53 >
892

893
  constant pow2sb : int = 9007199254740992
894 895 896 897
  constant max_real : real = 0x1.FFFFFFFFFFFFFp1023

  clone export GenericFloat with
    type t = t,
898 899
    constant eb = t'eb,
    constant sb = t'sb,
900
    constant max_real = max_real,
901
    constant pow2sb = pow2sb,
902 903
    function to_real = t'real,
    predicate is_finite = t'isFinite,
904 905 906
    goal eb_gt_1,
    goal sb_gt_1,
    goal max_int,
907 908
    goal pow2sb,
    axiom . (* TODO: "lemma"? "goal"? *)
909

910 911
  lemma round_bound_ne :
    forall x:real [round RNE x].
912 913
      no_overflow RNE x ->
        x - 0x1p-53 * Abs.abs(x) - 0x1p-1075 <= round RNE x <= x + 0x1p-53 * Abs.abs(x) + 0x1p-1075
914 915

  lemma round_bound :
916
    forall m:mode, x:real [round m x].
917
      no_overflow m x ->
918
      x - 0x1p-52 * Abs.abs(x) - 0x1p-1074 <= round m x <= x + 0x1p-52 * Abs.abs(x) + 0x1p-1074
919 920
end

921 922
(** {2 Conversions between float formats} *)

923
module FloatConverter
924

925 926
  use Float64 as Float64
  use Float32 as Float32
927 928 929 930 931 932 933 934 935 936 937

  use export RoundingMode

  function to_float64 mode Float32.t : Float64.t
  function to_float32 mode Float64.t : Float32.t

  lemma round_double_single :
  forall m1 m2:mode, x:real.
    Float64.round m1 (Float32.round m2 x) = Float32.round m2 x

  lemma to_float64_exact:
938 939 940
    forall m:mode, x:Float32.t. Float32.t'isFinite x ->
      Float64.t'isFinite (to_float64 m x)
   /\ Float64.t'real (to_float64 m x) = Float32.t'real x
941 942

  lemma to_float32_conv:
943 944 945 946
    forall m:mode, x:Float64.t. Float64.t'isFinite x ->
      Float32.no_overflow m (Float64.t'real x) ->
        Float32.t'isFinite (to_float32 m x)
     /\ Float32.t'real (to_float32 m x) = Float32.round m (Float64.t'real x)
947 948 949

end

950
module Float32_BV_Converter
951
  use Float32
952 953 954

  clone export Float_BV_Converter with
    type t = t,
955
    predicate is_finite = t'isFinite,
956
    predicate le = (.<=),
957
    function to_real = t'real,
958 959
    function round = round,
    axiom . (* TODO: "lemma"? "goal"? *)
960 961

  axiom of_ubv32_to_real : forall m, x.
962
    t'real (of_ubv32 m x) = round m (FromInt.from_int (BV32.t'int x))
963 964
end

965
module Float64_BV_Converter
966
  use Float64
967 968 969

  clone export Float_BV_Converter with
    type t = t,
970
    predicate is_finite = t'isFinite,
971
    predicate le = (.<=),
972
    function to_real = t'real,
973 974
    function round = round,
    axiom . (* TODO: "lemma"? "goal"? *)
975 976

  axiom of_ubv32_to_real : forall m, x.
977
    t'real (of_ubv32 m x) = FromInt.from_int (BV32.t'int x)
978
end