ieee_float.mlw 32.1 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
29
  (** {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
     </ul} *)
30
31

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

34
module GenericFloat
35
36

  use import int.Int
37
  use import bv.Pow2int
38
39
40
41
42
43
  use real.Abs
  use real.FromInt
  use real.Truncate
  use import real.RealInfix
  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
  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

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

MARCHE Claude's avatar
MARCHE Claude committed
79
80
81
82
  function abs         t   : t   (** Absolute value *)
  function neg         t   : t   (** Opposite *)
  function fma  mode t t t : t   (** Fused multiply-add: x * y + z *)
  function sqrt mode   t   : t   (** Square root *)
83
84
85
86
87
88

  function (.-_) (x:t)   : t = neg x
  function (.+)  (x y:t) : t = add RNE x y
  function (.-)  (x y:t) : t = sub RNE x y
  function (.*)  (x y:t) : t = mul RNE x y
  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

  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

MARCHE Claude's avatar
MARCHE Claude committed
103
   *)
104

MARCHE Claude's avatar
MARCHE Claude committed
105
  (** {3 Comparisons} *)
106

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

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

121

122
  (** {3 Classification of numbers} *)
123
124
125

  predicate is_normal    t
  predicate is_subnormal t
126
  predicate is_zero      t
127
128
129
130
131
  predicate is_infinite  t
  predicate is_nan       t
  predicate is_positive  t
  predicate is_negative  t

132
  (** helper predicate for zeros, normals and subnormals. not defined
133
134
135
136
137
138
139
140
141
142
     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

143
144
145
146
  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)
147

MARCHE Claude's avatar
MARCHE Claude committed
148
149


150
  (** {3 Conversions from other sorts} *)
151
152
153
154
155
156
157
158
159

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

160
  (** {3 Conversions to other sorts} *)
161
162
163
164
165
166

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

167
  (** {2 Part II - Private Axiomatisation} *)
168

169
  (** {3 Constructors and Constants} *)
170

171
  axiom zeroF_is_positive : is_positive zeroF
172
  axiom zeroF_is_zero     : is_zero zeroF
173

174
175
  axiom zero_to_real : forall x [is_zero x].
    is_zero x <-> is_finite x /\ to_real x = 0.0
176

177
  (** {3 Conversions from other sorts} *)
178
179
180
181

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

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

190
  (** {3 Conversions to other sorts} *)
191

192
193
  (* Intended semantics for to_int are the same as (_ fp.to_sbv) with a    *)
  (* suitably sized bitvector. Safe minimum sizes are given below:         *)
