floating_point.mlw 21.4 KB
Newer Older
MARCHE Claude's avatar
MARCHE Claude committed
1

2
(** {1 Formalization of Floating-Point Arithmetic}
MARCHE Claude's avatar
MARCHE Claude committed
3

MARCHE Claude's avatar
MARCHE Claude committed
4 5 6
This formalization follows the {h <a
href="http://ieeexplore.ieee.org/servlet/opac?punumber=4610933">IEEE-754
standard</a>}.
MARCHE Claude's avatar
MARCHE Claude committed
7 8

*)
MARCHE Claude's avatar
MARCHE Claude committed
9

10 11
(** {2 Definition of IEEE-754 rounding modes} *)

12
module Rounding
MARCHE Claude's avatar
MARCHE Claude committed
13

14
  type mode = NearestTiesToEven | ToZero | Up | Down | NearestTiesToAway
MARCHE Claude's avatar
MARCHE Claude committed
15
  (** nearest ties to even, to zero, upward, downward, nearest ties to away *)
MARCHE Claude's avatar
MARCHE Claude committed
16 17 18

end

19
(** {2 Handling of IEEE-754 special values}
MARCHE Claude's avatar
MARCHE Claude committed
20

21
Special values are +infinity, -infinity, NaN, +0, -0.  These are
22
handled as described in {h <a
MARCHE Claude's avatar
MARCHE Claude committed
23 24
href="http://www.lri.fr/~marche/ayad10ijcar.pdf">[Ayad, Marché, IJCAR,
2010]</a>}.
25 26

*)
27

28
module SpecialValues
Andrei Paskevich's avatar
Andrei Paskevich committed
29 30

  type class = Finite | Infinite | NaN
MARCHE Claude's avatar
MARCHE Claude committed
31 32 33

  type sign = Neg | Pos

34
  use real.Real
MARCHE Claude's avatar
MARCHE Claude committed
35

36
  inductive same_sign_real sign real =
37 38
    | Neg_case: forall x:real. x < 0.0 -> same_sign_real Neg x
    | Pos_case: forall x:real. x > 0.0 -> same_sign_real Pos x
MARCHE Claude's avatar
MARCHE Claude committed
39

40
(*** useful ???
MARCHE Claude's avatar
MARCHE Claude committed
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71

lemma sign_not_pos_neg : forall x:t.
      sign(x) <> Pos -> sign(x) = Neg
lemma sign_not_neg_pos : forall x:t.
      sign(x) <> Neg -> sign(x) = Pos

lemma sign_not_pos_neg : forall x:double.
      sign(x) <> Pos -> sign(x) = Neg
lemma sign_not_neg_pos : forall x:double.
      sign(x) <> Neg -> sign(x) = Pos
*)

lemma same_sign_real_zero1 :
  forall b:sign. not same_sign_real b 0.0

lemma same_sign_real_zero2 :
  forall x:real.
    same_sign_real Neg x /\ same_sign_real Pos x -> false

lemma same_sign_real_zero3 :
  forall b:sign, x:real. same_sign_real b x -> x <> 0.0

lemma same_sign_real_correct2 :
  forall b:sign, x:real.
    same_sign_real b x -> (x < 0.0 <-> b = Neg)

lemma same_sign_real_correct3 :
  forall b:sign, x:real.
    same_sign_real b x -> (x > 0.0 <-> b = Pos)


MARCHE Claude's avatar
MARCHE Claude committed
72 73
end

74 75
(** {2 Generic theory of floats}

76 77
The theory is parametrized by the real constant `max` which denotes
the maximal representable number, the real constant `min` which denotes
MARCHE Claude's avatar
MARCHE Claude committed
78
the minimal number whose rounding is not null, and the integer constant
79 80
`max_representable_integer` which is the maximal integer `n` such that
every integer between `0` and `n` are representable.
81

82
*)
83

84
module GenFloat
MARCHE Claude's avatar
MARCHE Claude committed
85

