stdlib/ieee_float.mlw

   1 (** {1 Formalization of Floating-Point Arithmetic}
   2
   3   Full float theory (with infinities and NaN).
   4
   5   A note on intended semantics: we use the same idea as the SMTLIB floating
   6   point theory, that defers any inconsistencies to the "parent" document.
   7   Hence, in doubt, the correct axiomatisation is one that implements
   8   [BTRW14] "An Automatable Formal Semantics for IEEE-754 Floating-Point
   9   Arithmetic", which in turn defers any inconsistencies to IEEE-754.
  10
  11   This theory is split into two parts: the first part talks about IEEE
  12   operations and this is what you should use as a user, the second part is
  13   internal only and is an axiomatisation for the provers that do not
  14   natively support floating point. You should not use any symbols you find
  15   there in your verification conditions as solvers with native floating
  16   point support will leave them uninterpreted.
  17 *)
  18
  19 (** {2 Rounding Modes} *)
  20
  21 module RoundingMode
  22   type mode = RNE | RNA | RTP | RTN | RTZ
  23   (** {h <ul>
  24      <li>RNE : Round Nearest ties to Even
  25      <li>RNA : Round Nearest ties to Away
  26      <li>RTP : Round Towards Positive
  27      <li>RTN : Round Towards Negative
  28      <li>RTZ : Round Towards Zero
  29      </ul>} *)
  30
  31   predicate to_nearest (m:mode) = m = RNE \/ m = RNA
  32 end
  33
  34 module GenericFloat
  35
  36   use int.Int
  37   use bv.Pow2int
  38   use real.Abs as Abs
  39   use real.FromInt as FromInt
  40   use real.Truncate as Truncate
  41   use real.RealInfix
  42   use export RoundingMode
  43
  44   (** {2 Part I - Public Interface}   *)
  45
  46   constant eb : int
  47   (** the number of bits in the exponent.  *)
  48
  49   constant sb : int
  50   (** the number of bits in the significand, including the hidden bit. *)
  51
  52   axiom eb_gt_1: 1 < eb
  53   axiom sb_gt_1: 1 < sb
  54
  55   (** {3 Sorts} *)
  56
  57   type t
  58   (** abstract type to denote floating-point numbers, including the special values
  59       for infinities and NaNs *)
  60
  61   (** {3 Constructors and Constants} *)
  62
  63   val constant zeroF     : t   (** +0.0 *)
  64   (* exp_bias = 2^(eb - 1) - 1 *)
  65   (* max_finite_exp = 2^sb - 2 - exp_bias = exp_bias *)
  66   (* 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   (* max_value = [2^(eb-1); 2^sb - 2] *)
  70
  71   (** {3 Operators} *)
  72
  73   val function add mode t t : t
  74   val function sub mode t t : t
  75   val function mul mode t t : t
  76   val function div mode t t : t
  77     (** The four basic operations, rounded in the given mode *)
  78
  79   val function abs         t   : t   (** Absolute value *)
  80   val function neg         t   : t   (** Opposite *)
  81   val function fma  mode t t t : t   (** Fused multiply-add: x * y + z *)
  82   val function sqrt mode   t   : t   (** Square root *)
  83
  84   let function (.-_) (x:t)   : t = neg x
  85   let function (.+)  (x y:t) : t = add RNE x y
  86   let function (.-)  (x y:t) : t = sub RNE x y
  87   let function (.*)  (x y:t) : t = mul RNE x y
  88   let function (./)  (x y:t) : t = div RNE x y
  89     (** Notations for operations in the default mode RNE *)
  90
  91   val function roundToIntegral mode t : t
  92     (** Rounding to an integer *)
  93
  94   function min       t t : t
  95   function max       t t : t
  96   (** Minimum and Maximum
  97
  98      Note that we have to follow IEEE-754 and SMTLIB here. Two things to
  99      note in particular:
 100
 101      1) min(-0, 0) is either 0 or -0, there is a choice
 102
 103      2) if either argument is NaN then the other argument is returned
 104
 105      Due to the unclear status of min and max in IEEE norm, we
 106      intentionally not provide these function as "val"s to be used in
 107      programs
 108
 109    *)
 110
 111   (** {3 Comparisons} *)
 112
 113   val predicate le t t
 114   val predicate lt t t
 115   let predicate ge (x:t) (y:t) = le y x
 116   let predicate gt (x:t) (y:t) = lt y x
 117   val predicate eq t t
 118   (** equality on floats, different from = since not (eq NaN NaN) *)
 119
 120   let predicate (.<=) (x:t) (y:t) = le x y
 121   let predicate (.<)  (x:t) (y:t) = lt x y
 122   let predicate (.>=) (x:t) (y:t) = ge x y
 123   let predicate (.>)  (x:t) (y:t) = gt x y
 124   let predicate (.=)  (x:t) (y:t) = eq x y
 125   (** Notations *)
 126
 127
 128   (** {3 Classification of numbers} *)
 129
 130   predicate is_normal    t
 131   predicate is_subnormal t
 132   val predicate is_zero      t
 133   val predicate is_infinite  t
 134   val predicate is_nan       t
 135   val predicate is_positive  t
 136   val predicate is_negative  t
 137
 138   (** helper predicate for zeros, normals and subnormals. not defined
 139      so that the axiomatisation below can use it without talking about
 140      subnormals *)
 141   predicate is_finite    t
 142
 143   predicate is_plus_infinity  (x:t) = is_infinite x /\ is_positive x
 144   predicate is_minus_infinity (x:t) = is_infinite x /\ is_negative x
 145   predicate is_plus_zero      (x:t) = is_zero x /\ is_positive x
 146   predicate is_minus_zero     (x:t) = is_zero x /\ is_negative x
 147   predicate is_not_nan        (x:t) = is_finite x \/ is_infinite x
 148
 149   axiom is_not_nan: forall x:t. is_not_nan x <-> not (is_nan x)
 150
 151   axiom is_not_finite: forall x:t.
 152     not (is_finite x) <-> (is_infinite x \/ is_nan x)
 153
 154
 155
 156   (** {3 Conversions from other sorts} *)
 157
 158   (* from bitvec binary interchange                           *)
 159   (*   partly done with from_binary (for literals only)       *)
 160   (* from another fp                 - see FloatConverter     *)
 161   (* from real                                                *)
 162   (*   partly done with (!) (for literals only)               *)
 163   (* from unsigned integer bitvector - see Float_BV_Converter *)
 164   (* from signed integer bitvector                            *)
 165
 166   (** {3 Conversions to other sorts} *)
 167
 168   (* to unsigned integer bitvector   - see Float_BV_Converter *)
 169   (* to signed integer bitvector                              *)
 170   (* to real *)
 171   function to_real   t    : real
 172
 173   (** {2 Part II - Private Axiomatisation} *)
 174
 175   (** {3 Constructors and Constants} *)
 176
 177   axiom zeroF_is_positive : is_positive zeroF
 178   axiom zeroF_is_zero     : is_zero zeroF
 179
 180   axiom zero_to_real : forall x [is_zero x].
 181     is_zero x <-> is_finite x /\ to_real x = 0.0
 182
 183   (** {3 Conversions from other sorts} *)
 184
 185   (* with mathematical int *)
 186   (* note that these conversions do not feature in SMTLIB *)
 187
 188   (* intended semantics for of_int are the same as (_ to_fp eb sb) with a   *)
 189   (* suitably sized bitvector, large enough to hold x                       *)
 190   (* note values >= than the below should result in infinities              *)
 191   (*    float32 : 0x1ffffff * 2^103                                         *)
 192   (*    float64 : 0x3fffffffffffff * 2^970                                  *)
 193   (* also note that this function never yields a subnormal, or a NaN, or -0 *)
 194   function of_int (m:mode) (x:int) : t
 195
 196   (** {3 Conversions to other sorts} *)
 197
 198   (* Intended semantics for to_int are the same as (_ fp.to_sbv) with a    *)
 199   (* suitably sized bitvector. Safe minimum sizes are given below:         *)
 200   (*    float32 : 129                                                      *)
 201   (*    float64 : 1025                                                     *)
 202   (* In particular this function should be uninterpreted for infinities    *)
 203   (* and NaN. Take care that no conclusion can be made on the result based *)
 204   (* on the size of the bitvector chosen in those cases, i.e. this should  *)
 205   (* not hold:                                                             *)
 206   (*    to_int +INF < 2 ** 2048     // nope                                *)
 207   (*    to_int +INF > 0             // nope                                *)
 208   function to_int (m:mode) (x:t) : int
 209
 210   axiom zero_of_int : forall m. zeroF = of_int m 0
 211
 212   (** {3 Arithmetic} *)
 213
 214   (* The intended meaning for round is the rounding for floating point as  *)
 215   (* described on p17 of IEEE-754. For results where the corresponding     *)
 216   (* floating point result would be infinity or NaN this function should   *)
 217   (* be uninterpreted.                                                     *)
 218   (*                                                                       *)
 219   (* Note that this means round (+INF) > 0 is not true.                    *)
 220   (* Note also this means round (2*INF) > round (INF) is not true either.  *)
 221   function round mode real : real
 222
 223   constant max_real : real (* defined when cloning *)
 224   constant max_int  : int  (* defined when cloning *)
 225
 226   constant emax : int = pow2 (eb - 1)
 227
 228   axiom max_int_spec : max_int = pow2 emax - pow2 (emax - sb)
 229   axiom max_real_int: max_real = FromInt.from_int max_int
 230
 231   predicate in_range (x:real) = -. max_real <=. x <=. max_real
 232
 233   predicate in_int_range (i:int) = - max_int <= i <= max_int
 234
 235   axiom is_finite: forall x:t. is_finite x -> in_range (to_real x)
 236
 237   (* used as a condition to propagate is_finite *)
 238   predicate no_overflow (m:mode) (x:real) = in_range (round m x)
 239
 240   (* Axioms on round *)
 241
 242   axiom Bounded_real_no_overflow [@W:non_conservative_extension:N] :
 243     forall m:mode, x:real. in_range x -> no_overflow m x
 244
 245   axiom Round_monotonic :
 246     forall m:mode, x y:real. x <=. y -> round m x <=. round m y
 247
 248   axiom Round_idempotent :
 249     forall m1 m2:mode, x:real. round m1 (round m2 x) = round m2 x
 250
 251   axiom Round_to_real :
 252     forall m:mode, x:t. is_finite x -> round m (to_real x) = to_real x
 253
 254   (** rounding up and down *)
 255   axiom Round_down_le:
 256     forall x:real. round RTN x <=. x
 257   axiom Round_up_ge:
 258     forall x:real. round RTP x >=. x
 259   axiom Round_down_neg:
 260     forall x:real. round RTN (-.x) = -. round RTP x
 261   axiom Round_up_neg:
 262     forall x:real. round RTP (-.x) = -. round RTN x
 263
 264   (* The biggest representable integer whose predecessor (i.e. -1) is  representable *)
 265   constant pow2sb : int (* defined when cloning *)
 266   axiom pow2sb: pow2sb = pow2 sb
 267
 268   (* range in which every integer is representable *)
 269   predicate in_safe_int_range (i: int) = - pow2sb <= i <= pow2sb
 270
 271   (* round and integers *)
 272
 273   axiom Exact_rounding_for_integers:
 274     forall m:mode, i:int.
 275       in_safe_int_range i ->
 276         round m (FromInt.from_int i) = FromInt.from_int i
 277
 278   (** {3 Comparisons} *)
 279
 280   (** Comparison predicates *)
 281
 282   predicate same_sign (x y : t) =
 283     (is_positive x /\ is_positive y) \/ (is_negative x /\ is_negative y)
 284   predicate diff_sign (x y : t) =
 285     (is_positive x /\ is_negative y) \/ (is_negative x /\ is_positive y)
 286
 287   axiom feq_eq: forall x y.
 288     is_finite x -> is_finite y -> not (is_zero x) -> x .= y -> x = y
 289
 290   axiom eq_feq: forall x y.
 291     is_finite x -> is_finite y -> x = y -> x .= y
 292
 293   axiom eq_refl: forall x. is_finite x -> x .= x
 294
 295   axiom eq_sym :
 296     forall x y. x .= y -> y .= x
 297
 298   axiom eq_trans :
 299     forall x y z. x .= y -> y .= z -> x .= z
 300
 301   axiom eq_zero: zeroF .= (.- zeroF)
 302
 303   axiom eq_to_real_finite: forall x y.
 304     is_finite x /\ is_finite y -> (x .= y <-> to_real x = to_real y)
 305
 306   axiom eq_special: forall x y. x .= y ->
 307        (is_not_nan x  /\ is_not_nan y
 308     /\ ((is_finite x /\ is_finite y)
 309         \/ (is_infinite x /\ is_infinite y /\ same_sign x y)))
 310
 311   axiom lt_finite: forall x y [lt x y].
 312     is_finite x /\ is_finite y -> (lt x y <-> to_real x <. to_real y)
 313
 314   axiom le_finite: forall x y [le x y].
 315     is_finite x /\ is_finite y -> (le x y <-> to_real x <=. to_real y)
 316
 317   lemma le_lt_trans:
 318     forall x y z:t. x .<= y /\ y .< z -> x .< z
 319
 320   lemma lt_le_trans:
 321     forall x y z:t. x .< y /\ y .<= z -> x .< z
 322
 323   lemma le_ge_asym:
 324     forall x y:t. x .<= y /\ x .>= y -> x .= y
 325
 326   lemma not_lt_ge: forall x y:t.
 327     not (x .< y) /\ is_not_nan x /\ is_not_nan y -> x .>= y
 328
 329   lemma not_gt_le: forall x y:t.
 330     not (x .> y) /\ is_not_nan x /\ is_not_nan y -> x .<= y
 331
 332   axiom le_special: forall x y [le x y]. le x y ->
 333       ((is_finite x /\ is_finite y)
 334    \/ ((is_minus_infinity x /\ is_not_nan y)
 335    \/  (is_not_nan x /\ is_plus_infinity y)))
 336
 337   axiom lt_special: forall x y [lt x y]. lt x y ->
 338       ((is_finite x /\ is_finite y)
 339    \/ ((is_minus_infinity x /\ is_not_nan y /\ not (is_minus_infinity y))
 340    \/  (is_not_nan x /\ not (is_plus_infinity x) /\ is_plus_infinity y)))
 341
 342   axiom lt_lt_finite: forall x y z. lt x y -> lt y z -> is_finite y
 343
 344   (*  lemmas on sign *)
 345   axiom positive_to_real: forall x[is_positive x|to_real x >=. 0.0].
 346     is_finite x -> is_positive x -> to_real x >=. 0.0
 347   axiom to_real_positive: forall x[is_positive x].
 348     is_finite x -> to_real x >. 0.0 -> is_positive x
 349
 350   axiom negative_to_real: forall x [is_negative x|to_real x <=. 0.0].
 351     is_finite x -> is_negative x -> to_real x <=. 0.0
 352   axiom to_real_negative: forall x [is_negative x].
 353     is_finite x -> to_real x <. 0.0 -> is_negative x
 354
 355   axiom negative_xor_positive: forall x.
 356     not (is_positive x /\ is_negative x)
 357   axiom negative_or_positive: forall x.
 358     is_not_nan x -> is_positive x \/ is_negative x
 359
 360   lemma diff_sign_trans:
 361     forall x y z:t. (diff_sign x y /\ diff_sign y z) -> same_sign x z
 362
 363   lemma diff_sign_product:
 364     forall x y:t.
 365       (is_finite x /\ is_finite y /\ to_real x *. to_real y <. 0.0) ->
 366         diff_sign x y
 367
 368   lemma same_sign_product:
 369     forall x y:t.
 370       (is_finite x /\ is_finite y /\ same_sign x y) ->
 371         to_real x *. to_real y >=. 0.0
 372
 373   predicate product_sign (z x y: t) =
 374     (same_sign x y -> is_positive z) /\ (diff_sign x y -> is_negative z)
 375
 376   (** {3 Overflow} *)
 377
 378   (* This predicate is used to tell what is the result of a rounding
 379      in case of overflow in the axioms specifying add/sub/mul and fma
 380      *)
 381   predicate overflow_value (m:mode) (x:t) =
 382     match m with
 383     | RTN -> if is_positive x then is_finite x /\ to_real x = max_real
 384                               else is_infinite x
 385     | RTP -> if is_positive x then is_infinite x
 386                               else is_finite x /\ to_real x = -. max_real
 387     | RTZ -> if is_positive x then is_finite x /\ to_real x = max_real
 388                               else is_finite x /\ to_real x = -. max_real
 389     | (RNA | RNE) -> is_infinite x
 390     end
 391
 392   (* This predicate is used to tell what is the sign of zero in the
 393      axioms specifying add and sub *)
 394   predicate sign_zero_result (m:mode) (x:t) =
 395     is_zero x ->
 396     match m with
 397     | RTN -> is_negative x
 398     | _   -> is_positive x
 399     end
 400
 401   (** {3 binary operations} *)
 402
 403   axiom add_finite: forall m:mode, x y:t [add m x y].
 404     is_finite x -> is_finite y -> no_overflow m (to_real x +. to_real y) ->
 405     is_finite (add m x y) /\
 406     to_real (add m x y) = round m (to_real x +. to_real y)
 407
 408   lemma add_finite_rev: forall m:mode, x y:t [add m x y].
 409     is_finite (add m x y) ->
 410     is_finite x /\ is_finite y
 411
 412   lemma add_finite_rev_n: forall m:mode, x y:t [add m x y].
 413     to_nearest m ->
 414     is_finite (add m x y) ->
 415     no_overflow m (to_real x +. to_real y) /\
 416     to_real (add m x y) = round m (to_real x +. to_real y)
 417
 418   axiom sub_finite: forall m:mode, x y:t [sub m x y].
 419     is_finite x -> is_finite y -> no_overflow m (to_real x -. to_real y) ->
 420     is_finite (sub m x y) /\
 421     to_real (sub m x y) = round m (to_real x -. to_real y)
 422
 423   lemma sub_finite_rev: forall m:mode, x y:t [sub m x y].
 424     is_finite (sub m x y) ->
 425     is_finite x /\ is_finite y
 426
 427   lemma sub_finite_rev_n: forall m:mode, x y:t [sub m x y].
 428     to_nearest m ->
 429     is_finite (sub m x y) ->
 430     no_overflow m (to_real x -. to_real y) /\
 431     to_real (sub m x y) = round m (to_real x -. to_real y)
 432
 433   axiom mul_finite: forall m:mode, x y:t [mul m x y].
 434     is_finite x -> is_finite y -> no_overflow m (to_real x *. to_real y) ->
 435     is_finite (mul m x y) /\
 436     to_real (mul m x y) = round m (to_real x *. to_real y)
 437
 438   lemma mul_finite_rev: forall m:mode, x y:t [mul m x y].
 439     is_finite (mul m x y) ->
 440     is_finite x /\ is_finite y
 441
 442   lemma mul_finite_rev_n: forall m:mode, x y:t [mul m x y].
 443     to_nearest m ->
 444     is_finite (mul m x y) ->
 445     no_overflow m (to_real x *. to_real y) /\
 446     to_real (mul m x y) = round m (to_real x *. to_real y)
 447
 448   axiom div_finite: forall m:mode, x y:t [div m x y].
 449     is_finite x -> is_finite y ->
 450     not is_zero y -> no_overflow m (to_real x /. to_real y) ->
 451     is_finite (div m x y) /\
 452     to_real (div m x y) = round m (to_real x /. to_real y)
 453
 454   lemma div_finite_rev: forall m:mode, x y:t [div m x y].
 455     is_finite (div m x y) ->
 456     (is_finite x /\ is_finite y /\ not is_zero y) \/
 457     (is_finite x /\ is_infinite y /\ to_real (div m x y) = 0.)
 458
 459   lemma div_finite_rev_n: forall m:mode, x y:t [div m x y].
 460     to_nearest m ->
 461     is_finite (div m x y) -> is_finite y ->
 462     no_overflow m (to_real x /. to_real y) /\
 463     to_real (div m x y) = round m (to_real x /. to_real y)
 464
 465   axiom neg_finite: forall x:t [neg x].
 466     is_finite x ->
 467     is_finite (neg x) /\
 468     to_real (neg x) = -. to_real x
 469
 470   lemma neg_finite_rev: forall x:t [neg x].
 471     is_finite (neg x) ->
 472     is_finite x /\
 473     to_real (neg x) = -. to_real x
 474
 475   axiom abs_finite: forall x:t [abs x].
 476     is_finite x ->
 477     is_finite (abs x) /\
 478     to_real (abs x) = Abs.abs (to_real x) /\
 479     is_positive (abs x)
 480
 481   lemma abs_finite_rev: forall x:t [abs x].
 482     is_finite (abs x) ->
 483     is_finite x /\
 484     to_real (abs x) = Abs.abs (to_real x)
 485
 486   axiom abs_universal : forall x:t [abs x]. not (is_negative (abs x))
 487
 488   axiom fma_finite: forall m:mode, x y z:t [fma m x y z].
 489     is_finite x -> is_finite y -> is_finite z ->
 490     no_overflow m (to_real x *. to_real y +. to_real z) ->
 491     is_finite (fma m x y z) /\
 492     to_real (fma m x y z) = round m (to_real x *. to_real y +. to_real z)
 493
 494   lemma fma_finite_rev: forall m:mode, x y z:t [fma m x y z].
 495     is_finite (fma m x y z) ->
 496     is_finite x /\ is_finite y /\ is_finite z
 497
 498   lemma fma_finite_rev_n: forall m:mode, x y z:t [fma m x y z].
 499     to_nearest m ->
 500     is_finite (fma m x y z) ->
 501     no_overflow m (to_real x *. to_real y +. to_real z) /\
 502     to_real (fma m x y z) = round m (to_real x *. to_real y +. to_real z)
 503
 504   use real.Square as S
 505
 506   axiom sqrt_finite: forall m:mode, x:t [sqrt m x].
 507     is_finite x -> to_real x >=. 0. ->
 508     is_finite (sqrt m x) /\
 509     to_real (sqrt m x) = round m (S.sqrt (to_real x))
 510
 511   lemma sqrt_finite_rev: forall m:mode, x:t [sqrt m x].
 512     is_finite (sqrt m x) ->
 513     is_finite x /\ to_real x >=. 0. /\
 514     to_real (sqrt m x) = round m (S.sqrt (to_real x))
 515
 516   predicate same_sign_real (x:t) (r:real) =
 517     (is_positive x /\ r >. 0.0) \/ (is_negative x /\ r <. 0.0)
 518
 519   axiom add_special: forall m:mode, x y:t [add m x y].
 520     let r = add m x y in
 521     (is_nan x \/ is_nan y -> is_nan r)
 522     /\
 523     (is_finite x /\ is_infinite y -> is_infinite r /\ same_sign r y)
 524     /\
 525     (is_infinite x /\ is_finite y -> is_infinite r /\ same_sign r x)
 526     /\
 527     (is_infinite x /\ is_infinite y /\ same_sign x y
 528       -> is_infinite r /\ same_sign r x)
 529     /\
 530     (is_infinite x /\ is_infinite y /\ diff_sign x y -> is_nan r)
 531     /\
 532     (is_finite x /\ is_finite y /\ not no_overflow m (to_real x +. to_real y)
 533       -> same_sign_real r (to_real x +. to_real y) /\ overflow_value m r)
 534     /\
 535     (is_finite x /\ is_finite y
 536       -> if same_sign x y then same_sign r x else sign_zero_result m r)
 537
 538   axiom sub_special: forall m:mode, x y:t [sub m x y].
 539     let r = sub m x y in
 540     (is_nan x \/ is_nan y -> is_nan r)
 541     /\
 542     (is_finite x /\ is_infinite y -> is_infinite r /\ diff_sign r y)
 543     /\
 544     (is_infinite x /\ is_finite y -> is_infinite r /\ same_sign r x)
 545     /\
 546     (is_infinite x /\ is_infinite y /\ same_sign x y -> is_nan r)
 547     /\
 548     (is_infinite x /\ is_infinite y /\ diff_sign x y
 549       -> is_infinite r /\ same_sign r x)
 550     /\
 551     (is_finite x /\ is_finite y /\ not no_overflow m (to_real x -. to_real y)
 552       -> same_sign_real r (to_real x -. to_real y) /\ overflow_value m r)
 553     /\
 554     (is_finite x /\ is_finite y
 555       -> if diff_sign x y then same_sign r x else sign_zero_result m r)
 556
 557   axiom mul_special: forall m:mode, x y:t [mul m x y].
 558     let r = mul m x y in
 559        (is_nan x \/ is_nan y -> is_nan r)
 560     /\ (is_zero x /\ is_infinite y -> is_nan r)
 561     /\ (is_finite x /\ is_infinite y /\ not (is_zero x)
 562          -> is_infinite r)
 563     /\ (is_infinite x /\ is_zero y -> is_nan r)
 564     /\ (is_infinite x /\ is_finite y /\ not (is_zero  y)
 565          -> is_infinite r)
 566     /\ (is_infinite x /\ is_infinite y -> is_infinite r)
 567     /\ (is_finite x /\ is_finite y /\ not no_overflow m (to_real x *. to_real y)
 568          -> overflow_value m r)
 569     /\ (not is_nan r -> product_sign r x y)
 570
 571   axiom div_special: forall m:mode, x y:t [div m x y].
 572     let r = div m x y in
 573        (is_nan x \/ is_nan y -> is_nan r)
 574     /\ (is_finite x /\ is_infinite y -> is_zero r)
 575     /\ (is_infinite x /\ is_finite y -> is_infinite r)
 576     /\ (is_infinite x /\ is_infinite y -> is_nan r)
 577     /\ (is_finite x /\ is_finite y /\ not (is_zero y) /\
 578         not no_overflow m (to_real x /. to_real y)
 579          -> overflow_value m r)
 580     /\ (is_finite x /\ is_zero y /\ not (is_zero x)
 581          -> is_infinite r)
 582     /\ (is_zero x /\ is_zero y -> is_nan r)
 583     /\ (not is_nan r -> product_sign r x y)
 584
 585   axiom neg_special: forall x:t [neg x].
 586        (is_nan x -> is_nan (neg x))
 587     /\ (is_infinite x -> is_infinite (neg x))
 588     /\ (not is_nan x -> diff_sign x (neg x))
 589
 590   axiom abs_special: forall x:t [abs x].
 591        (is_nan x -> is_nan (abs x))
 592     /\ (is_infinite x -> is_infinite (abs x))
 593     /\ (not is_nan x -> is_positive (abs x))
 594
 595   axiom fma_special: forall m:mode, x y z:t [fma m x y z].
 596     let r = fma m x y z in
 597        (is_nan x \/ is_nan y \/ is_nan z -> is_nan r)
 598     /\ (is_zero x /\ is_infinite y -> is_nan r)
 599     /\ (is_infinite x /\ is_zero y -> is_nan r)
 600     /\ (is_finite x /\ not (is_zero x) /\ is_infinite y /\ is_finite z
 601         -> is_infinite r /\ product_sign r x y)
 602     /\ (is_finite x /\ not (is_zero x) /\ is_infinite y /\ is_infinite z
 603         -> (if product_sign z x y then is_infinite r /\ same_sign r z
 604             else is_nan r))
 605     /\ (is_infinite x /\ is_finite y /\ not (is_zero y) /\ is_finite z
 606         -> is_infinite r /\ product_sign r x y)
 607     /\ (is_infinite x /\ is_finite y /\ not (is_zero y) /\ is_infinite z
 608         -> (if product_sign z x y then is_infinite r /\ same_sign r z
 609             else is_nan r))
 610     /\ (is_infinite x /\ is_infinite y /\ is_finite z
 611         -> is_infinite r /\ product_sign r x y)
 612     /\ (is_finite x /\ is_finite y /\ is_infinite z
 613         -> is_infinite r /\ same_sign r z)
 614     /\ (is_infinite x /\ is_infinite y /\ is_infinite z
 615         -> (if product_sign z x y then is_infinite r /\ same_sign r z
 616             else is_nan r))
 617     /\ (is_finite x /\ is_finite y /\ is_finite z /\
 618         not no_overflow m (to_real x *. to_real y +. to_real z)
 619         -> same_sign_real r (to_real x *. to_real y +. to_real z)
 620         /\ overflow_value m r)
 621     /\ (is_finite x /\ is_finite y /\ is_finite z
 622         -> if product_sign z x y then same_sign r z
 623            else (to_real x *. to_real y +. to_real z = 0.0 ->
 624                  if m = RTN then is_negative r else is_positive r))
 625
 626   axiom sqrt_special: forall m:mode, x:t [sqrt m x].
 627     let r = sqrt m x in
 628        (is_nan x -> is_nan r)
 629     /\ (is_plus_infinity x -> is_plus_infinity r)
 630     /\ (is_minus_infinity x -> is_nan r)
 631     /\ (is_finite x /\ to_real x <. 0.0 -> is_nan r)
 632     /\ (is_zero x -> same_sign r x)
 633     /\ (is_finite x /\ to_real x >. 0.0 -> is_positive r)
 634
 635   (* exact arithmetic with integers *)
 636
 637   axiom of_int_add_exact: forall m n, i j.
 638     in_safe_int_range i -> in_safe_int_range j ->
 639     in_safe_int_range (i + j) -> eq (of_int m (i + j)) (add n (of_int m i) (of_int m j))
 640
 641   axiom of_int_sub_exact: forall m n, i j.
 642     in_safe_int_range i -> in_safe_int_range j ->
 643     in_safe_int_range (i - j) -> eq (of_int m (i - j)) (sub n (of_int m i) (of_int m j))
 644
 645   axiom of_int_mul_exact: forall m n, i j.
 646     in_safe_int_range i -> in_safe_int_range j ->
 647     in_safe_int_range (i * j) -> eq (of_int m (i * j)) (mul n (of_int m i) (of_int m j))
 648
 649   (* min and max *)
 650
 651   lemma Min_r : forall x y:t. y .<= x -> (min x y) .= y
 652   lemma Min_l : forall x y:t. x .<= y -> (min x y) .= x
 653   lemma Max_r : forall x y:t. y .<= x -> (max x y) .= x
 654   lemma Max_l : forall x y:t. x .<= y -> (max x y) .= y
 655
 656   (* _____________ *)
 657
 658   use real.Truncate as Truncate
 659
 660   (* This predicate specify if a float is finite and is an integer *)
 661   predicate is_int (x:t)
 662
 663   (** characterizing integers *)
 664
 665   (* by construction *)
 666   axiom zeroF_is_int: is_int zeroF
 667
 668   axiom of_int_is_int: forall m, x.
 669     in_int_range x -> is_int (of_int m x)
 670
 671   axiom big_float_is_int: forall m i.
 672     is_finite i ->
 673     i .<= neg (of_int m pow2sb) \/ (of_int m pow2sb) .<= i ->
 674       is_int i
 675
 676   axiom roundToIntegral_is_int: forall m:mode, x:t. is_finite x ->
 677     is_int (roundToIntegral m x)
 678
 679   (* by propagation *)
 680   axiom eq_is_int: forall x y. eq x y -> is_int x -> is_int y
 681
 682   axiom add_int: forall x y m. is_int x -> is_int y ->
 683     is_finite (add m x y) -> is_int (add m x y)
 684
 685   axiom sub_int: forall x y m. is_int x -> is_int y ->
 686     is_finite (sub m x y) -> is_int (sub m x y)
 687
 688   axiom mul_int: forall x y m. is_int x -> is_int y ->
 689     is_finite (mul m x y) -> is_int (mul m x y)
 690
 691   axiom fma_int: forall x y z m. is_int x -> is_int y -> is_int z ->
 692     is_finite (fma m x y z) -> is_int (fma m x y z)
 693
 694   axiom neg_int: forall x. is_int x -> is_int (neg x)
 695
 696   axiom abs_int: forall x. is_int x -> is_int (abs x)
 697
 698   (** basic properties of float integers *)
 699
 700   axiom is_int_of_int: forall x m m'.
 701     is_int x -> eq x (of_int m' (to_int m x))
 702
 703   axiom is_int_to_int: forall m x.
 704     is_int x -> in_int_range (to_int m x)
 705
 706   axiom is_int_is_finite: forall x.
 707     is_int x -> is_finite x
 708
 709   axiom int_to_real: forall m x.
 710     is_int x -> to_real x = FromInt.from_int (to_int m x)
 711
 712 (*  axiom int_mode: forall m1 m2 x.
 713     is_int x -> to_int m1 x = to_int m2 x  etc ...*)
 714
 715   (** rounding ints *)
 716
 717   axiom truncate_int: forall m:mode, i:t. is_int i ->
 718     roundToIntegral m i .= i
 719
 720   (** truncate *)
 721
 722   axiom truncate_neg: forall x:t.
 723     is_finite x -> is_negative x -> roundToIntegral RTZ x = roundToIntegral RTP x
 724
 725   axiom truncate_pos: forall x:t.
 726     is_finite x -> is_positive x -> roundToIntegral RTZ x = roundToIntegral RTN x
 727
 728   (** ceil *)
 729
 730   axiom ceil_le: forall x:t. is_finite x -> x .<= (roundToIntegral RTP x)
 731
 732   axiom ceil_lest: forall x y:t. x .<= y /\ is_int y -> (roundToIntegral RTP x) .<= y
 733
 734   axiom ceil_to_real: forall x:t.
 735     is_finite x ->
 736       to_real (roundToIntegral RTP x) = FromInt.from_int (Truncate.ceil (to_real x))
 737
 738   axiom ceil_to_int: forall m:mode, x:t.
 739     is_finite x ->
 740       to_int m (roundToIntegral RTP x) = Truncate.ceil (to_real x)
 741
 742   (** floor *)
 743
 744   axiom floor_le: forall x:t. is_finite x -> (roundToIntegral RTN x) .<= x
 745
 746   axiom floor_lest: forall x y:t. y .<= x /\ is_int y -> y .<= (roundToIntegral RTN x)
 747
 748   axiom floor_to_real: forall x:t.
 749     is_finite x ->
 750       to_real (roundToIntegral RTN x) = FromInt.from_int (Truncate.floor (to_real x))
 751
 752   axiom floor_to_int: forall m:mode, x:t.
 753     is_finite x ->
 754       to_int m (roundToIntegral RTN x) = Truncate.floor (to_real x)
 755
 756   (* Rna *)
 757
 758   axiom RNA_down:
 759     forall x:t. (x .- (roundToIntegral RTN x)) .< ((roundToIntegral RTP x) .- x) ->
 760       roundToIntegral RNA x = roundToIntegral RTN x
 761
 762   axiom RNA_up:
 763     forall x:t. (x .- (roundToIntegral RTN x)) .> ((roundToIntegral RTP x) .- x) ->
 764       roundToIntegral RNA x = roundToIntegral RTP x
 765
 766   axiom RNA_down_tie:
 767     forall x:t. (x .- (roundToIntegral RTN x)) .= ((roundToIntegral RTP x) .- x) ->
 768       is_negative x -> roundToIntegral RNA x = roundToIntegral RTN x
 769
 770   axiom RNA_up_tie:
 771     forall x:t. ((roundToIntegral RTP x) .- x) .= (x .- (roundToIntegral RTN x)) ->
 772     is_positive x -> roundToIntegral RNA x = roundToIntegral RTP x
 773
 774   (* to_int *)
 775   axiom to_int_roundToIntegral: forall m:mode, x:t.
 776     to_int m x = to_int m (roundToIntegral m x)
 777
 778   axiom to_int_monotonic: forall m:mode, x y:t.
 779     is_finite x -> is_finite y -> le x y -> to_int m x <= to_int m y
 780
 781   axiom to_int_of_int: forall m:mode, i:int.
 782     in_safe_int_range i ->
 783       to_int m (of_int m i) = i
 784
 785   axiom eq_to_int: forall m, x y. is_finite x -> x .= y ->
 786     to_int m x = to_int m y
 787
 788   axiom neg_to_int: forall m x.
 789     is_int x -> to_int m (neg x) = - (to_int m x)
 790
 791   axiom roundToIntegral_is_finite  : forall m:mode, x:t. is_finite x ->
 792     is_finite (roundToIntegral m x)
 793 end
 794
 795 (** {2 Conversions to/from bitvectors} *)
 796
 797 module Float_BV_Converter
 798   use bv.BV8 as BV8
 799   use bv.BV16 as BV16
 800   use bv.BV32 as BV32
 801   use bv.BV64 as BV64
 802   use RoundingMode
 803
 804   type t (* the underlying float type, to be cloned *)
 805   predicate is_finite t
 806   predicate le t t
 807   function to_real         t : real
 808   function round   mode real : real
 809
 810   (* convert from signed bitvector *)
 811   val function of_sbv8  mode BV8.t  : t
 812   val function of_sbv16 mode BV16.t : t
 813   val function of_sbv32 mode BV32.t : t
 814   val function of_sbv64 mode BV64.t : t
 815
 816   (* convert to signed bitvector *)
 817   val function to_sbv8  mode BV8.t  : t
 818   val function to_sbv16 mode BV16.t : t
 819   val function to_sbv32 mode BV32.t : t
 820   val function to_sbv64 mode BV64.t : t
 821
 822   (* convert from unsigned bitvector *)
 823   val function of_ubv8  mode BV8.t  : t
 824   val function of_ubv16 mode BV16.t : t
 825   val function of_ubv32 mode BV32.t : t
 826   val function of_ubv64 mode BV64.t : t
 827
 828   (* convert to unsigned bitvector *)
 829   val function to_ubv8  mode t : BV8.t
 830   val function to_ubv16 mode t : BV16.t
 831   val function to_ubv32 mode t : BV32.t
 832   val function to_ubv64 mode t : BV64.t
 833
 834   use real.RealInfix
 835   use real.FromInt as FromInt
 836
 837   (** of unsigned bv axioms  *)
 838   (* only true for big enough floats...  *)
 839
 840   axiom of_ubv8_is_finite : forall m, x. is_finite (of_ubv8  m x)
 841   axiom of_ubv16_is_finite: forall m, x. is_finite (of_ubv16 m x)
 842   axiom of_ubv32_is_finite: forall m, x. is_finite (of_ubv32 m x)
 843   axiom of_ubv64_is_finite: forall m, x. is_finite (of_ubv64 m x)
 844
 845   axiom of_ubv8_monotonic :
 846     forall m, x y. BV8.ule  x y -> le (of_ubv8  m x) (of_ubv8  m y)
 847   axiom of_ubv16_monotonic:
 848     forall m, x y. BV16.ule x y -> le (of_ubv16 m x) (of_ubv16 m y)
 849   axiom of_ubv32_monotonic:
 850     forall m, x y. BV32.ule x y -> le (of_ubv32 m x) (of_ubv32 m y)
 851   axiom of_ubv64_monotonic:
 852     forall m, x y. BV64.ule x y -> le (of_ubv64 m x) (of_ubv64 m y)
 853
 854   axiom of_ubv8_to_real : forall m, x.
 855     to_real (of_ubv8  m x) = FromInt.from_int (BV8.t'int x)
 856   axiom of_ubv16_to_real: forall m, x.
 857     to_real (of_ubv16 m x) = FromInt.from_int (BV16.t'int x)
 858   (* of_ubv32_to_real is defined at cloning *)
 859   axiom of_ubv64_to_real: forall m, x.
 860     to_real (of_ubv64 m x) = round m (FromInt.from_int (BV64.t'int x))
 861 end
 862
 863 (** {2 Standard simple precision floats (32 bits)} *)
 864
 865 module Float32
 866   use int.Int as Int
 867   use real.Real
 868
 869   type t = < float 8 24 >
 870
 871   constant pow2sb : int = 16777216
 872   constant max_int : int = 0xFFFF_FF00_0000_0000_0000_0000_0000_0000
 873   constant max_real : real = 0x1.FFFFFEp127
 874
 875   clone export GenericFloat with
 876     type t = t,
 877     constant eb = t'eb,
 878     constant sb = t'sb,
 879     constant pow2sb = pow2sb,
 880     constant max_int = max_int,
 881     constant max_real = max_real,
 882     function to_real = t'real,
 883     predicate is_finite = t'isFinite,
 884     (* the following are lemmas and not goals, because we want to
 885        prove them in the realizations. See also
 886        [https://gitlab.inria.fr/why3/why3/-/issues/664] *)
 887     lemma eb_gt_1,
 888     lemma sb_gt_1,
 889     lemma max_int_spec,
 890     lemma max_real_int,
 891     lemma pow2sb,
 892     axiom . (* TODO: "lemma"? "goal"? *)
 893
 894   lemma round_bound_ne :
 895     forall x:real [round RNE x].
 896       no_overflow RNE x ->
 897         x - 0x1p-24 * Abs.abs(x) - 0x1p-150 <= round RNE x <= x + 0x1p-24 * Abs.abs(x) + 0x1p-150
 898
 899   lemma round_bound :
 900     forall m:mode, x:real [round m x].
 901       no_overflow m x ->
 902       x - 0x1p-23 * Abs.abs(x) - 0x1p-149 <= round m x <= x + 0x1p-23 * Abs.abs(x) + 0x1p-149
 903 end
 904
 905 (** {2 Standard double precision floats (64 bits)} *)
 906
 907 module Float64
 908   use int.Int as Int
 909   use real.Real
 910
 911   type t = < float 11 53 >
 912
 913   constant pow2sb : int = 9007199254740992  (* 242 *)
 914   constant max_int : int = 0xFFFF_FFFF_FFFF_F800_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000_0000
 915   constant max_real : real = 0x1.FFFFFFFFFFFFFp1023
 916
 917   clone export GenericFloat with
 918     type t = t,
 919     constant eb = t'eb,
 920     constant sb = t'sb,
 921     constant pow2sb = pow2sb,
 922     constant max_int = max_int,
 923     constant max_real = max_real,
 924     function to_real = t'real,
 925     predicate is_finite = t'isFinite,
 926     (* the following are lemmas and not goals, because we want to
 927        prove them in the realizations. See also
 928        [https://gitlab.inria.fr/why3/why3/-/issues/664] *)
 929     lemma eb_gt_1,
 930     lemma sb_gt_1,
 931     lemma max_int_spec,
 932     lemma max_real_int,
 933     lemma pow2sb,
 934     axiom . (* TODO: "lemma"? "goal"? *)
 935
 936   lemma round_bound_ne :
 937     forall x:real [round RNE x].
 938       no_overflow RNE x ->
 939         x - 0x1p-53 * Abs.abs(x) - 0x1p-1075 <= round RNE x <= x + 0x1p-53 * Abs.abs(x) + 0x1p-1075
 940
 941   lemma round_bound :
 942     forall m:mode, x:real [round m x].
 943       no_overflow m x ->
 944       x - 0x1p-52 * Abs.abs(x) - 0x1p-1074 <= round m x <= x + 0x1p-52 * Abs.abs(x) + 0x1p-1074
 945 end
 946
 947 (** {2 Conversions between float formats} *)
 948
 949 module FloatConverter
 950
 951   use Float64 as Float64
 952   use Float32 as Float32
 953
 954   use export RoundingMode
 955
 956   function to_float64 mode Float32.t : Float64.t
 957   function to_float32 mode Float64.t : Float32.t
 958
 959   lemma round_double_single :
 960   forall m1 m2:mode, x:real.
 961     Float64.round m1 (Float32.round m2 x) = Float32.round m2 x
 962
 963   lemma to_float64_exact:
 964     forall m:mode, x:Float32.t. Float32.t'isFinite x ->
 965       Float64.t'isFinite (to_float64 m x)
 966    /\ Float64.t'real (to_float64 m x) = Float32.t'real x
 967
 968   lemma to_float32_conv:
 969     forall m:mode, x:Float64.t. Float64.t'isFinite x ->
 970       Float32.no_overflow m (Float64.t'real x) ->
 971         Float32.t'isFinite (to_float32 m x)
 972      /\ Float32.t'real (to_float32 m x) = Float32.round m (Float64.t'real x)
 973
 974 end
 975
 976 module Float32_BV_Converter
 977   use Float32
 978
 979   clone export Float_BV_Converter with
 980     type t = t,
 981     predicate is_finite = t'isFinite,
 982     predicate le = (.<=),
 983     function to_real = t'real,
 984     function round = round,
 985     axiom . (* TODO: "lemma"? "goal"? *)
 986
 987   axiom of_ubv32_to_real : forall m, x.
 988     t'real (of_ubv32 m x) = round m (FromInt.from_int (BV32.t'int x))
 989 end
 990
 991 module Float64_BV_Converter
 992   use Float64
 993
 994   clone export Float_BV_Converter with
 995     type t = t,
 996     predicate is_finite = t'isFinite,
 997     predicate le = (.<=),
 998     function to_real = t'real,
 999     function round = round,
1000     axiom . (* TODO: "lemma"? "goal"? *)
1001
1002   axiom of_ubv32_to_real : forall m, x.
1003     t'real (of_ubv32 m x) = FromInt.from_int (BV32.t'int x)
1004 end