ieeeScript.sml

1
2(* ========================================================================= *)
3(* Formalization of IEEE-754 Standard for binary floating-point arithmetic.  *)
4(* ========================================================================= *)
5
6(*---------------------------------------------------------------------------*
7 * First, make standard environment available.                               *
8 *---------------------------------------------------------------------------*)
9Theory ieee
10Ancestors
11  real pred_set arithmetic
12Libs
13  Num_conv
14
15(*---------------------------------------------------------------------------*
16 * Next, bring in extra tools used.                                          *
17 *---------------------------------------------------------------------------*)
18
19(*---------------------------------------------------------------------------*
20 * Create the theory.                                                        *
21 *---------------------------------------------------------------------------*)
22
23val _ = ParseExtras.temp_loose_equality()
24
25(* ------------------------------------------------------------------------- *)
26(* Derived parameters for floating point formats.                            *)
27(* ------------------------------------------------------------------------- *)
28
29Definition expwidth[nocompute]:
30  expwidth (ew:num,fw:num) = ew
31End
32
33Definition fracwidth[nocompute]:
34  fracwidth (ew:num,fw:num) = fw
35End
36
37Definition wordlength[nocompute]:
38  wordlength(X: (num#num)) = (expwidth(X) + fracwidth(X) + (1:num))
39End
40
41Definition emax[nocompute]:
42  emax(X: (num#num)) = (((2:num) EXP (expwidth (X))) - (1:num))
43End
44
45Definition bias[nocompute]:
46  bias(X: (num#num)) = ((2:num) EXP (expwidth(X) - (1:num)) - (1:num))
47End
48
49(* ------------------------------------------------------------------------- *)
50(* Predicates for the four IEEE formats.                                     *)
51(* ------------------------------------------------------------------------- *)
52
53Definition is_single[nocompute]:
54  is_single (X: (num#num)) = (expwidth(X) = (8:num)) /\ (wordlength(X) = (32:num))
55End
56
57Definition is_double[nocompute]:
58  is_double(X: (num#num)) = (expwidth(X) = (11:num)) /\ (wordlength(X) = (64:num))
59End
60
61Definition is_single_extended[nocompute]:
62  is_single_extended(X) = expwidth(X) >= (11:num) /\ wordlength(X) >= (43:num)
63End
64
65Definition is_double_extended[nocompute]:
66  is_double_extended(X) = expwidth(X) >= (15:num) /\ wordlength(X) >= (79:num)
67End
68
69(* ------------------------------------------------------------------------- *)
70(* Extractors for fields.                                                    *)
71(* ------------------------------------------------------------------------- *)
72
73Definition sign[nocompute]:
74  sign ((s:num),(e:num),(f:num)) = (s:num)
75End
76
77Definition exponent[nocompute]:
78  exponent ((s:num),(e:num),(f:num)) = (e:num)
79End
80
81Definition fraction[nocompute]:
82  fraction ((s:num),(e:num),(f:num)) = (f:num)
83End
84
85(* ------------------------------------------------------------------------- *)
86(* Partition of numbers into disjoint classes.                               *)
87(* ------------------------------------------------------------------------- *)
88
89Definition is_nan[nocompute]:
90  is_nan(X: (num#num)) (a:num#num#num) =
91  (exponent (a:num#num#num) = emax(X)) /\ ~(fraction (a:num#num#num) = (0:num))
92End
93
94Definition is_infinity[nocompute]:
95  is_infinity((X: num#num)) (a:num#num#num) = (exponent a = emax(X)) /\ (fraction a = (0:num))
96End
97
98Definition is_normal[nocompute]:
99  is_normal(X: (num#num)) (a:num#num#num) = (((0:num) < exponent a) /\ (exponent a < emax(X)))
100End
101
102Definition is_denormal[nocompute]:
103  is_denormal(X:num#num) (a:num#num#num) = ((exponent a = (0:num)) /\ ~(fraction a = (0:num)))
104End
105
106Definition is_zero[nocompute]:
107  is_zero (X:num#num) (a:num#num#num) = ((exponent a = 0) /\ (fraction a = 0))
108End
109
110(* ------------------------------------------------------------------------- *)
111(* Other useful predicates.                                                  *)
112(* ------------------------------------------------------------------------- *)
113
114Definition is_valid[nocompute]:
115  is_valid(X:num#num) (s:num,e:num,f:num) =
116  ((s < SUC(SUC 0)) /\ (e < 2 EXP expwidth(X)) /\ (f < 2 EXP fracwidth(X)))
117End
118
119Definition is_finite[nocompute]:
120  is_finite(X:num#num) (a : num#num#num) =
121  ((is_valid (X) a) /\ ((is_normal(X:num#num) a) \/ (is_denormal(X:num#num) a) \/ (is_zero (X:num#num) a)))
122End
123
124(* ------------------------------------------------------------------------- *)
125(* Some special values.                                                      *)
126(* ------------------------------------------------------------------------- *)
127
128Definition plus_infinity[nocompute]:
129  plus_infinity(X:num#num) = (0:num,emax(X),0:num)
130End
131
132Definition minus_infinity[nocompute]:
133  minus_infinity(X:num#num) = (1:num,emax(X),0:num)
134End
135
136Definition plus_zero[nocompute]:
137  plus_zero(X:num#num) = (0:num,0:num,0:num)
138End
139
140Definition minus_zero[nocompute]:
141  minus_zero(X:num#num) = (1:num,0:num,0:num)
142End
143
144Definition topfloat[nocompute]:
145  topfloat(X) = (0:num, (emax (X:num#num) - (1:num)) , (2 EXP fracwidth(X) - (1:num)))
146End
147
148Definition bottomfloat[nocompute]:
149  bottomfloat (X:num#num) = ((1:num), (emax(X) - 1) , (2 EXP fracwidth(X) - 1))
150End
151
152(* ------------------------------------------------------------------------- *)
153(* Negation operation on floating point values.                              *)
154(* ------------------------------------------------------------------------- *)
155
156Definition minus[nocompute]:
157  minus(X:num#num) (a : num#num#num)= ((1 - sign(a)),exponent(a),fraction(a))
158End
159
160(* ------------------------------------------------------------------------- *)
161(* Concrete encodings (at least valid for single and double).                *)
162(* ------------------------------------------------------------------------- *)
163
164Definition encoding[nocompute]:
165  encoding(X:num#num) (s:num,e:num,f:num) =
166  ((s * 2 EXP (wordlength(X) - 1)) + (e * 2 EXP fracwidth(X)) + f)
167End
168
169(* ------------------------------------------------------------------------- *)
170(* Real number valuations.                                                   *)
171(* ------------------------------------------------------------------------- *)
172
173Definition valof[nocompute]:
174  valof (X:num#num) (s:num,e:num,f:num) =
175  (if (e = 0) then  ~(&1) pow s * (&2 / &2 pow bias(X)) * (&f / &2 pow (fracwidth(X)))
176    else  ~(&1) pow s * (&2 pow e / &2 pow bias(X)) * (&1 + &f / &2 pow fracwidth(X)))
177End
178
179(* ------------------------------------------------------------------------- *)
180(* A few handy values.                                                       *)
181(* ------------------------------------------------------------------------- *)
182
183Definition largest[nocompute]:
184  largest(X:num#num) = (&2 pow (emax(X) - 1) / &2 pow bias(X)) * (&2 - inv(&2 pow fracwidth(X)))
185End
186
187Definition threshold[nocompute]:
188  threshold(X:num#num) = (&2 pow (emax(X) - 1) / &2 pow bias(X)) * (&2 - inv(&2 pow SUC(fracwidth(X))))
189End
190
191Definition ulp[nocompute]:
192  ulp(X:num#num) (a :num#num#num) = valof(X) (0,exponent(a),1) - valof(X) (0,exponent(a),0)
193End
194
195(* ------------------------------------------------------------------------- *)
196(* Enumerated type for rounding modes.                                       *)
197(* ------------------------------------------------------------------------- *)
198
199Datatype:
200   roundmode = To_nearest
201  | float_To_zero
202  | To_pinfinity
203  | To_ninfinity
204End
205
206(* ------------------------------------------------------------------------- *)
207(* Characterization of best approximation from a set of abstract values.     *)
208(* ------------------------------------------------------------------------- *)
209
210Definition is_closest[nocompute]:
211  is_closest (v) (s) (x) (a) = (a IN s) /\ (!b. (b IN s) ==> abs(v(a) - x) <= abs(v(b) - x))
212End
213
214(* ------------------------------------------------------------------------- *)
215(* Best approximation with a deciding preference for multiple possibilities. *)
216(* ------------------------------------------------------------------------- *)
217
218Definition closest[nocompute]:
219  closest (v) (p) (s) (x) = @(a). is_closest v s x a /\
220  ((?b. is_closest v s x b /\ p(b)) ==> p(a))
221End
222
223(* ------------------------------------------------------------------------- *)
224(* Rounding to floating point formats.                                       *)
225(* ------------------------------------------------------------------------- *)
226
227Definition round_def:   (round(X:num#num) (To_nearest) (x:real) =
228  (if (x <= ~(threshold(X))) then (minus_infinity(X))
229   else if (x >= threshold(X)) then (plus_infinity(X))
230   else (closest (valof(X)) (\a. EVEN(fraction(a)))
231  { a | is_finite(X) a } x))) /\
232
233  (round(X) (float_To_zero) x =
234  (if (x < ~(largest(X))) then (bottomfloat(X))
235   else if (x > largest(X)) then (topfloat(X))
236   else (closest (valof(X)) (\x. T)
237  { a | is_finite(X) a /\ abs(valof(X) a) <= abs(x) } x))) /\
238
239  (round(X) (To_pinfinity) x =
240   if x < ~(largest(X)) then bottomfloat(X)
241   else if (x > largest(X)) then plus_infinity(X)
242   else closest (valof(X)) (\x. T)
243  { a | is_finite(X) a /\ valof(X) a >= x } x) /\
244
245  (round(X) (To_ninfinity) x =
246   if x < ~(largest(X)) then minus_infinity(X)
247   else if (x > largest(X)) then topfloat(X)
248   else closest (valof(X)) (\x. T)
249  { a | is_finite(X) a /\ valof(X) a <= x } x)
250End
251
252(* ------------------------------------------------------------------------- *)
253(* Rounding to integer values in floating point format.                      *)
254(* ------------------------------------------------------------------------- *)
255
256Definition is_integral[nocompute]:
257  is_integral(X:num#num) (a:(num#num#num)) = is_finite(X) a /\ ?n. abs(valof(X) a) = &n
258End
259
260Definition intround_def:   (intround(X:num#num) (To_nearest) (x:real) =
261  (if (x <= ~(threshold(X))) then (minus_infinity(X))
262   else if (x >= threshold(X)) then (plus_infinity(X))
263   else (closest (valof(X)) (\a. (?n. (EVEN n) /\ (abs(valof(X) a) = &n)))
264  { a | is_integral(X) a} x))) /\
265
266  (intround(X) float_To_zero x =
267  (if (x < ~(largest(X))) then (bottomfloat(X))
268   else if (x > largest(X)) then (topfloat(X))
269   else (closest (valof(X)) (\x. T)
270  { a | is_integral(X) a /\ abs(valof(X) a) <= abs(x) } x))) /\
271
272  (intround(X) To_pinfinity x =
273  (if (x < ~(largest(X))) then (bottomfloat(X))
274   else if (x > largest(X)) then (plus_infinity(X))
275   else (closest (valof(X)) (\x. T)
276  { a | is_integral(X) a /\ valof(X) a >= x } x))) /\
277
278  (intround(X) To_ninfinity x =
279   if (x < ~(largest(X))) then (minus_infinity(X))
280   else if (x > largest(X)) then (topfloat(X))
281   else (closest (valof(X)) (\x. T)
282  { a | is_integral(X) a /\ valof(X) a <= x } x))
283End
284
285(* ------------------------------------------------------------------------- *)
286(* A hack for our (non-standard) treatment of NaNs.                          *)
287(* ------------------------------------------------------------------------- *)
288
289Definition some_nan[nocompute]:
290  some_nan(X:num#num) = @(a:num#num#num). is_nan(X) a
291End
292
293(* ------------------------------------------------------------------------- *)
294(* Coercion for signs of zero results.                                       *)
295(* ------------------------------------------------------------------------- *)
296
297Definition zerosign[nocompute]:
298  zerosign (X:num#num) (s:num) (a:num#num#num) = (if (is_zero(X) a) then
299  (if (s = 0) then plus_zero(X) else minus_zero(X)) else a)
300End
301
302(* ------------------------------------------------------------------------- *)
303(* Useful mathematical operations not already in the HOL Light core.         *)
304(* ------------------------------------------------------------------------- *)
305
306val rem = new_infixl_definition (
307  "rem", Term`$rem x y = let n = closest I (\x. ?n. EVEN(n) /\ (abs(x) = &n))
308  { x | ?n. abs(x) = &n } (x / y) in x - n * y`, 500);
309
310(* ------------------------------------------------------------------------- *)
311(* Definitions of the arithmetic operations.                                 *)
312(* ------------------------------------------------------------------------- *)
313
314Definition fintrnd[nocompute]:
315  fintrnd(X:num#num) (m:roundmode) (a:num#num#num) =
316  if is_nan(X) a then some_nan(X)
317    else if is_infinity(X) a then a
318      else zerosign(X) (sign(a)) (intround(X) m (valof(X) a))
319End
320
321Definition fadd[nocompute]:
322  fadd(X:num#num) (m:roundmode) (a:num#num#num) (b:num#num#num) =
323  if (is_nan(X) a) \/ (is_nan(X) b) \/ ((is_infinity(X) a) /\ (is_infinity(X) b) /\ (~(sign(a) = sign(b)))) then (some_nan(X))
324  else if is_infinity(X) a then a
325  else if is_infinity(X) b then b
326  else zerosign(X) (if is_zero(X) a /\ is_zero(X) b /\ (sign(a) = sign(b)) then sign(a) else if (m = To_ninfinity) then 1 else 0) (round(X) m (valof(X) a + valof(X) b))
327End
328
329Definition fsub[nocompute]:
330  fsub(X:num#num) (m:roundmode) (a:num#num#num) (b:num#num#num) =
331  (if is_nan(X) a \/ is_nan(X) b \/ (is_infinity(X) a /\ is_infinity(X) b /\ (sign(a) = sign(b))) then some_nan(X)
332   else if is_infinity(X) a then a
333   else if is_infinity(X) b then minus(X) b
334   else zerosign(X) (if is_zero(X) a /\ is_zero(X) b /\ ~(sign(a) = sign(b)) then sign(a) else if m = To_ninfinity then 1 else 0) (round(X) m (valof(X) a - valof(X) b)))
335End
336
337Definition fmul[nocompute]:
338  fmul(X:num#num) (m:roundmode) (a:num#num#num) (b:num#num#num) =
339  (if is_nan(X) a \/ is_nan(X) b \/ is_zero(X) a /\ is_infinity(X) b \/ is_infinity(X) a /\ is_zero(X) b then some_nan(X)
340   else if is_infinity(X) a \/ is_infinity(X) b then (if sign(a) = sign(b) then plus_infinity(X) else minus_infinity(X))
341   else zerosign(X) (if sign(a) = sign(b) then 0 else 1) (round(X) m (valof(X) a * valof(X) b)))
342End
343
344Definition fdiv[nocompute]:
345  fdiv(X:num#num) (m:roundmode) (a:num#num#num) (b:num#num#num) =
346  (if is_nan(X) a \/ is_nan(X) b \/ is_zero(X) a /\ is_zero(X) b \/ is_infinity(X) a /\ is_infinity(X) b then some_nan(X)
347   else if is_infinity(X) a \/ is_zero(X) b then (if sign(a) = sign(b) then plus_infinity(X) else minus_infinity(X))
348   else if is_infinity(X) b then (if sign(a) = sign(b) then plus_zero(X) else minus_zero(X))
349   else zerosign(X) (if sign(a) = sign(b) then 0 else 1) (round(X) m (valof(X) a / valof(X) b)))
350End
351
352Definition fsqrt[nocompute]:
353  fsqrt (X:num#num) (m:roundmode) (a:num#num#num) =
354  (if is_nan(X) a then some_nan(X)
355   else if is_zero(X) a \/ is_infinity(X) a /\ (sign(a) = 0) then a
356   else if (sign(a) = 1) then some_nan(X)
357   else zerosign(X) (sign(a)) (round(X) m (sqrt(valof(X) a))))
358End
359
360
361Definition frem[nocompute]:
362  frem(X:num#num) (m:roundmode) (a:num#num#num) (b:num#num#num) =
363  (if is_nan(X) a \/ is_nan(X) b \/ is_infinity(X) a \/ is_zero(X) b then some_nan(X)
364   else if is_infinity(X) b then a
365   else zerosign(X) (sign(a)) (round(X) m (valof(X) a rem valof(X) b)))
366End
367
368(* ------------------------------------------------------------------------- *)
369(* Negation is specially simple.                                             *)
370(* ------------------------------------------------------------------------- *)
371
372Definition fneg[nocompute]:
373  fneg(X:num#num) (m:roundmode) (a:num#num#num) = (((1:num)-sign(a)),(exponent(a)),(fraction(a)))
374End
375
376(* ------------------------------------------------------------------------- *)
377(* Comparison codes.                                                         *)
378(* ------------------------------------------------------------------------- *)
379
380Datatype:
381  ccode = Gt | Lt | Eq | Un
382End
383
384(* ------------------------------------------------------------------------- *)
385(* Comparison operations.                                                    *)
386(* ------------------------------------------------------------------------- *)
387
388Definition fcompare[nocompute]:
389  fcompare(X) a b =
390  (if is_nan(X) a \/ is_nan(X) b then Un
391   else if is_infinity(X) a /\ (sign(a) = 1) then (if is_infinity(X) b /\ (sign(b) = 1) then Eq else Lt)
392   else if is_infinity(X) a /\ (sign(a) = 0) then (if is_infinity(X) b /\ (sign(b) = 0) then Eq else Gt)
393   else if is_infinity(X) b /\ (sign(b) = 1) then Gt
394   else if is_infinity(X) b /\ (sign(b) = 0) then Lt
395   else if valof(X) a < valof(X) b then Lt
396   else if valof(X) a = valof(X) b then Eq
397   else Gt)
398End
399
400Definition flt[nocompute]:
401  flt(X) a b = (fcompare(X) a b = Lt)
402End
403
404Definition fle[nocompute]:
405  fle(X) a b = (fcompare(X) a b = Lt) \/ (fcompare(X) a b = Eq)
406End
407
408Definition fgt[nocompute]:
409  fgt(X) a b = (fcompare(X) a b = Gt)
410End
411
412Definition fge[nocompute]:
413  fge(X) a b = (fcompare(X) a b = Gt) \/ (fcompare(X) a b = Eq)
414End
415
416Definition feq[nocompute]:
417  feq (X) a b = (fcompare(X) a b = Eq)
418End
419
420(* ------------------------------------------------------------------------- *)
421(* Actual float type with round-to-even.                                     *)
422(* ------------------------------------------------------------------------- *)
423
424Definition float_format[nocompute]:
425  float_format = ((8:num),(23:num))
426End
427  (* for double
428  “float_format = ((11:num),(52:num))”);
429  *)
430
431val float_tybij = define_new_type_bijections {
432  name="float_tybij",
433  ABS="float",
434  REP="defloat",
435  tyax =new_type_definition ("float",
436  Q.prove (`(?a. (is_valid float_format a))`,
437  EXISTS_TAC (“0:num,0:num,0:num” ) THEN
438  REWRITE_TAC[float_format, is_valid, GSYM NOT_LESS_EQUAL,
439  LE, num_CONV“2:num”, NOT_EXP_0, GSYM SUC_NOT]))};
440
441Definition Val[nocompute]:
442  Val a = valof(float_format) (defloat a)
443End
444
445Definition Float[nocompute]:
446  Float(x) = float (round(float_format) To_nearest x)
447End
448
449Definition Sign[nocompute]:
450  Sign(a) = sign(defloat a)
451End
452
453Definition Exponent[nocompute]:
454  Exponent(a) = exponent(defloat a)
455End
456
457Definition Fraction[nocompute]:
458  Fraction(a) = fraction(defloat a)
459End
460
461Definition Ulp[nocompute]:
462  Ulp(a) = ulp(float_format) (defloat a)
463End
464
465(* ------------------------------------------------------------------------- *)
466(* Lifting of the discriminator functions.                                   *)
467(* ------------------------------------------------------------------------- *)
468
469Definition Isnan[nocompute]:
470  Isnan(a) = is_nan(float_format) (defloat a)
471End
472
473Definition Infinity[nocompute]:
474  Infinity(a) = is_infinity(float_format) (defloat a)
475End
476
477Definition Isnormal[nocompute]:
478  Isnormal(a) = is_normal(float_format) (defloat a)
479End
480
481Definition Isdenormal[nocompute]:
482  Isdenormal(a) = is_denormal(float_format) (defloat a)
483End
484
485Definition Iszero[nocompute]:
486  Iszero(a) = is_zero(float_format) (defloat a)
487End
488
489Definition Finite[nocompute]:
490  Finite(a) = Isnormal(a) \/ Isdenormal(a) \/ Iszero(a)
491End
492
493Definition Isintegral[nocompute]:
494  Isintegral(a) = is_integral(float_format) (defloat a)
495End
496
497(* ------------------------------------------------------------------------- *)
498(* Basic operations on floats.                                               *)
499(* ------------------------------------------------------------------------- *)
500
501Definition Plus_zero[nocompute]:
502  Plus_zero = float (plus_zero(float_format))
503End
504
505Definition Minus_zero[nocompute]:
506  Minus_zero = float (minus_zero(float_format))
507End
508
509Definition Minus_infinity[nocompute]:
510  Minus_infinity = float (minus_infinity(float_format))
511End
512
513Definition Plus_infinity[nocompute]:
514  Plus_infinity = float (plus_infinity(float_format))
515End
516
517val float_add = new_infixl_definition (
518  "float_add",
519  Term`$float_add a b = float (fadd(float_format) To_nearest (defloat a) (defloat b))`, 500);
520
521Overload "+" = Term`$float_add`
522
523val float_sub =
524    new_infixl_definition
525    ("float_sub",
526     Term`$float_sub a b = float (fsub(float_format) To_nearest (defloat a) (defloat b))`,
527     500);
528
529Overload "-" = Term`$float_sub`
530
531val float_mul = new_infixl_definition (
532  "float_mul",
533  Term`$float_mul a b = float (fmul(float_format) To_nearest (defloat a) (defloat b))`, 500);
534
535Overload "*" = Term`$float_mul`
536
537val float_div = new_infixl_definition (
538  "float_div",
539  Term`$float_div a b = float (fdiv(float_format) To_nearest (defloat a) (defloat b))`, 500);
540
541Overload "/" = Term`$float_div`
542
543val float_rem = new_infixl_definition (
544  "float_rem",
545  Term`$float_rem a b = float (frem(float_format) To_nearest (defloat a) (defloat b))`, 500);
546
547Definition float_sqrt[nocompute]:
548  float_sqrt(a) = float (fsqrt(float_format) To_nearest (defloat a))
549End
550
551Definition ROUNDFLOAT[nocompute]:
552  ROUNDFLOAT(a) = float (fintrnd(float_format) To_nearest (defloat a))
553End
554
555Definition float_lt[nocompute]:
556  float_lt a b = flt(float_format) (defloat a) (defloat b)
557End
558Overload "<" = “$float_lt”
559
560Definition float_le[nocompute]:
561  float_le a b = fle(float_format) (defloat a) (defloat b)
562End
563Overload "<=" = “$float_le”
564
565Definition float_gt[nocompute]:
566  float_gt a b = fgt(float_format) (defloat a) (defloat b)
567End
568
569Overload ">" = “$float_gt”
570
571Definition float_ge[nocompute]:
572  float_ge a b = fge(float_format) (defloat a) (defloat b)
573End
574Overload ">=" = “$float_ge”
575
576
577Definition float_eq[nocompute]:
578  float_eq (a:float) (b:float) = feq(float_format) (defloat a) (defloat b)
579End
580Overload "==" = Term`$float_eq`
581val _ = set_fixity "==" (Infix(NONASSOC,450))
582
583Definition float_neg[nocompute]:
584  float_neg (a:float) = float (fneg (float_format) To_nearest (defloat (a:float)))
585End
586
587Overload "~" = Term`$float_neg`
588
589Definition float_abs[nocompute]:
590  float_abs a = (if a >= Plus_zero then a else (float_neg a))
591End
592
593(*---------------------------------------------------------------------------*
594 * Write the theory to disk.                                                 *
595 *---------------------------------------------------------------------------*)