86 87 88 89
  use Rounding
  use real.Real
  use real.Abs
  use real.FromInt
90
  use int.Int as Int
MARCHE Claude's avatar
MARCHE Claude committed
91 92 93

  type t

Andrei Paskevich's avatar
Andrei Paskevich committed
94
  function round mode real : real
MARCHE Claude's avatar
MARCHE Claude committed
95

Andrei Paskevich's avatar
Andrei Paskevich committed
96 97 98
  function value t : real
  function exact t : real
  function model t : real
MARCHE Claude's avatar
MARCHE Claude committed
99

Andrei Paskevich's avatar
Andrei Paskevich committed
100 101
  function round_error (x : t) : real = abs (value x - exact x)
  function total_error (x : t) : real = abs (value x - model x)
MARCHE Claude's avatar
MARCHE Claude committed
102

103
  constant max : real
MARCHE Claude's avatar
MARCHE Claude committed
104

Andrei Paskevich's avatar
Andrei Paskevich committed
105
  predicate no_overflow (m:mode) (x:real) = abs (round m x) <= max
MARCHE Claude's avatar
MARCHE Claude committed
106

107
  axiom Bounded_real_no_overflow :
108
    forall m:mode, x:real. abs x <= max -> no_overflow m x
MARCHE Claude's avatar
MARCHE Claude committed
109 110

  axiom Round_monotonic :
111
    forall m:mode, x y:real. x <= y -> round m x <= round m y
MARCHE Claude's avatar
MARCHE Claude committed
112

113 114 115 116 117 118 119 120 121
  axiom Round_idempotent :
    forall m1 m2:mode, x:real. round m1 (round m2 x) = round m2 x

  axiom Round_value :
    forall m:mode, x:t. round m (value x) = value x

  axiom Bounded_value :
    forall x:t. abs (value x) <= max

122
  constant max_representable_integer : int
MARCHE Claude's avatar
MARCHE Claude committed
123 124 125

  axiom Exact_rounding_for_integers:
    forall m:mode,i:int.
Andrei Paskevich's avatar
Andrei Paskevich committed
126
      Int.(<=) (Int.(-_) max_representable_integer) i /\
127 128
      Int.(<=) i max_representable_integer ->
        round m (from_int i) = from_int i
MARCHE Claude's avatar
MARCHE Claude committed
129

MARCHE Claude's avatar
MARCHE Claude committed
130
  (** rounding up and down *)
131

Andrei Paskevich's avatar
Andrei Paskevich committed
132
  axiom Round_down_le:
133
    forall x:real. round Down x <= x
Andrei Paskevich's avatar
Andrei Paskevich committed
134
  axiom Round_up_ge:
135
    forall x:real. round Up x >= x
Andrei Paskevich's avatar
Andrei Paskevich committed
136
  axiom Round_down_neg:
137 138 139
    forall x:real. round Down (-x) = - round Up x
  axiom Round_up_neg:
    forall x:real. round Up (-x) = - round Down x
MARCHE Claude's avatar
MARCHE Claude committed
140

141
  (** rounding into a float instead of a real *)
142

143 144 145 146 147 148
  function round_logic mode real : t
  axiom Round_logic_def:
    forall m:mode, x:real.
      no_overflow m x ->
      value (round_logic m x) = round m x

MARCHE Claude's avatar
MARCHE Claude committed
149 150
end

151 152 153 154 155 156 157

(** {2 Format of Single precision floats}

A.k.a. binary32 numbers.

*)

158
module SingleFormat
159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179

  type single

  constant max_single : real = 0x1.FFFFFEp127
  constant max_int : int = 16777216 (* 2^24 *)

(*
  clone export GenFloatSpecStrict with
    type t = single,
    constant max = max_single,
    constant max_representable_integer = max_int
*)

end

(** {2 Format of Double precision floats}

A.k.a. binary64 numbers.

*)

