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
  (* exp_bias = 2^(eb - 1) - 1 *)
65
66
  (* 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
  (** Minimum and Maximum
Guillaume Melquiond's avatar
Guillaume Melquiond committed
97
98

     Note that we have to follow IEEE-754 and SMTLIB here. Two things to
MARCHE Claude's avatar
MARCHE Claude committed
99
     note in particular:
Guillaume Melquiond's avatar
Guillaume Melquiond committed
100

MARCHE Claude's avatar
MARCHE Claude committed
101
     1) min(-0, 0) is either 0 or -0, there is a choice
102

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

105
106
107
108
     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
109
   *)
110

MARCHE Claude's avatar
MARCHE Claude committed
111
  (** {3 Comparisons} *)
112

113
114
115
116
117
  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
118
  (** equality on floats, different from = since not (eq NaN NaN) *)
119

120
121
122
123
124
  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
125
126
  (** Notations *)

127

128
  (** {3 Classification of numbers} *)
129
130
131

  predicate is_normal    t
  predicate is_subnormal t
132
133
134
135
136
  val predicate is_zero      t
  val predicate is_infinite  t
  val predicate is_nan       t
  val predicate is_positive  t
  val predicate is_negative  t
137

138
  (** helper predicate for zeros, normals and subnormals. not defined
139
140
141
142
143
144
145
146
147
148
     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

149
150
151
152
  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)
153

MARCHE Claude's avatar
MARCHE Claude committed
154
155


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

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

166
  (** {3 Conversions to other sorts} *)
167
168
169
170
171
172

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

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
181
  axiom zero_to_real : forall x [is_zero x].
    is_zero x <-> is_finite x /\ to_real x = 0.0
182

183
  (** {3 Conversions from other sorts} *)
184
185
186
187

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

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

196
  (** {3 Conversions to other sorts} *)
197

198
199
  (* Intended semantics for to_int are the same as (_ fp.to_sbv) with a    *)
  (* suitably sized bitvector. Safe minimum sizes are given below:         *)
200
201
202
203
204
205
206
207
208
209
210
211
  (*    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

212
  (** {3 Arithmetic} *)
213

214
215
  (* The intended meaning for round is the rounding for floating point as  *)
  (* described on p17 of IEEE-754. For results where the corresponding     *)
216
217
218
219
220
221
222
223
  (* 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 *)
224
225
  constant max_int  : int

226
  constant emax : int = pow2 (eb - 1)
Clément Fumex's avatar
Clément Fumex committed
227

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

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

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
238
239
240
241

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

  (* Axioms on round *)

242
  axiom Bounded_real_no_overflow [@W:non_conservative_extension:N] :
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
    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

264
265
266
  (* The biggest representable integer whose predecessor (i.e. -1) is  representable *)
  constant pow2sb : int (* defined when cloning *)
  axiom pow2sb: pow2sb = pow2 sb
Clément Fumex's avatar
Clément Fumex committed
267

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

  (* round and integers *)

  axiom Exact_rounding_for_integers:
    forall m:mode, i:int.
275
      in_safe_int_range i ->
276
277
278
279
280
281
282
283
284
285
286
287
        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.
288
    is_finite x -> is_finite y -> not (is_zero x) -> x .= y -> x = y
289
290
291
292
293
294

  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

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

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

  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
308
309
    /\ ((is_finite x /\ is_finite y)
        \/ (is_infinite x /\ is_infinite y /\ same_sign x y)))
310

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

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

  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

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

337
  axiom lt_special: forall x y [lt x y]. lt x y ->
338
339
340
341
      ((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)))

Clément Fumex's avatar
Clément Fumex committed
342
343
  axiom lt_lt_finite: forall x y z. lt x y -> lt y z -> is_finite y