194
195
196
197
198
199
200
201
202
203
204
205
  (*    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

206
  (** {3 Arithmetic} *)
207

208
209
  (* The intended meaning for round is the rounding for floating point as  *)
  (* described on p17 of IEEE-754. For results where the corresponding     *)
210
211
212
213
214
215
216
217
  (* 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 *)
218
219
  constant max_int  : int

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

222
  axiom max_int : max_int = pow2 emax - pow2 (emax - sb)
223
  axiom max_real_int: max_real = FromInt.from_int max_int
224
225
226

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

227
  predicate in_int_range (i:int) = - max_int <= i <= max_int
228

229
  axiom is_finite: forall x:t. is_finite x -> in_range (to_real x)
230
231
232
233
234
235

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

  (* Axioms on round *)

236
  axiom Bounded_real_no_overflow [@W:non_conservative_extension:N] :
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
    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

258
259
260
  (* 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
261

262
263
  (* range in which every integer is representable *)
  predicate in_safe_int_range (i: int) = - pow2sb <= i <= pow2sb
264
265
266
267
268

  (* round and integers *)

  axiom Exact_rounding_for_integers:
    forall m:mode, i:int.
269
      in_safe_int_range i ->
270
271
272
273
274
275
276
277
278
279
280
281
        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.
282
    is_finite x -> is_finite y -> not (is_zero x) -> x .= y -> x = y
283
284
285
286
287
288

  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

289
  axiom eq_sym :
MARCHE Claude's avatar
MARCHE Claude committed
290
    forall x y. x .= y -> y .= x
291

292
  axiom eq_trans :
MARCHE Claude's avatar
MARCHE Claude committed
293
    forall x y z. x .= y -> y .= z -> x .= z
294
295
296
297
298
299
300
301

  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
302
303
    /\ ((is_finite x /\ is_finite y)
        \/ (is_infinite x /\ is_infinite y /\ same_sign x y)))
304

305
306
  axiom lt_finite: forall x y [lt x y].
    is_finite x /\ is_finite y -> (lt x y <-> to_real x <. to_real y)
307

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

  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

326
  axiom le_special: forall x y [le x y]. le x y ->
327
328
329
330
      ((is_finite x /\ is_finite y)
   \/ ((is_minus_infinity x /\ is_not_nan y)
   \/  (is_not_nan x /\ is_plus_infinity y)))

331
  axiom lt_special: forall x y [lt x y]. lt x y ->
332
333
334
335
      ((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
336
337
  axiom lt_lt_finite: forall x y z. lt x y -> lt y z -> is_finite y

338
339
340
341
342
343
344
345
346
347
348
  (*  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

349
350
351
352
  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
353
354
355
356
357
358
359
360
361
362
363
364
365
366

  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

367
368
369
  predicate product_sign (z x y: t) =
    (same_sign x y -> is_positive z) /\ (diff_sign x y -> is_negative z)

370
  (** {3 Overflow} *)
371

372
373
374
  (* 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
     *)
375
376
377
378
379
380
381
382
383
384
385
  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

386
387
  (* This predicate is used to tell what is the sign of zero in the
     axioms specifying add and sub *)
388
389
390
391
392
393
394
395
396
397
  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].
398
    is_finite x -> is_finite y -> no_overflow m (to_real x +. to_real y) ->
399
    is_finite (add m x y) /\
400
    to_real (add m x y) = round m (to_real x +. to_real y)
401

402
403
  lemma add_finite_rev: forall m:mode, x y:t [add m x y].
    is_finite (add m x y) ->
404
    is_finite x /\ is_finite y
405

406
407
408
409
410
411
  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)

412
413
  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) ->
414
    is_finite (sub m x y) /\
415
    to_real (sub m x y) = round m (to_real x -. to_real y)
416

417
418
  lemma sub_finite_rev: forall m:mode, x y:t [sub m x y].
    is_finite (sub m x y) ->
419
    is_finite x /\ is_finite y
420

421
422
423
424
425
426
  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)

427
428
  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) ->
429
430
431
    is_finite (mul m x y) /\
    to_real (mul m x y) = round m (to_real x *. to_real y)

432
433
  lemma mul_finite_rev: forall m:mode, x y:t [mul m x y].
    is_finite (mul m x y) ->
434
    is_finite x /\ is_finite y
435

436
437
438
439
440
441
  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)

442
443
  axiom div_finite: forall m:mode, x y:t [div m x y].
    is_finite x -> is_finite y ->
444
    not is_zero y -> no_overflow m (to_real x /. to_real y) ->
445
446
447
    is_finite (div m x y) /\
    to_real (div m x y) = round m (to_real x /. to_real y)

448
449
  lemma div_finite_rev: forall m:mode, x y:t [div m x y].
    is_finite (div m x y) ->
450
    (is_finite x /\ is_finite y /\ not is_zero y) \/
451
452
    (is_finite x /\ is_infinite y /\ to_real (div m x y) = 0.)

453
454
455
456
457
458
  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)

459
460
461
  axiom neg_finite: forall x:t [neg x].
    is_finite x ->
    is_finite (neg x) /\
462
463
    to_real (neg x) = -. to_real x

464
465
466
467
468
  lemma neg_finite_rev: forall x:t [neg x].
    is_finite (neg x) ->
    is_finite x /\
    to_real (neg x) = -. to_real x

469
470
  axiom abs_finite: forall x:t [abs x].
    is_finite x ->
471
472
473
474
    is_finite (abs x) /\
    to_real (abs x) = Abs.abs (to_real x) /\
    is_positive (abs x)

475
476
477
478
479
  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)

480
481
  axiom abs_universal : forall x:t [abs x]. not (is_negative (abs x))

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