180
module DoubleFormat
181 182 183 184 185 186 187 188 189 190 191 192 193 194 195

  type double

  constant max_double : real = 0x1.FFFFFFFFFFFFFp1023
  constant max_int : int = 9007199254740992 (* 2^53 *)

(*
  clone export GenFloatSpecStrict with
    type t = double,
    constant max = max_double,
    constant max_representable_integer = max_int
*)

end

196
module GenFloatSpecStrict
197

198 199
  use Rounding
  use real.Real
200
  clone export GenFloat with axiom .
201

202 203 204 205

  predicate of_real_post (m:mode) (x:real) (res:t) =
    value res = round m x /\ exact res = x /\ model res = x

206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234
  (** {3 binary operations} *)

  predicate add_post (m:mode) (x y res:t) =
    value res = round m (value x + value y) /\
    exact res = exact x + exact y /\
    model res = model x + model y

  predicate sub_post (m:mode) (x y res:t) =
    value res = round m (value x - value y) /\
    exact res = exact x - exact y /\
    model res = model x - model y

  predicate mul_post (m:mode) (x y res:t) =
    value res = round m (value x * value y) /\
    exact res = exact x * exact y /\
    model res = model x * model y

  predicate div_post (m:mode) (x y res:t) =
    value res = round m (value x / value y) /\
    exact res = exact x / exact y /\
    model res = model x / model y

  predicate neg_post (x res:t) =
    value res = - value x /\
    exact res = - exact x /\
    model res = - model x

  (** {3 Comparisons} *)

235 236 237 238
  predicate lt (x:t) (y:t) = value x < value y

  predicate gt (x:t) (y:t) = value x > value y

239 240 241

end

242

243
module Single
MARCHE Claude's avatar
MARCHE Claude committed
244

245
  use export SingleFormat
MARCHE Claude's avatar
MARCHE Claude committed
246

247
  clone export GenFloatSpecStrict with
MARCHE Claude's avatar
MARCHE Claude committed
248
    type t = single,
249
    constant max = max_single,
250
    constant max_representable_integer = max_int,
251 252
    lemma Bounded_real_no_overflow,
    axiom . (* TODO: "lemma"? "goal"? *)
MARCHE Claude's avatar
MARCHE Claude committed
253 254 255

end

256

257
module Double
MARCHE Claude's avatar
MARCHE Claude committed
258

259
  use export DoubleFormat
MARCHE Claude's avatar
MARCHE Claude committed
260

261
  clone export GenFloatSpecStrict with
MARCHE Claude's avatar
MARCHE Claude committed
262
    type t = double,
263
    constant max = max_double,
264
    constant max_representable_integer = max_int,
265 266
    lemma Bounded_real_no_overflow,
    axiom . (* TODO: "lemma"? "goal"? *)
MARCHE Claude's avatar
MARCHE Claude committed
267 268 269

end

270 271


272 273
(** {2 Generic Full theory of floats}

274
This theory extends the generic theory above by adding the special values.
275

276
*)
MARCHE Claude's avatar
MARCHE Claude committed
277

278
module GenFloatFull
MARCHE Claude's avatar
MARCHE Claude committed
279

280 281 282
  use SpecialValues
  use Rounding
  use real.Real
MARCHE Claude's avatar
MARCHE Claude committed
283

284
  clone export GenFloat with axiom .
MARCHE Claude's avatar
MARCHE Claude committed
285

286
  (** {3 special values} *)
287

Andrei Paskevich's avatar
Andrei Paskevich committed
288
  function class t : class
MARCHE Claude's avatar
MARCHE Claude committed
289

MARCHE Claude's avatar
MARCHE Claude committed
290 291 292 293
  predicate is_finite (x:t) = class x = Finite
  predicate is_infinite (x:t) = class x = Infinite
  predicate is_NaN (x:t) = class x = NaN
  predicate is_not_NaN (x:t) = is_finite x \/ is_infinite x