344
345
346
347
348
349
350
351
352
353
354
  (*  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

355
356
357
358
  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
359
360
361
362
363
364
365
366
367
368
369
370
371
372

  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

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

376
  (** {3 Overflow} *)
377

378
379
380
  (* 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
     *)
381
382
383
384
385
386
387
388
389
390
391
  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

392
393
  (* This predicate is used to tell what is the sign of zero in the
     axioms specifying add and sub *)
394
395
396
397
398
399
400
401
402
403
  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].
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
409
  lemma add_finite_rev: forall m:mode, x y:t [add m x y].
    is_finite (add m x y) ->
410
    is_finite x /\ is_finite y
411

412
413
414
415
416
417
  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)

418
419
  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) ->
420
    is_finite (sub m x y) /\
421
    to_real (sub m x y) = round m (to_real x -. to_real y)
422

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

427
428
429
430
431
432
  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)

433
434
  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) ->
435
436
437
    is_finite (mul m x y) /\
    to_real (mul m x y) = round m (to_real x *. to_real y)

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

442
443
444
445
446
447
  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)

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

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

459
460
461
462
463
464
  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)

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

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

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

481
482
483
484
485
  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)

486
487
  axiom abs_universal : forall x:t [abs x]. not (is_negative (abs x))

488
489
490
491
492
  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)
493

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

498
499
500
501
502
503
  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)

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
512
513
  lemma sqrt_finite_rev: forall m:mode, x:t [sqrt m x].
    is_finite (sqrt m x) ->
    is_finite x /\ to_real x >=. 0. /\
514
    to_real (sqrt m x) = round m (S.sqrt (to_real x))
515

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

519
  axiom add_special: forall m:mode, x y:t [add m x y].
520
521
522
523
524
525
526
527
528
529
530
531
532
533
    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)
534
535
536
    /\
    (is_finite x /\ is_finite y
      -> 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
540
541
542
543
544
545
546
547
548
549
550
551
    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)
552
      -> same_sign_real r (to_real x -. to_real y) /\ overflow_value m r)
553
554
555
    /\
    (is_finite x /\ is_finite y
      -> 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
559
560
561
562
563
564
565
566
567
568
    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)
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
573
574
575
576
577
    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) /\
578
579
        not no_overflow m (to_real x /. to_real y)
         -> 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
584
585
586
587
588
    /\ (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))
589

590
591
592
593
  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))
594

595
596
  axiom fma_special: forall m:mode, x y z:t [fma m x y z].
    let r = fma m x y z in
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
       (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)
621
622
623
624
625
626
627
628
629
    /\ (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)
630
631
    /\ (is_minus_infinity x -> is_nan r)
    /\ (is_finite x /\ to_real x <. 0.0 -> is_nan r)
632
633
    /\ (is_zero x -> same_sign r x)
    /\ (is_finite x /\ to_real x >. 0.0 -> is_positive r)
634

Clément Fumex's avatar
Clément Fumex committed
635
636
637
  (* exact arithmetic with integers *)

  axiom of_int_add_exact: forall m n, i j.
638
639
    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))
Clément Fumex's avatar
Clément Fumex committed
640
641

  axiom of_int_sub_exact: forall m n, i j.
642
643
    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))
Clément Fumex's avatar
Clément Fumex committed
644
645

  axiom of_int_mul_exact: forall m n, i j.
646
647
    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))
Clément Fumex's avatar
Clément Fumex committed
648

649
  (* min and max *)
650
651
652
653
654
655
656
657

  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

  (* _____________ *)

658
  use real.Truncate as Truncate
659

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

663
664
665
  (** characterizing integers *)

  (* by construction *)
Clément Fumex's avatar
Clément Fumex committed
666
667
  axiom zeroF_is_int: is_int zeroF

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

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

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

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
713
714
  (* 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 ...*)

715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
  (** 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

Clément Fumex's avatar
Clément Fumex committed
734
735
  axiom ceil_to_real: forall x:t.
    is_finite x ->
736
737
738
      to_real (roundToIntegral RTP x) = FromInt.from_int (Truncate.ceil (to_real x))

  axiom ceil_to_int: forall m:mode, x:t.
Clément Fumex's avatar
Clément Fumex committed
739
    is_finite x ->
740
741
742
743
744
745
746
747
      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)