488
489
  lemma fma_finite_rev: forall m:mode, x y z:t [fma m x y z].
    is_finite (fma m x y z) ->
490
    is_finite x /\ is_finite y /\ is_finite z
491

492
493
494
495
496
497
  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)

498
499
  use real.Square

500
  axiom sqrt_finite: forall m:mode, x:t [sqrt m x].
501
    is_finite x -> to_real x >=. 0. ->
502
503
    is_finite (sqrt m x) /\
    to_real (sqrt m x) = round m (Square.sqrt (to_real x))
504

505
506
507
508
509
  lemma sqrt_finite_rev: forall m:mode, x:t [sqrt m x].
    is_finite (sqrt m x) ->
    is_finite x /\ to_real x >=. 0. /\
    to_real (sqrt m x) = round m (Square.sqrt (to_real x))

510
511
512
  predicate same_sign_real (x:t) (r:real) =
    (is_positive x /\ r >. 0.0) \/ (is_negative x /\ r <. 0.0)

513
  axiom add_special: forall m:mode, x y:t [add m x y].
514
515
516
517
518
519
520
521
522
523
524
525
526
527
    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)
528
529
530
    /\
    (is_finite x /\ is_finite y
      -> if same_sign x y then same_sign r x else sign_zero_result m r)
531

532
  axiom sub_special: forall m:mode, x y:t [sub m x y].
533
534
535
536
537
538
539
540
541
542
543
544
545
    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)
546
      -> same_sign_real r (to_real x -. to_real y) /\ overflow_value m r)
547
548
549
    /\
    (is_finite x /\ is_finite y
      -> if diff_sign x y then same_sign r x else sign_zero_result m r)
550

551
  axiom mul_special: forall m:mode, x y:t [mul m x y].
552
553
554
555
556
557
558
559
560
561
562
    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)
563
    /\ (not is_nan r -> product_sign r x y)
564

565
  axiom div_special: forall m:mode, x y:t [div m x y].
566
567
568
569
570
571
    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) /\
572
573
        not no_overflow m (to_real x /. to_real y)
         -> overflow_value m r)
574
    /\ (is_finite x /\ is_zero y /\ not (is_zero x)
575
         -> is_infinite r)
576
    /\ (is_zero x /\ is_zero y -> is_nan r)
577
578
579
580
581
582
    /\ (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))
583

584
585
586
587
  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))
588

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

Clément Fumex's avatar
Clément Fumex committed
629
630
631
  (* exact arithmetic with integers *)

  axiom of_int_add_exact: forall m n, i j.
632
633
    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
634
635

  axiom of_int_sub_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)) (sub n (of_int m i) (of_int m j))
Clément Fumex's avatar
Clément Fumex committed
638
639

  axiom of_int_mul_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)) (mul n (of_int m i) (of_int m j))
Clément Fumex's avatar
Clément Fumex committed
642

643
  (* min and max *)
644
645
646
647
648
649
650
651
652
653

  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

  (* _____________ *)

  use real.Truncate

654
  (* This predicate specify if a float is finite and is an integer *)
655
656
  predicate is_int (x:t)

657
658
659
  (** characterizing integers *)

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

662
  axiom of_int_is_int: forall m, x.
663
    in_int_range x -> is_int (of_int m x)
664

665
666
  axiom big_float_is_int: forall m i.
    is_finite i ->
667
    i .<= neg (of_int m pow2sb) \/ (of_int m pow2sb) .<= i ->
668
      is_int i
669
670
671
672

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

673
674
675
676
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
  (* 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 ...*)