294 295 296 297 298

  lemma is_not_NaN: forall x:t. is_not_NaN x <-> not (is_NaN x)

  function sign t  : sign

MARCHE Claude's avatar
MARCHE Claude committed
299 300
  predicate same_sign_real (x:t) (y:real) =
    SpecialValues.same_sign_real (sign x) y
301
  predicate same_sign (x:t) (y:t)  = sign x = sign y
302 303 304 305 306 307 308 309
  predicate diff_sign (x:t) (y:t)  = sign x <> sign y

  predicate sign_zero_result (m:mode) (x:t) =
    value x = 0.0 ->
    match m with
    | Down -> sign x = Neg
    | _ -> sign x = Pos
    end
310

MARCHE Claude's avatar
MARCHE Claude committed
311 312 313 314 315 316
  predicate is_minus_infinity (x:t) = is_infinite x /\ sign x = Neg
  predicate is_plus_infinity (x:t) = is_infinite x /\ sign x = Pos
  predicate is_gen_zero (x:t) = is_finite x /\ value x = 0.0
  predicate is_gen_zero_plus (x:t) = is_gen_zero x /\ sign x = Pos
  predicate is_gen_zero_minus (x:t) = is_gen_zero x /\ sign x = Neg

MARCHE Claude's avatar
MARCHE Claude committed
317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351
  (** {3 Useful lemmas on sign} *)

  (* non-zero finite gen_float has the same sign as its float_value *)
  lemma finite_sign :
    forall x:t.
      class x = Finite /\ value x <> 0.0 -> same_sign_real x (value x)

  lemma finite_sign_pos1:
    forall x:t.
      class x = Finite /\ value x > 0.0 -> sign x = Pos

  lemma finite_sign_pos2:
    forall x:t.
      class x = Finite /\ value x <> 0.0 /\ sign x = Pos -> value x > 0.0

  lemma finite_sign_neg1:
    forall x:t. class x = Finite /\ value x < 0.0 -> sign x = Neg

  lemma finite_sign_neg2:
    forall x:t.
      class x = Finite /\ value x <> 0.0 /\ sign x = Neg -> value x < 0.0

  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.
      class x = Finite /\ class y = Finite /\ value x * value y < 0.0
      -> diff_sign x y

  lemma same_sign_product:
    forall x y:t.
      class x = Finite /\ class y = Finite /\ same_sign x y
      -> value x * value y >= 0.0

352 353 354 355 356 357 358 359 360 361 362 363 364 365 366
  (** {3 overflow} *)

  (** This predicate tells what is the result of a rounding in case
      of overflow *)
  predicate overflow_value (m:mode) (x:t) =
    match m, sign x with
    | Down, Neg -> is_infinite x
    | Down, Pos -> is_finite x /\ value x = max
    | Up, Neg -> is_finite x /\ value x = - max
    | Up, Pos -> is_infinite x
    | ToZero, Neg -> is_finite x /\ value x = - max
    | ToZero, Pos -> is_finite x /\ value x = max
    | (NearestTiesToAway | NearestTiesToEven), _ -> is_infinite x
    end

MARCHE Claude's avatar
MARCHE Claude committed
367

MARCHE Claude's avatar
MARCHE Claude committed
368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387
  (** rounding in logic *)
  lemma round1 :
    forall m:mode, x:real.
      no_overflow m x ->
        is_finite (round_logic m x) /\ value(round_logic m x) = round m x

  lemma round2 :
    forall m:mode, x:real.
      not no_overflow m x ->
        same_sign_real (round_logic m x) x /\
        overflow_value m (round_logic m x)

  lemma round3 : forall m:mode, x:real. exact (round_logic m x) = x

  lemma round4 : forall m:mode, x:real. model (round_logic m x) = x

  (** rounding of zero *)

  lemma round_of_zero : forall m:mode. is_gen_zero (round_logic m 0.0)

388
  use real.Abs