Clément Fumex's avatar
Clément Fumex committed
748
749
  axiom floor_to_real: forall x:t.
    is_finite x ->
750
751
752
      to_real (roundToIntegral RTN x) = FromInt.from_int (Truncate.floor (to_real x))

  axiom floor_to_int: forall m:mode, x:t.
Clément Fumex's avatar
Clément Fumex committed
753
    is_finite x ->
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
      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.
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
784
785
786
787
      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

Clément Fumex's avatar
Clément Fumex committed
788
789
790
  axiom neg_to_int: forall m x.
    is_int x -> to_int m (neg x) = - (to_int m x)

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

795
796
(** {2 Conversions to/from bitvectors} *)

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

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

807
808
809
810
  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
811

812
813
814
815
  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
816

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

  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.
843
    to_real (of_ubv8  m x) = FromInt.from_int (BV8.t'int x)
844
  axiom of_ubv16_to_real: forall m, x.
845
    to_real (of_ubv16 m x) = FromInt.from_int (BV16.t'int x)
846
847
  (* of_ubv32_to_real is defined at cloning *)
  axiom of_ubv64_to_real: forall m, x.
848
    to_real (of_ubv64 m x) = round m (FromInt.from_int (BV64.t'int x))
849
850
end

851
852
(** {2 Standard simple precision floats (32 bits)} *)

853
module Float32
854
  use int.Int as Int
855
  use real.Real
856

857
  type t = < float 8 24 >
858

859
  constant pow2sb : int = 16777216
860
861
862
863
  constant max_real : real = 0x1.FFFFFEp127

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

876
877
  lemma round_bound_ne :
    forall x:real [round RNE x].
Clément Fumex's avatar
Clément Fumex committed
878
879
      no_overflow RNE x ->
        x - 0x1p-24 * Abs.abs(x) - 0x1p-150 <= round RNE x <= x + 0x1p-24 * Abs.abs(x) + 0x1p-150
880
881

  lemma round_bound :
882
    forall m:mode, x:real [round m x].
Clément Fumex's avatar
Clément Fumex committed
883
      no_overflow m x ->
884
      x - 0x1p-23 * Abs.abs(x) - 0x1p-149 <= round m x <= x + 0x1p-23 * Abs.abs(x) + 0x1p-149
885
886
end

887
888
(** {2 Standard double precision floats (64 bits)} *)

889
module Float64
890
  use int.Int as Int
891
  use real.Real
892

893
  type t = < float 11 53 >
894

895
  constant pow2sb : int = 9007199254740992
896
897
898
899
  constant max_real : real = 0x1.FFFFFFFFFFFFFp1023

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

912
913
  lemma round_bound_ne :
    forall x:real [round RNE x].
Clément Fumex's avatar
Clément Fumex committed
914
915
      no_overflow RNE x ->
        x - 0x1p-53 * Abs.abs(x) - 0x1p-1075 <= round RNE x <= x + 0x1p-53 * Abs.abs(x) + 0x1p-1075
916
917

  lemma round_bound :
918
    forall m:mode, x:real [round m x].
Clément Fumex's avatar
Clément Fumex committed
919
      no_overflow m x ->
920
      x - 0x1p-52 * Abs.abs(x) - 0x1p-1074 <= round m x <= x + 0x1p-52 * Abs.abs(x) + 0x1p-1074
921
922
end

923
924
(** {2 Conversions between float formats} *)

925
module FloatConverter
926

927
928
  use Float64 as Float64
  use Float32 as Float32
929
930
931
932
933
934
935
936
937
938
939

  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:
940
941
942
    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
943
944

  lemma to_float32_conv:
945
946
947
948
    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)
949
950
951

end

952
module Float32_BV_Converter
953
  use Float32
954
955
956

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

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

967
module Float64_BV_Converter
968
  use Float64
969
970
971

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

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