709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
  (** 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
728
729
  axiom ceil_to_real: forall x:t.
    is_finite x ->
730
731
732
      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
733
    is_finite x ->
734
735
736
737
738
739
740
741
      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
742
743
  axiom floor_to_real: forall x:t.
    is_finite x ->
744
745
746
      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
747
    is_finite x ->
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
      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.
773
    is_finite x -> is_finite y -> le x y -> to_int m x <= to_int m y
774
775

  axiom to_int_of_int: forall m:mode, i:int.
776
    in_safe_int_range i ->
777
778
779
780
781
      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
782
783
784
  axiom neg_to_int: forall m x.
    is_int x -> to_int m (neg x) = - (to_int m x)

785
786
787
788
  axiom roundToIntegral_is_finite  : forall m:mode, x:t. is_finite x ->
    is_finite (roundToIntegral m x)
end

789
790
(** {2 Conversions to/from bitvectors} *)

791
module Float_BV_Converter
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
  use bv.BV8
  use bv.BV16
  use bv.BV32
  use bv.BV64
  use import RoundingMode

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

  function of_ubv8  mode BV8.t  : t
  function of_ubv16 mode BV16.t : t
  function of_ubv32 mode BV32.t : t
  function of_ubv64 mode BV64.t : t

  function to_ubv8  mode t : BV8.t
  function to_ubv16 mode t : BV16.t
  function to_ubv32 mode t : BV32.t
  function to_ubv64 mode t : BV64.t

  use import real.RealInfix
  use real.FromInt

  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.
837
    to_real (of_ubv8  m x) = FromInt.from_int (BV8.t'int x)
838
  axiom of_ubv16_to_real: forall m, x.
839
    to_real (of_ubv16 m x) = FromInt.from_int (BV16.t'int x)
840
841
  (* of_ubv32_to_real is defined at cloning *)
  axiom of_ubv64_to_real: forall m, x.
842
    to_real (of_ubv64 m x) = round m (FromInt.from_int (BV64.t'int x))
843
844
end

845
846
(** {2 Standard simple precision floats (32 bits)} *)

847
module Float32
848
849
850
  use int.Int
  use import real.Real

851
  type t = < float 8 24 >
852

853
  constant pow2sb : int = 16777216
854
855
856
857
  constant max_real : real = 0x1.FFFFFEp127

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

870
871
  lemma round_bound_ne :
    forall x:real [round RNE x].
Clément Fumex's avatar
Clément Fumex committed
872
873
      no_overflow RNE x ->
        x - 0x1p-24 * Abs.abs(x) - 0x1p-150 <= round RNE x <= x + 0x1p-24 * Abs.abs(x) + 0x1p-150
874
875

  lemma round_bound :
876
    forall m:mode, x:real [round m x].
Clément Fumex's avatar
Clément Fumex committed
877
      no_overflow m x ->
878
      x - 0x1p-23 * Abs.abs(x) - 0x1p-149 <= round m x <= x + 0x1p-23 * Abs.abs(x) + 0x1p-149
879
880
end

881
882
(** {2 Standard double precision floats (64 bits)} *)

883
module Float64
884
885
886
  use int.Int
  use import real.Real

887
  type t = < float 11 53 >
888

889
  constant pow2sb : int = 9007199254740992
890
891
892
893
  constant max_real : real = 0x1.FFFFFFFFFFFFFp1023

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

906
907
  lemma round_bound_ne :
    forall x:real [round RNE x].
Clément Fumex's avatar
Clément Fumex committed
908
909
      no_overflow RNE x ->
        x - 0x1p-53 * Abs.abs(x) - 0x1p-1075 <= round RNE x <= x + 0x1p-53 * Abs.abs(x) + 0x1p-1075
910
911

  lemma round_bound :
912
    forall m:mode, x:real [round m x].
Clément Fumex's avatar
Clément Fumex committed
913
      no_overflow m x ->
914
      x - 0x1p-52 * Abs.abs(x) - 0x1p-1074 <= round m x <= x + 0x1p-52 * Abs.abs(x) + 0x1p-1074
915
916
end

917
918
(** {2 Conversions between float formats} *)

919
module FloatConverter
920
921
922
923
924
925
926
927
928
929
930
931
932
933

  use Float64
  use Float32

  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:
934
935
936
    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
937
938

  lemma to_float32_conv:
939
940
941
942
    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)
943
944
945

end

946
module Float32_BV_Converter
947
948
949
950
  use import Float32

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

  axiom of_ubv32_to_real : forall m, x.
958
    t'real (of_ubv32 m x) = round m (FromInt.from_int (BV32.t'int x))
959
960
end

961
module Float64_BV_Converter
962
963
964
965
  use import Float64

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

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