MARCHE Claude's avatar
MARCHE Claude committed
389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408

  lemma round_logic_le : forall m:mode, x:real.
    is_finite (round_logic m x) -> abs (value (round_logic m x)) <= max

  lemma round_no_overflow : forall m:mode, x:real.
    abs x <= max ->
    is_finite (round_logic m x) /\ value (round_logic m x) = round m x

  constant min : real

  lemma positive_constant : forall m:mode, x:real.
    min <= x <= max ->
    is_finite (round_logic m x) /\ value (round_logic m x) > 0.0 /\
    sign (round_logic m x) = Pos

  lemma negative_constant : forall m:mode, x:real.
    - max <= x <= - min ->
    is_finite (round_logic m x) /\ value (round_logic m x) < 0.0 /\
    sign (round_logic m x) = Neg

409
  (** lemmas on `gen_zero` *)
MARCHE Claude's avatar
MARCHE Claude committed
410 411 412 413 414 415 416 417

  lemma is_gen_zero_comp1 : forall x y:t.
    is_gen_zero x /\ value x = value y /\ is_finite y -> is_gen_zero y

  lemma is_gen_zero_comp2 : forall x y:t.
    is_finite x /\ not (is_gen_zero x) /\ value x = value y
    -> not (is_gen_zero y)

MARCHE Claude's avatar
MARCHE Claude committed
418 419
end

420
module GenFloatSpecFull
421

422 423 424 425
  use real.Real
  use real.Square
  use Rounding
  use SpecialValues
426
  clone export GenFloatFull with axiom .
427 428 429 430 431 432 433 434 435 436

  (** {3 binary operations} *)

  predicate add_post (m:mode) (x:t) (y:t) (r:t) =
    (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)
    /\
437 438 439
    (is_infinite x /\ is_infinite y ->
     if same_sign x y then is_infinite r /\ same_sign r x
     else is_NaN r)
440 441 442 443
    /\
    (is_finite x /\ is_finite y /\ no_overflow m (value x + value y)
      -> is_finite r /\
        value r = round m (value x + value y) /\
444 445
        (if same_sign x y then same_sign r x
         else sign_zero_result m r))
446 447 448 449 450 451 452 453 454 455 456 457 458 459
    /\
    (is_finite x /\ is_finite y /\ not no_overflow m (value x + value y)
      -> SpecialValues.same_sign_real (sign r) (value x + value y)
         /\ overflow_value m r)
    /\ exact r = exact x + exact y
    /\ model r = model x + model y

  predicate sub_post (m:mode) (x:t) (y:t) (r:t) =
    (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)
    /\
460 461 462
    (is_infinite x /\ is_infinite y ->
     if diff_sign x y then is_infinite r /\ same_sign r x
     else is_NaN r)
463 464 465 466
    /\
    (is_finite x /\ is_finite y /\ no_overflow m (value x - value y)
      -> is_finite r /\
        value r = round m (value x - value y) /\
467 468
        (if diff_sign x y then same_sign r x
         else sign_zero_result m r))
469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491
    /\
    (is_finite x /\ is_finite y /\ not no_overflow m (value x - value y)
      -> SpecialValues.same_sign_real (sign r) (value x - value y)
         /\ overflow_value m r)
    /\ exact r = exact x - exact y
    /\ model r = model x - model y

  predicate product_sign (z x y: t) =
    (same_sign x y -> sign z = Pos) /\ (diff_sign x y -> sign z = Neg)

  predicate mul_post (m:mode) (x:t) (y:t) (r:t) =
    (is_NaN x \/ is_NaN y -> is_NaN r)
    /\ (is_gen_zero x /\ is_infinite y -> is_NaN r)
    /\ (is_finite x /\ is_infinite y /\ value x <> 0.0
         -> is_infinite r)
    /\ (is_infinite x /\ is_gen_zero y -> is_NaN r)
    /\ (is_infinite x /\ is_finite y /\ value y <> 0.0
         -> is_infinite r)
    /\ (is_infinite x /\ is_infinite y -> is_infinite r)
    /\ (is_finite x /\ is_finite y /\ no_overflow m (value x * value y)
         -> is_finite r /\ value r = round m (value x * value y))
    /\ (is_finite x /\ is_finite y /\ not no_overflow m (value x * value y)
         -> overflow_value m r)
492
    /\ (not is_NaN r -> product_sign r x y)
493 494 495 496 497 498 499
    /\ exact r = exact x * exact y
    /\ model r = model x * model y

  predicate neg_post (x:t) (r:t) =
    (is_NaN x -> is_NaN r)
    /\ (is_infinite x -> is_infinite r)
    /\ (is_finite x -> is_finite r /\ value r = - value x)
500
    /\ (not is_NaN r -> diff_sign r x)
501 502 503
    /\ exact r = - exact x
    /\ model r = - model x

504 505 506 507 508 509 510 511 512 513 514 515 516
  predicate div_post (m:mode) (x:t) (y:t) (r:t) =
    (is_NaN x \/ is_NaN y -> is_NaN r)
    /\ (is_finite x /\ is_infinite y -> is_gen_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 /\ value y <> 0.0 /\
        no_overflow m (value x / value y)
        -> is_finite r /\ value r = round m (value x / value y))
    /\ (is_finite x /\ is_finite y /\ value y <> 0.0 /\
        not no_overflow m (value x / value y)
        -> overflow_value m r)
    /\ (is_finite x /\ is_gen_zero y /\ value x <> 0.0
        -> is_infinite r)
517
    /\ (is_gen_zero x /\ is_gen_zero y -> is_NaN r)
518
    /\ (not is_NaN r -> product_sign r x y)
519 520
    /\ exact r = exact x / exact y
    /\ model r = model x / model y
521

522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547
  predicate fma_post (m:mode) (x y z:t) (r:t) =
    (is_NaN x \/ is_NaN y \/ is_NaN z -> is_NaN r)
    /\ (is_gen_zero x /\ is_infinite y -> is_NaN r)
    /\ (is_infinite x /\ is_gen_zero y -> is_NaN r)
    /\ (is_finite x /\ value x <> 0.0 /\ is_infinite y /\ is_finite z
        -> is_infinite r /\ product_sign r x y)
    /\ (is_finite x /\ value x <> 0.0 /\ 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 /\ value y <> 0.0 /\ is_finite z
        -> is_infinite r /\ product_sign r x y)
    /\ (is_infinite x /\ is_finite y /\ value y <> 0.0 /\ 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 /\
        no_overflow m (value x * value y + value z)
        -> is_finite r /\
        value r = round m (value x * value y + value z) /\
        (if product_sign z x y then same_sign r z
548 549
         else (value x * value y + value z = 0.0 ->
               if m = Down then sign r = Neg else sign r = Pos)))
550 551 552 553 554 555 556 557 558
    /\ (is_finite x /\ is_finite y /\ is_finite z /\
        not no_overflow m (value x * value y + value z)
        -> SpecialValues.same_sign_real (sign r) (value x * value y + value z)
        /\ overflow_value m r)
    /\ exact r = exact x * exact y + exact z
    /\ model r = model x * model y + model z

  predicate sqrt_post (m:mode) (x:t) (r:t) =
    (is_NaN x -> is_NaN r)
559
    /\ (is_plus_infinity x -> is_plus_infinity r)
560 561 562 563 564 565 566 567 568
    /\ (is_minus_infinity x -> is_NaN r)
    /\ (is_finite x /\ value x < 0.0 -> is_NaN r)
    /\ (is_finite x /\ value x = 0.0
        -> is_finite r /\ value r = 0.0 /\ same_sign r x)
    /\ (is_finite x /\ value x > 0.0
        -> is_finite r /\ value r = round m (sqrt (value x)) /\ sign r = Pos)
    /\ exact r = sqrt (exact x)
    /\ model r = sqrt (model x)

569 570 571 572 573
  predicate of_real_exact_post (x:real) (r:t) =
    is_finite r /\ value r = x /\ exact r = x /\ model r = x

  (** {3 Comparisons} *)

574 575 576 577 578 579 580 581 582 583 584 585 586 587 588
  predicate le (x:t) (y:t) =
    (is_finite x /\ is_finite y /\ value x <= value y)
    \/ (is_minus_infinity x /\ is_not_NaN y)
    \/ (is_not_NaN x /\ is_plus_infinity y)

  predicate lt (x:t) (y:t) =
    (is_finite x /\ is_finite y /\ value x < value 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)

  predicate ge (x:t) (y:t) = le y x

  predicate gt (x:t) (y:t) = lt y x

  predicate eq (x:t) (y:t) =
589 590
    (is_finite x /\ is_finite y /\ value x = value y) \/
    (is_infinite x /\ is_infinite y /\ same_sign x y)
591 592 593 594 595 596 597 598 599 600 601 602 603 604

  predicate ne (x:t) (y:t) = not (eq x y)

  lemma le_lt_trans:
    forall x y z:t. le x y /\ lt y z -> lt x z

  lemma lt_le_trans:
    forall x y z:t. lt x y /\ le y z -> lt x z

  lemma le_ge_asym:
    forall x y:t. le x y /\ ge x y -> eq x y

  lemma not_lt_ge: forall x y:t.
    not (lt x y) /\ is_not_NaN x /\ is_not_NaN y -> ge x y
605

606 607
  lemma not_gt_le: forall x y:t.
    not (gt x y) /\ is_not_NaN x /\ is_not_NaN y -> le x y
608 609 610 611

end


612
(** {2 Full theory of single precision floats} *)
MARCHE Claude's avatar
MARCHE Claude committed
613

614
module SingleFull
MARCHE Claude's avatar
MARCHE Claude committed
615

616
  use export SingleFormat
MARCHE Claude's avatar
MARCHE Claude committed
617

MARCHE Claude's avatar
MARCHE Claude committed
618 619
  constant min_single : real = 0x1p-149

620
  clone export GenFloatSpecFull with
Andrei Paskevich's avatar
Andrei Paskevich committed
621
    type t = single,
622
    constant min = min_single,
623
    constant max = max_single,
624
    constant max_representable_integer = max_int,
625 626
    lemma Bounded_real_no_overflow,
    axiom . (* TODO: "lemma"? "goal"? *)
MARCHE Claude's avatar
MARCHE Claude committed
627 628 629

end

630
(** {2 Full theory of double precision floats} *)
MARCHE Claude's avatar
MARCHE Claude committed
631

632
module DoubleFull
MARCHE Claude's avatar
MARCHE Claude committed
633

634
  use export DoubleFormat
MARCHE Claude's avatar
MARCHE Claude committed
635

MARCHE Claude's avatar
MARCHE Claude committed
636 637
  constant min_double : real = 0x1p-1074

638
  clone export GenFloatSpecFull with
Andrei Paskevich's avatar
Andrei Paskevich committed
639
    type t = double,
640
    constant min = min_double,
641
    constant max = max_double,
642
    constant max_representable_integer = max_int,
643 644
    lemma Bounded_real_no_overflow,
    axiom . (* TODO: "lemma"? "goal"? *)
MARCHE Claude's avatar
MARCHE Claude committed
645 646

end
647 648 649



650
module GenFloatSpecMultiRounding
651

652 653 654
  use Rounding
  use real.Real
  use real.Abs
655
  clone export GenFloat with axiom .
656 657 658 659 660 661 662 663 664

  (** {3 binary operations} *)

  constant min_normalized : real
  constant eps_normalized : real
  constant eta_normalized : real

  predicate add_post (m:mode) (x y res:t) =
    ((* CASE 1 : normalized numbers *)
665
     (abs(value x + value y) >= min_normalized ->
666
      abs(value res - (value x + value y)) <= eps_normalized * abs(value x + value y))
667
     /\
668
     (*CASE 2: denormalized numbers *)
669
     (abs(value x + value y) <= min_normalized ->
670 671 672 673 674 675 676 677 678
      abs(value res - (value x + value y)) <= eta_normalized))
   /\
    exact res = exact x + exact y
   /\
    model res = model x + model y


  predicate sub_post (m:mode) (x y res:t) =
    ((* CASE 1 : normalized numbers *)
679
     (abs(value x - value y) >= min_normalized ->
680
      abs(value res - (value x - value y)) <= eps_normalized * abs(value x - value y))
681
     /\
682
     (*CASE 2: denormalized numbers *)
683
     (abs(value x - value y) <= min_normalized ->
684 685 686 687 688 689 690 691 692
      abs(value res - (value x - value y)) <= eta_normalized))
   /\
    exact res = exact x - exact y
   /\
    model res = model x - model y


  predicate mul_post (m:mode) (x y res:t) =
    ((* CASE 1 : normalized numbers *)
693 694 695
     (abs(value x * value y) >= min_normalized ->
      abs(value res - (value x * value y)) <= eps_normalized * abs(value x * value y))
     /\
696
     (*CASE 2: denormalized numbers *)
697 698
     (abs(value x * value y) <= min_normalized ->
      abs(value res - (value x * value y)) <= eta_normalized))
699 700 701 702 703
   /\ exact res = exact x * exact y
   /\ model res = model x * model y


  predicate neg_post (m:mode) (x res:t) =
704 705 706 707 708
   ((abs(value x) >= min_normalized ->
     abs(value res - (- value x)) <= eps_normalized * abs(- value x))
    /\
    (abs(value x) <= min_normalized ->
     abs(value res - (- value x)) <= eta_normalized))
709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726
   /\ exact res = - exact x
   /\ model res = - model x

  predicate of_real_post (m:mode) (x:real) (res:t) =
    value res = round m x /\ exact res = x /\ model res = x

  predicate of_real_exact_post (x:real) (r:t) =
    value r = x /\ exact r = x /\ model r = x

  (** {3 Comparisons} *)

  predicate lt (x:t) (y:t) = value x < value y

  predicate gt (x:t) (y:t) = value x > value y

end


727
(*** TODO: find constants for type single *)
728

729
(***
730
module SingleMultiRounding
731

732 733 734
  constant min_normalized_single : real =
  constant eps_normalized_single : real =
  constant eta_normalized_single : real =
735 736 737 738 739 740 741

  clone export GenFloatSpecMultiRounding with
    type t = single,
    constant max = max_single,
    constant max_representable_integer = max_int,
    constant min_normalized = min_normalized_single,
    constant eps_normalized = eps_normalized_single,
742
    constant eta_normalized = eta_normalized_single,
743 744 745
    lemma Bounded_real_no_overflow,
    axiom . (* TODO: "lemma"? "goal"? *)

746 747 748 749 750

end

*)

751
module DoubleMultiRounding
752 753 754 755 756 757 758 759 760 761 762 763 764

  use export DoubleFormat

  constant min_normalized_double : real = 0x1p-1022
  constant eps_normalized_double : real = 0x1.004p-53
  constant eta_normalized_double : real = 0x1.002p-1075

  clone export GenFloatSpecMultiRounding with
    type t = double,
    constant max = max_double,
    constant max_representable_integer = max_int,
    constant min_normalized = min_normalized_double,
    constant eps_normalized = eps_normalized_double,
765
    constant eta_normalized = eta_normalized_double,
766 767
    lemma Bounded_real_no_overflow,
    axiom . (* TODO: "lemma"? "goal"? *)
768

MARCHE Claude's avatar
MARCHE Claude committed
769
end