?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
drd/tests/swapcontext: Improve the portability of this test further
[valgrind.git] / none / tests / ppc32 / round.c
blob5a4b91c56fe9fb9431fc7f271519daf4bab192ee
1
2 /* Copyright (C) 2006 Dave Nomura
3 dcnltc@us.ibm.com
4
5 This program is free software; you can redistribute it and/or
6 modify it under the terms of the GNU General Public License as
7 published by the Free Software Foundation; either version 2 of the
8 License, or (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful, but
11 WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, see <http://www.gnu.org/licenses/>.
17
18 The GNU General Public License is contained in the file COPYING.
19 */
20
21 #include <stdio.h>
22 #include <stdlib.h>
23 #include <limits.h>
24
25 typedef enum { FALSE=0, TRUE } bool_t;
26
27 typedef enum {
28 FADDS, FSUBS, FMULS, FDIVS,
29 FMADDS, FMSUBS, FNMADDS, FNMSUBS,
30 FADD, FSUB, FMUL, FDIV, FMADD,
31 FMSUB, FNMADD, FNMSUB, FSQRT
32 } flt_op_t;
33
34 typedef enum {
35 TO_NEAREST=0, TO_ZERO, TO_PLUS_INFINITY, TO_MINUS_INFINITY } round_mode_t;
36 char *round_mode_name[] = { "near", "zero", "+inf", "-inf" };
37
38 const char *flt_op_names[] = {
39 "fadds", "fsubs", "fmuls", "fdivs",
40 "fmadds", "fmsubs", "fnmadds", "fnmsubs",
41 "fadd", "fsub", "fmul", "fdiv", "fmadd", "fmsub", "fnmadd",
42 "fnmsub", "fsqrt"
43 };
44
45 typedef unsigned int fpscr_t;
46
47 typedef union {
48 float flt;
49 struct {
50 #if defined(VGP_ppc64le_linux)
51 unsigned int frac:23;
52 unsigned int exp:8;
53 unsigned int sign:1;
54 #else
55 unsigned int sign:1;
56 unsigned int exp:8;
57 unsigned int frac:23;
58 #endif
59 } layout;
60 } flt_overlay;
61
62 typedef union {
63 double dbl;
64 struct {
65 #if defined(VGP_ppc64le_linux)
66 unsigned int frac_lo:32;
67 unsigned int frac_hi:20;
68 unsigned int exp:11;
69 unsigned int sign:1;
70 #else
71 unsigned int sign:1;
72 unsigned int exp:11;
73 unsigned int frac_hi:20;
74 unsigned int frac_lo:32;
75 #endif
76 } layout;
77 struct {
78 unsigned int hi;
79 unsigned int lo;
80 } dbl_pair;
81 } dbl_overlay;
82
83 void assert_fail(const char *msg,
84 const char* expr, const char* file, int line, const char*fn);
85
86 #define STRING(__str) #__str
87 #define assert(msg, expr) \
88 ((void) ((expr) ? 0 : \
89 (assert_fail (msg, STRING(expr), \
90 __FILE__, __LINE__, \
91 __PRETTY_FUNCTION__), 0)))
92 float denorm_small;
93 double dbl_denorm_small;
94 float norm_small;
95 bool_t debug = FALSE;
96 bool_t long_is_64_bits = sizeof(long) == 8;
97
98 void assert_fail (msg, expr, file, line, fn)
99 const char* msg;
100 const char* expr;
101 const char* file;
102 int line;
103 const char*fn;
104 {
105 printf( "\n%s: %s:%d (%s): Assertion `%s' failed.\n",
106 msg, file, line, fn, expr );
107 exit( 1 );
108 }
109 void set_rounding_mode(round_mode_t mode)
110 {
111 switch(mode) {
112 case TO_NEAREST:
113 asm volatile("mtfsfi 7, 0");
114 break;
115 case TO_ZERO:
116 asm volatile("mtfsfi 7, 1");
117 break;
118 case TO_PLUS_INFINITY:
119 asm volatile("mtfsfi 7, 2");
120 break;
121 case TO_MINUS_INFINITY:
122 asm volatile("mtfsfi 7, 3");
123 break;
124 }
125 }
126
127 void print_double(char *msg, double dbl)
128 {
129 dbl_overlay D;
130 D.dbl = dbl;
131
132 printf("%15s : dbl %-20a = %c(%4d, %05x%08x)\n",
133 msg, D.dbl, (D.layout.sign == 0 ? '+' : '-'),
134 D.layout.exp, D.layout.frac_hi, D.layout.frac_lo);
135 }
136
137 void print_single(char *msg, float *flt)
138 {
139 flt_overlay F;
140 F.flt = *flt;
141
142 /* NOTE: for the purposes of comparing the fraction of a single with
143 ** a double left shift the .frac so that hex digits are grouped
144 ** from left to right. this is necessary because the size of a
145 ** single mantissa (23) bits is not a multiple of 4
146 */
147 printf("%15s : flt %-20a = %c(%4d, %06x)\n",
148 msg, F.flt, (F.layout.sign == 0 ? '+' : '-'), F.layout.exp, F.layout.frac << 1);
149 }
150
151 int check_dbl_to_flt_round(round_mode_t mode, double dbl, float *expected)
152 {
153 int status = 0;
154 flt_overlay R, E;
155 char *result;
156
157 set_rounding_mode(mode);
158
159 E.flt = *expected;
160 R.flt = (float)dbl;
161
162 if ((R.layout.sign != E.layout.sign) ||
163 (R.layout.exp != E.layout.exp) ||
164 (R.layout.frac != E.layout.frac)) {
165 result = "FAILED";
166 status = 1;
167 } else {
168 result = "PASSED";
169 status = 0;
170 }
171 printf("%s:%s:(double)(%-20a) = %20a",
172 round_mode_name[mode], result, R.flt, dbl);
173 if (status) {
174 print_single("\n\texpected", &E.flt);
175 print_single("\n\trounded ", &R.flt);
176 }
177 putchar('\n');
178 return status;
179 }
180
181 int test_dbl_to_float_convert(char *msg, float *base)
182 {
183 int status = 0;
184 double half = (double)denorm_small/2;
185 double qtr = half/2;
186 double D_hi = (double)*base + half + qtr;
187 double D_lo = (double)*base + half - qtr;
188 float F_lo = *base;
189 float F_hi = F_lo + denorm_small;
190
191
192 /*
193 ** .....+-----+-----+-----+-----+---....
194 ** ^F_lo ^ ^ ^
195 ** D_lo
196 ** D_hi
197 ** F_hi
198 ** F_lo and F_hi are two consecutive single float model numbers
199 ** denorm_small distance apart. D_lo and D_hi are two numbers
200 ** within that range that are not representable as single floats
201 ** and will be rounded to either F_lo or F_hi.
202 */
203 printf("-------------------------- %s --------------------------\n", msg);
204 if (debug) {
205 print_double("D_lo", D_lo);
206 print_double("D_hi", D_hi);
207 print_single("F_lo", &F_lo);
208 print_single("F_hi", &F_hi);
209 }
210
211 /* round to nearest */
212 status |= check_dbl_to_flt_round(TO_NEAREST, D_hi, &F_hi);
213 status |= check_dbl_to_flt_round(TO_NEAREST, D_lo, &F_lo);
214
215 /* round to zero */
216 status |= check_dbl_to_flt_round(TO_ZERO, D_hi, (D_hi > 0 ? &F_lo : &F_hi));
217 status |= check_dbl_to_flt_round(TO_ZERO, D_lo, (D_hi > 0 ? &F_lo : &F_hi));
218
219 /* round to +inf */
220 status |= check_dbl_to_flt_round(TO_PLUS_INFINITY, D_hi, &F_hi);
221 status |= check_dbl_to_flt_round(TO_PLUS_INFINITY, D_lo, &F_hi);
222
223 /* round to -inf */
224 status |= check_dbl_to_flt_round(TO_MINUS_INFINITY, D_hi, &F_lo);
225 status |= check_dbl_to_flt_round(TO_MINUS_INFINITY, D_lo, &F_lo);
226 return status;
227 }
228
229 void
230 init()
231 {
232 flt_overlay F;
233 dbl_overlay D;
234
235 /* small is the smallest denormalized single float number */
236 F.layout.sign = 0;
237 F.layout.exp = 0;
238 F.layout.frac = 1;
239 denorm_small = F.flt; /* == 2^(-149) */
240 if (debug) {
241 print_single("float small", &F.flt);
242 }
243
244 D.layout.sign = 0;
245 D.layout.exp = 0;
246 D.layout.frac_hi = 0;
247 D.layout.frac_lo = 1;
248 dbl_denorm_small = D.dbl; /* == 2^(-1022) */
249 if (debug) {
250 print_double("double small", D.dbl);
251 }
252
253 /* n_small is the smallest normalized single precision float */
254 F.layout.exp = 1;
255 norm_small = F.flt;
256 }
257
258 int check_int_to_flt_round(round_mode_t mode, long L, float *expected)
259 {
260 int status = 0;
261 int I = L;
262 char *int_name = "int";
263 flt_overlay R, E;
264 char *result;
265 int iter;
266
267 set_rounding_mode(mode);
268 E.flt = *expected;
269
270 for (iter = 0; iter < 2; iter++) {
271 int stat = 0;
272 R.flt = (iter == 0 ? (float)I : (float)L);
273
274 if ((R.layout.sign != E.layout.sign) ||
275 (R.layout.exp != E.layout.exp) ||
276 (R.layout.frac != E.layout.frac)) {
277 result = "FAILED";
278 stat = 1;
279 } else {
280 result = "PASSED";
281 stat = 0;
282 }
283 printf("%s:%s:(float)(%4s)%9d = %11.1f",
284 round_mode_name[mode], result, int_name, I, R.flt);
285 if (stat) {
286 print_single("\n\texpected: %.1f ", &E.flt);
287 print_single("\n\trounded ", &R.flt);
288 }
289 putchar('\n');
290 status |= stat;
291
292 if (!long_is_64_bits) break;
293 int_name = "long";
294 }
295 return status;
296 }
297
298 int check_long_to_dbl_round(round_mode_t mode, long L, double *expected)
299 {
300 int status = 0;
301 dbl_overlay R, E;
302 char *result;
303
304 set_rounding_mode(mode);
305 E.dbl = *expected;
306
307 R.dbl = (double)L;
308
309 if ((R.layout.sign != E.layout.sign) ||
310 (R.layout.exp != E.layout.exp) ||
311 (R.layout.frac_lo != E.layout.frac_lo) ||
312 (R.layout.frac_hi != E.layout.frac_hi)) {
313 result = "FAILED";
314 status = 1;
315 } else {
316 result = "PASSED";
317 status = 0;
318 }
319 printf("%s:%s:(double)(%18ld) = %20.1f",
320 round_mode_name[mode], result, L, R.dbl);
321 if (status) {
322 printf("\n\texpected %.1f : ", E.dbl);
323 }
324 putchar('\n');
325 return status;
326 }
327
328 int test_int_to_float_convert(char *msg)
329 {
330 int status = 0;
331 int int24_hi = 0x03ff0fff;
332 int int24_lo = 0x03ff0ffd;
333 float pos_flt_lo = 67047420.0;
334 float pos_flt_hi = 67047424.0;
335 float neg_flt_lo = -67047420.0;
336 float neg_flt_hi = -67047424.0;
337
338 printf("-------------------------- %s --------------------------\n", msg);
339 status |= check_int_to_flt_round(TO_NEAREST, int24_lo, &pos_flt_lo);
340 status |= check_int_to_flt_round(TO_NEAREST, int24_hi, &pos_flt_hi);
341 status |= check_int_to_flt_round(TO_ZERO, int24_lo, &pos_flt_lo);
342 status |= check_int_to_flt_round(TO_ZERO, int24_hi, &pos_flt_lo);
343 status |= check_int_to_flt_round(TO_PLUS_INFINITY, int24_lo, &pos_flt_hi);
344 status |= check_int_to_flt_round(TO_PLUS_INFINITY, int24_hi, &pos_flt_hi);
345 status |= check_int_to_flt_round(TO_MINUS_INFINITY, int24_lo, &pos_flt_lo);
346 status |= check_int_to_flt_round(TO_MINUS_INFINITY, int24_hi, &pos_flt_lo);
347
348 status |= check_int_to_flt_round(TO_NEAREST, -int24_lo, &neg_flt_lo);
349 status |= check_int_to_flt_round(TO_NEAREST, -int24_hi, &neg_flt_hi);
350 status |= check_int_to_flt_round(TO_ZERO, -int24_lo, &neg_flt_lo);
351 status |= check_int_to_flt_round(TO_ZERO, -int24_hi, &neg_flt_lo);
352 status |= check_int_to_flt_round(TO_PLUS_INFINITY, -int24_lo, &neg_flt_lo);
353 status |= check_int_to_flt_round(TO_PLUS_INFINITY, -int24_hi, &neg_flt_lo);
354 status |= check_int_to_flt_round(TO_MINUS_INFINITY, -int24_lo, &neg_flt_hi);
355 status |= check_int_to_flt_round(TO_MINUS_INFINITY, -int24_hi, &neg_flt_hi);
356 return status;
357 }
358
359 #ifdef __powerpc64__
360 int test_long_to_double_convert(char *msg)
361 {
362 int status = 0;
363 long long55_hi = 0x07ff0ffffffffff;
364 long long55_lo = 0x07ff0fffffffffd;
365 double pos_dbl_lo = 36012304344547324.0;
366 double pos_dbl_hi = 36012304344547328.0;
367 double neg_dbl_lo = -36012304344547324.0;
368 double neg_dbl_hi = -36012304344547328.0;
369
370 printf("-------------------------- %s --------------------------\n", msg);
371 status |= check_long_to_dbl_round(TO_NEAREST, long55_lo, &pos_dbl_lo);
372 status |= check_long_to_dbl_round(TO_NEAREST, long55_hi, &pos_dbl_hi);
373 status |= check_long_to_dbl_round(TO_ZERO, long55_lo, &pos_dbl_lo);
374 status |= check_long_to_dbl_round(TO_ZERO, long55_hi, &pos_dbl_lo);
375 status |= check_long_to_dbl_round(TO_PLUS_INFINITY, long55_lo, &pos_dbl_hi);
376 status |= check_long_to_dbl_round(TO_PLUS_INFINITY, long55_hi, &pos_dbl_hi);
377 status |= check_long_to_dbl_round(TO_MINUS_INFINITY, long55_lo, &pos_dbl_lo);
378 status |= check_long_to_dbl_round(TO_MINUS_INFINITY, long55_hi, &pos_dbl_lo);
379
380 status |= check_long_to_dbl_round(TO_NEAREST, -long55_lo, &neg_dbl_lo);
381 status |= check_long_to_dbl_round(TO_NEAREST, -long55_hi, &neg_dbl_hi);
382 status |= check_long_to_dbl_round(TO_ZERO, -long55_lo, &neg_dbl_lo);
383 status |= check_long_to_dbl_round(TO_ZERO, -long55_hi, &neg_dbl_lo);
384 status |= check_long_to_dbl_round(TO_PLUS_INFINITY, -long55_lo, &neg_dbl_lo);
385 status |= check_long_to_dbl_round(TO_PLUS_INFINITY, -long55_hi, &neg_dbl_lo);
386 status |= check_long_to_dbl_round(TO_MINUS_INFINITY, -long55_lo, &neg_dbl_hi);
387 status |= check_long_to_dbl_round(TO_MINUS_INFINITY, -long55_hi, &neg_dbl_hi);
388 return status;
389 }
390 #endif
391
392 int check_single_arithmetic_op(flt_op_t op)
393 {
394 char *result;
395 int status = 0;
396 dbl_overlay R, E;
397 double qtr, half, fA, fB, fD;
398 round_mode_t mode;
399 int q, s;
400 bool_t two_args = TRUE;
401 float whole = denorm_small;
402
403 #define BINOP(op) \
404 __asm__ volatile( \
405 op" %0, %1, %2\n\t" \
406 : "=f"(fD) : "f"(fA) , "f"(fB));
407 #define UNOP(op) \
408 __asm__ volatile( \
409 op" %0, %1\n\t" \
410 : "=f"(fD) : "f"(fA));
411
412 half = (double)whole/2;
413 qtr = half/2;
414
415 if (debug) {
416 print_double("qtr", qtr);
417 print_double("whole", whole);
418 print_double("2*whole", 2*whole);
419 }
420
421 for (mode = TO_NEAREST; mode <= TO_MINUS_INFINITY; mode++)
422 for (s = -1; s < 2; s += 2)
423 for (q = 1; q < 4; q += 2) {
424 double expected;
425 double lo = s*whole;
426 double hi = s*2*whole;
427
428 switch(op) {
429 case FADDS:
430 fA = s*whole;
431 fB = s*q*qtr;
432 break;
433 case FSUBS:
434 fA = s*2*whole;
435 fB = s*(q == 1 ? 3 : 1)*qtr;
436 break;
437 case FMULS:
438 fA = 0.5;
439 fB = s*(4+q)*half;
440 break;
441 case FDIVS:
442 fA = s*(4+q)*half;
443 fB = 2.0;
444 break;
445 default:
446 assert("check_single_arithmetic_op: unexpected op",
447 FALSE);
448 break;
449 }
450
451 switch(mode) {
452 case TO_NEAREST:
453 expected = (q == 1 ? lo : hi);
454 break;
455 case TO_ZERO:
456 expected = lo;
457 break;
458 case TO_PLUS_INFINITY:
459 expected = (s == 1 ? hi : lo);
460 break;
461 case TO_MINUS_INFINITY:
462 expected = (s == 1 ? lo : hi);
463 break;
464 }
465
466 set_rounding_mode(mode);
467
468 /*
469 ** do the double precision dual operation just for comparison
470 ** when debugging
471 */
472 switch(op) {
473 case FADDS:
474 BINOP("fadds");
475 R.dbl = fD;
476 BINOP("fadd");
477 break;
478 case FSUBS:
479 BINOP("fsubs");
480 R.dbl = fD;
481 BINOP("fsub");
482 break;
483 case FMULS:
484 BINOP("fmuls");
485 R.dbl = fD;
486 BINOP("fmul");
487 break;
488 case FDIVS:
489 BINOP("fdivs");
490 R.dbl = fD;
491 BINOP("fdiv");
492 break;
493 default:
494 assert("check_single_arithmetic_op: unexpected op",
495 FALSE);
496 break;
497 }
498 #undef UNOP
499 #undef BINOP
500
501 E.dbl = expected;
502
503 if ((R.layout.sign != E.layout.sign) ||
504 (R.layout.exp != E.layout.exp) ||
505 (R.layout.frac_lo != E.layout.frac_lo) ||
506 (R.layout.frac_hi != E.layout.frac_hi)) {
507 result = "FAILED";
508 status = 1;
509 } else {
510 result = "PASSED";
511 status = 0;
512 }
513
514 printf("%s:%s:%s(%-13a",
515 round_mode_name[mode], result, flt_op_names[op], fA);
516 if (two_args) printf(", %-13a", fB);
517 printf(") = %-13a", R.dbl);
518 if (status) printf("\n\texpected %a", E.dbl);
519 putchar('\n');
520
521 if (debug) {
522 print_double("hi", hi);
523 print_double("lo", lo);
524 print_double("expected", expected);
525 print_double("got", R.dbl);
526 print_double("double result", fD);
527 }
528 }
529
530 return status;
531 }
532
533 int check_single_guarded_arithmetic_op(flt_op_t op)
534 {
535 typedef struct {
536 int num, den, frac;
537 } fdivs_t;
538
539 char *result;
540 int status = 0;
541 flt_overlay A, B, Z;
542 dbl_overlay Res, Exp;
543 double fA, fB, fC, fD;
544 round_mode_t mode;
545 int g, s;
546 int arg_count;
547
548 fdivs_t divs_guard_cases[16] = {
549 { 105, 56, 0x700000 }, /* : 0 */
550 { 100, 57, 0x608FB8 }, /* : 1 */
551 { 000, 00, 0x000000 }, /* : X */
552 { 100, 52, 0x762762 }, /* : 3 */
553 { 000, 00, 0x000000 }, /* : X */
554 { 100, 55, 0x68BA2E }, /* : 5 */
555 { 000, 00, 0x000000 }, /* : X */
556 { 100, 51, 0x7AFAFA }, /* : 7 */
557 { 000, 00, 0x000000 }, /* : X */
558 { 100, 56, 0x649249 }, /* : 9 */
559 { 000, 00, 0x000000 }, /* : X */
560 { 100, 54, 0x6D097B }, /* : B */
561 { 000, 00, 0x000000 }, /* : X */
562 { 100, 59, 0x58F2FB }, /* : D */
563 { 000, 00, 0x000000 }, /* : X */
564 { 101, 52, 0x789D89 } /* : F */
565 };
566
567 /* 0x1.00000 00000000p-3 */
568 /* set up the invariant fields of B, the arg to cause rounding */
569 B.flt = 0.0;
570 B.layout.exp = 124; /* -3 */
571
572 /* set up args so result is always Z = 1.200000000000<g>p+0 */
573 Z.flt = 1.0;
574 Z.layout.sign = 0;
575
576 #define TERNOP(op) \
577 arg_count = 3; \
578 __asm__ volatile( \
579 op" %0, %1, %2, %3\n\t" \
580 : "=f"(fD) : "f"(fA) , "f"(fB), "f"(fC));
581 #define BINOP(op) \
582 arg_count = 2; \
583 __asm__ volatile( \
584 op" %0, %1, %2\n\t" \
585 : "=f"(fD) : "f"(fA) , "f"(fB));
586 #define UNOP(op) \
587 arg_count = 1; \
588 __asm__ volatile( \
589 op" %0, %1\n\t" \
590 : "=f"(fD) : "f"(fA));
591
592 for (mode = TO_NEAREST; mode <= TO_MINUS_INFINITY; mode++)
593 for (s = -1; s < 2; s += 2)
594 for (g = 0; g < 16; g += 1) {
595 double lo, hi, expected;
596 int LSB;
597 int guard = 0;
598 int z_sign = s;
599
600 /*
601 ** one argument will have exponent = 0 as will the result (by
602 ** design) so choose the other argument with exponent -3 to
603 ** force a 3 bit shift for scaling leaving us with 3 guard bits
604 ** and the LSB bit at the bottom of the manitssa.
605 */
606 switch(op) {
607 case FADDS:
608 /* 1p+0 + 1.00000<g>p-3 */
609 B.layout.frac = g;
610
611 fB = s*B.flt;
612 fA = s*1.0;
613
614 /* set up Z to be truncated result */
615
616 /* mask off LSB from resulting guard bits */
617 guard = g & 7;
618
619 Z.layout.frac = 0x100000 | (g >> 3);
620 break;
621 case FSUBS:
622 /* 1.200002p+0 - 1.000000000000<g>p-3 */
623 A.flt = 1.125;
624 /* add enough to avoid scaling of the result */
625 A.layout.frac |= 0x2;
626 fA = s*A.flt;
627
628 B.layout.frac = g;
629 fB = s*B.flt;
630
631 /* set up Z to be truncated result */
632 guard = (0x10-g);
633 Z.layout.frac = guard>>3;
634
635 /* mask off LSB from resulting guard bits */
636 guard &= 7;
637 break;
638 case FMULS:
639 /* 1 + g*2^-23 */
640 A.flt = 1.0;
641 A.layout.frac = g;
642 fA = s*A.flt;
643 fB = 1.125;
644
645 /* set up Z to be truncated result */
646 Z.flt = 1.0;
647 Z.layout.frac = 0x100000;
648 Z.layout.frac |= g + (g>>3);
649 guard = g & 7;
650 break;
651 case FDIVS:
652 /* g >> 3 == LSB, g & 7 == guard bits */
653 guard = g & 7;
654 if ((guard & 1) == 0) {
655 /* special case: guard bit X = 0 */
656 A.flt = denorm_small;
657 A.layout.frac = g;
658 fA = A.flt;
659 fB = s*8.0;
660 Z.flt = 0.0;
661 Z.layout.frac |= (g >> 3);
662 } else {
663 fA = s*divs_guard_cases[g].num;
664 fB = divs_guard_cases[g].den;
665
666 Z.flt = 1.0;
667 Z.layout.frac = divs_guard_cases[g].frac;
668 }
669 break;
670 case FMADDS:
671 case FMSUBS:
672 case FNMADDS:
673 case FNMSUBS:
674 /* 1 + g*2^-23 */
675 A.flt = 1.0;
676 A.layout.frac = g;
677 fA = s*A.flt;
678 fB = 1.125;
679
680 /* 1.000001p-1 */
681 A.flt = 0.5;
682 A.layout.frac = 1;
683 fC = (op == FMADDS || op == FNMADDS ? s : -s)*A.flt;
684
685 /* set up Z to be truncated result */
686 z_sign = (op == FNMADDS || op == FNMSUBS ? -s : s);
687 guard = ((g & 7) + 0x4) & 7;
688 Z.flt = 1.0;
689 Z.layout.frac = 0x500000;
690 Z.layout.frac |= g + (g>>3) + ((g & 7)>> 2 ? 1 : 0);
691 break;
692 default:
693 assert("check_single_arithmetic_op: unexpected op",
694 FALSE);
695 break;
696 }
697
698 /* get LSB for tie breaking */
699 LSB = Z.layout.frac & 1;
700
701 /* set up hi and lo */
702 lo = z_sign*Z.flt;
703 Z.layout.frac += 1;
704 hi = z_sign*Z.flt;
705
706 switch(mode) {
707 case TO_NEAREST:
708 /* look at 3 guard bits to determine expected rounding */
709 switch(guard) {
710 case 0:
711 case 1: case 2: case 3:
712 expected = lo;
713 break;
714 case 4: /* tie: round to even */
715 if (debug) printf("tie: LSB = %d\n", LSB);
716 expected = (LSB == 0 ? lo : hi);
717 break;
718 case 5: case 6: case 7:
719 expected = hi;
720 break;
721 default:
722 assert("check_single_guarded_arithmetic_op: unexpected guard",
723 FALSE);
724 }
725 break;
726 case TO_ZERO:
727 expected = lo;
728 break;
729 case TO_PLUS_INFINITY:
730 if (guard == 0) {
731 /* no rounding */
732 expected = lo;
733 } else {
734 expected = (s == 1 ? hi : lo);
735 }
736 break;
737 case TO_MINUS_INFINITY:
738 if (guard == 0) {
739 /* no rounding */
740 expected = lo;
741 } else {
742 expected = (s == 1 ? lo : hi);
743 }
744 break;
745 }
746
747 set_rounding_mode(mode);
748
749 /*
750 ** do the double precision dual operation just for comparison
751 ** when debugging
752 */
753 switch(op) {
754 case FADDS:
755 BINOP("fadds");
756 Res.dbl = fD;
757 break;
758 case FSUBS:
759 BINOP("fsubs");
760 Res.dbl = fD;
761 break;
762 case FMULS:
763 BINOP("fmuls");
764 Res.dbl = fD;
765 break;
766 case FDIVS:
767 BINOP("fdivs");
768 Res.dbl = fD;
769 break;
770 case FMADDS:
771 TERNOP("fmadds");
772 Res.dbl = fD;
773 break;
774 case FMSUBS:
775 TERNOP("fmsubs");
776 Res.dbl = fD;
777 break;
778 case FNMADDS:
779 TERNOP("fnmadds");
780 Res.dbl = fD;
781 break;
782 case FNMSUBS:
783 TERNOP("fnmsubs");
784 Res.dbl = fD;
785 break;
786 default:
787 assert("check_single_guarded_arithmetic_op: unexpected op",
788 FALSE);
789 break;
790 }
791 #undef UNOP
792 #undef BINOP
793 #undef TERNOP
794
795 Exp.dbl = expected;
796
797 if ((Res.layout.sign != Exp.layout.sign) ||
798 (Res.layout.exp != Exp.layout.exp) ||
799 (Res.layout.frac_lo != Exp.layout.frac_lo) ||
800 (Res.layout.frac_hi != Exp.layout.frac_hi)) {
801 result = "FAILED";
802 status = 1;
803 } else {
804 result = "PASSED";
805 status = 0;
806 }
807
808 /* There seems to be some noise in the lower bits. The value
809 * on the least significant digit seems to vary when printing
810 * based on the rounding mode of the compiler. Just trying
811 * to get rid of the noise in the least significant bits when
812 * printing the operand.
813 */
814
815 fA = ((long int)(fA*10000))/10000.0;
816 /* Change -0.0 to a positive 0.0. Some compilers print -0.0
817 * others do not. Make it consistent.
818 */
819 if (fA == -0.0)
820 fA = 0.0;
821
822 printf("%s:%s:%s(%-13.6f",
823 round_mode_name[mode], result, flt_op_names[op], fA);
824 if (arg_count > 1) printf(", %-13a", fB);
825 if (arg_count > 2) printf(", %-13a", fC);
826 printf(") = %-13a", Res.dbl);
827 if (status) printf("\n\texpected %a", Exp.dbl);
828 putchar('\n');
829
830 if (debug) {
831 print_double("hi", hi);
832 print_double("lo", lo);
833 print_double("expected", expected);
834 print_double("got", Res.dbl);
835 }
836 }
837
838 return status;
839 }
840
841 int check_double_guarded_arithmetic_op(flt_op_t op)
842 {
843 typedef struct {
844 int num, den, hi, lo;
845 } fdiv_t;
846 typedef struct {
847 double arg;
848 int exp, hi, lo;
849 } fsqrt_t;
850
851 char *result;
852 int status = 0;
853 dbl_overlay A, B, Z;
854 dbl_overlay Res, Exp;
855 double fA, fB, fC, fD;
856 round_mode_t mode;
857 int g, s;
858 int arg_count;
859 fdiv_t div_guard_cases[16] = {
860 { 62, 62, 0x00000, 0x00000000 }, /* 0 */
861 { 64, 62, 0x08421, 0x08421084 }, /* 1 */
862 { 66, 62, 0x10842, 0x10842108 }, /* 2 */
863 { 100, 62, 0x9ce73, 0x9ce739ce }, /* 3 */
864 { 100, 62, 0x9ce73, 0x9ce739ce }, /* X */
865 { 102, 62, 0xa5294, 0xa5294a52 }, /* 5 */
866 { 106, 62, 0xb5ad6, 0xb5ad6b5a }, /* 6 */
867 { 108, 62, 0xbdef7, 0xbdef7bde }, /* 7 */
868 { 108, 108, 0x00000, 0x00000000 }, /* 8 */
869 { 112, 62, 0xce739, 0xce739ce7 }, /* 9 */
870 { 114, 62, 0xd6b5a, 0xd6b5ad6b }, /* A */
871 { 116, 62, 0xdef7b, 0xdef7bdef }, /* B */
872 { 84, 62, 0x5ad6b, 0x5ad6b5ad }, /* X */
873 { 118, 62, 0xe739c, 0xe739ce73 }, /* D */
874 { 90, 62, 0x739ce, 0x739ce739 }, /* E */
875 { 92, 62, 0x7bdef, 0x7bdef7bd } /* F */
876 };
877
878
879 fsqrt_t sqrt_guard_cases[16] = {
880 { 0x1.08800p0, 0, 0x04371, 0xd9ab72fb}, /* :0 B8.8440 */
881 { 0x0.D2200p0, -1, 0xcfdca, 0xf353049e}, /* :1 A4.6910 */
882 { 0x1.A8220p0, 0, 0x49830, 0x2b49cd6d}, /* :2 E9.D411 */
883 { 0x1.05A20p0, 0, 0x02cd1, 0x3b44f3bf}, /* :3 B7.82D1 */
884 { 0x0.CA820p0, -1, 0xc7607, 0x3cec0937}, /* :4 A1.6541 */
885 { 0x1.DCA20p0, 0, 0x5d4f8, 0xd4e4c2b2}, /* :5 F7.EE51 */
886 { 0x1.02C80p0, 0, 0x01630, 0x9cde7483}, /* :6 B6.8164 */
887 { 0x0.DC800p0, -1, 0xdb2cf, 0xe686fe7c}, /* :7 A8.6E40 */
888 { 0x0.CF920p0, -1, 0xcd089, 0xb6860626}, /* :8 A3.67C9 */
889 { 0x1.1D020p0, 0, 0x0e1d6, 0x2e78ed9d}, /* :9 BF.8E81 */
890 { 0x0.E1C80p0, -1, 0xe0d52, 0x6020fb6b}, /* :A AA.70E4 */
891 { 0x0.C8000p0, -1, 0xc48c6, 0x001f0abf}, /* :B A0.6400 */
892 { 0x1.48520p0, 0, 0x21e9e, 0xd813e2e2}, /* :C CD.A429 */
893 { 0x0.F4C20p0, -1, 0xf4a1b, 0x09bbf0b0}, /* :D B1.7A61 */
894 { 0x0.CD080p0, -1, 0xca348, 0x79b907ae}, /* :E A2.6684 */
895 { 0x1.76B20p0, 0, 0x35b67, 0x81aed827} /* :F DB.BB59 */
896 };
897
898 /* 0x1.00000 00000000p-3 */
899 /* set up the invariant fields of B, the arg to cause rounding */
900 B.dbl = 0.0;
901 B.layout.exp = 1020;
902
903 /* set up args so result is always Z = 1.200000000000<g>p+0 */
904 Z.dbl = 1.0;
905 Z.layout.sign = 0;
906
907 #define TERNOP(op) \
908 arg_count = 3; \
909 __asm__ volatile( \
910 op" %0, %1, %2, %3\n\t" \
911 : "=f"(fD) : "f"(fA) , "f"(fB), "f"(fC));
912 #define BINOP(op) \
913 arg_count = 2; \
914 __asm__ volatile( \
915 op" %0, %1, %2\n\t" \
916 : "=f"(fD) : "f"(fA) , "f"(fB));
917 #define UNOP(op) \
918 arg_count = 1; \
919 __asm__ volatile( \
920 op" %0, %1\n\t" \
921 : "=f"(fD) : "f"(fA));
922
923 for (mode = TO_NEAREST; mode <= TO_MINUS_INFINITY; mode++)
924 for (s = (op != FSQRT ? -1 : 1); s < 2; s += 2)
925 for (g = 0; g < 16; g += 1) {
926 double lo, hi, expected;
927 int LSB;
928 int guard;
929 int z_sign = s;
930
931 /*
932 ** one argument will have exponent = 0 as will the result (by
933 ** design) so choose the other argument with exponent -3 to
934 ** force a 3 bit shift for scaling leaving us with 3 guard bits
935 ** and the LSB bit at the bottom of the manitssa.
936 */
937 switch(op) {
938 case FADD:
939 /* 1p+0 + 1.000000000000<g>p-3 */
940 B.layout.frac_lo = g;
941
942 fB = s*B.dbl;
943 fA = s*1.0;
944
945 /* set up Z to be truncated result */
946
947 /* mask off LSB from resulting guard bits */
948 guard = g & 7;
949
950 Z.layout.frac_hi = 0x20000;
951 Z.layout.frac_lo = g >> 3;
952
953 break;
954 case FSUB:
955 /* 1.2000000000002p+0 - 1.000000000000<g>p-3 */
956 A.dbl = 1.125;
957 /* add enough to avoid scaling of the result */
958 A.layout.frac_lo = 0x2;
959 fA = s*A.dbl;
960
961 B.layout.frac_lo = g;
962 fB = s*B.dbl;
963
964 /* set up Z to be truncated result */
965 guard = (0x10-g);
966 Z.layout.frac_hi = 0x0;
967 Z.layout.frac_lo = guard>>3;
968
969 /* mask off LSB from resulting guard bits */
970 guard &= 7;
971 break;
972 case FMUL:
973 /* 1 + g*2^-52 */
974 A.dbl = 1.0;
975 A.layout.frac_lo = g;
976 fA = s*A.dbl;
977 fB = 1.125;
978
979 /* set up Z to be truncated result */
980 Z.dbl = 1.0;
981 Z.layout.frac_hi = 0x20000;
982 Z.layout.frac_lo = g + (g>>3);
983 guard = g & 7;
984 break;
985 case FMADD:
986 case FMSUB:
987 case FNMADD:
988 case FNMSUB:
989 /* 1 + g*2^-52 */
990 A.dbl = 1.0;
991 A.layout.frac_lo = g;
992 fA = s*A.dbl;
993 fB = 1.125;
994
995 /* 1.0000000000001p-1 */
996 A.dbl = 0.5;
997 A.layout.frac_lo = 1;
998 fC = (op == FMADD || op == FNMADD ? s : -s)*A.dbl;
999
1000 /* set up Z to be truncated result */
1001 z_sign = (op == FNMADD || op == FNMSUB ? -s : s);
1002 guard = ((g & 7) + 0x4) & 7;
1003 Z.dbl = 1.0;
1004 Z.layout.frac_hi = 0xa0000;
1005 Z.layout.frac_lo = g + (g>>3) + ((g & 7)>> 2 ? 1 : 0);
1006 break;
1007 case FDIV:
1008 /* g >> 3 == LSB, g & 7 == guard bits */
1009 guard = g & 7;
1010 if (guard == 0x4) {
1011 /* special case guard bits == 4, inexact tie */
1012 fB = s*2.0;
1013 Z.dbl = 0.0;
1014 if (g >> 3) {
1015 fA = dbl_denorm_small + 2*dbl_denorm_small;
1016 Z.layout.frac_lo = 0x1;
1017 } else {
1018 fA = dbl_denorm_small;
1019 }
1020 } else {
1021 fA = s*div_guard_cases[g].num;
1022 fB = div_guard_cases[g].den;
1023
1024 printf("%d/%d\n",
1025 s*div_guard_cases[g].num,
1026 div_guard_cases[g].den);
1027 Z.dbl = 1.0;
1028 Z.layout.frac_hi = div_guard_cases[g].hi;
1029 Z.layout.frac_lo = div_guard_cases[g].lo;
1030 }
1031 break;
1032 case FSQRT:
1033 fA = s*sqrt_guard_cases[g].arg;
1034 Z.dbl = 1.0;
1035 Z.layout.exp = sqrt_guard_cases[g].exp + 1023;
1036 Z.layout.frac_hi = sqrt_guard_cases[g].hi;
1037 Z.layout.frac_lo = sqrt_guard_cases[g].lo;
1038 guard = g >> 1;
1039 if (g & 1) guard |= 1;
1040 /* don't have test cases for when X bit = 0 */
1041 if (guard == 0 || guard == 4) continue;
1042 break;
1043 default:
1044 assert("check_double_guarded_arithmetic_op: unexpected op",
1045 FALSE);
1046 break;
1047 }
1048
1049 /* get LSB for tie breaking */
1050 LSB = Z.layout.frac_lo & 1;
1051
1052 /* set up hi and lo */
1053 lo = z_sign*Z.dbl;
1054 Z.layout.frac_lo += 1;
1055 hi = z_sign*Z.dbl;
1056
1057 switch(mode) {
1058 case TO_NEAREST:
1059 /* look at 3 guard bits to determine expected rounding */
1060 switch(guard) {
1061 case 0:
1062 case 1: case 2: case 3:
1063 expected = lo;
1064 break;
1065 case 4: /* tie: round to even */
1066 if (debug) printf("tie: LSB = %d\n", LSB);
1067 expected = (LSB == 0 ? lo : hi);
1068 break;
1069 case 5: case 6: case 7:
1070 expected = hi;
1071 break;
1072 default:
1073 assert("check_double_guarded_arithmetic_op: unexpected guard",
1074 FALSE);
1075 }
1076 break;
1077 case TO_ZERO:
1078 expected = lo;
1079 break;
1080 case TO_PLUS_INFINITY:
1081 if (guard == 0) {
1082 /* no rounding */
1083 expected = lo;
1084 } else {
1085 expected = (s == 1 ? hi : lo);
1086 }
1087 break;
1088 case TO_MINUS_INFINITY:
1089 if (guard == 0) {
1090 /* no rounding */
1091 expected = lo;
1092 } else {
1093 expected = (s == 1 ? lo : hi);
1094 }
1095 break;
1096 }
1097
1098 set_rounding_mode(mode);
1099
1100 /*
1101 ** do the double precision dual operation just for comparison
1102 ** when debugging
1103 */
1104 switch(op) {
1105 case FADD:
1106 BINOP("fadd");
1107 Res.dbl = fD;
1108 break;
1109 case FSUB:
1110 BINOP("fsub");
1111 Res.dbl = fD;
1112 break;
1113 case FMUL:
1114 BINOP("fmul");
1115 Res.dbl = fD;
1116 break;
1117 case FMADD:
1118 TERNOP("fmadd");
1119 Res.dbl = fD;
1120 break;
1121 case FMSUB:
1122 TERNOP("fmsub");
1123 Res.dbl = fD;
1124 break;
1125 case FNMADD:
1126 TERNOP("fnmadd");
1127 Res.dbl = fD;
1128 break;
1129 case FNMSUB:
1130 TERNOP("fnmsub");
1131 Res.dbl = fD;
1132 break;
1133 case FDIV:
1134 BINOP("fdiv");
1135 Res.dbl = fD;
1136 break;
1137 case FSQRT:
1138 UNOP("fsqrt");
1139 Res.dbl = fD;
1140 break;
1141 default:
1142 assert("check_double_guarded_arithmetic_op: unexpected op",
1143 FALSE);
1144 break;
1145 }
1146 #undef UNOP
1147 #undef BINOP
1148 #undef TERNOP
1149
1150 Exp.dbl = expected;
1151
1152 if ((Res.layout.sign != Exp.layout.sign) ||
1153 (Res.layout.exp != Exp.layout.exp) ||
1154 (Res.layout.frac_lo != Exp.layout.frac_lo) ||
1155 (Res.layout.frac_hi != Exp.layout.frac_hi)) {
1156 result = "FAILED";
1157 status = 1;
1158 } else {
1159 result = "PASSED";
1160 status = 0;
1161 }
1162
1163 printf("%s:%s:%s(%-13a",
1164 round_mode_name[mode], result, flt_op_names[op], fA);
1165 if (arg_count > 1) printf(", %-13a", fB);
1166 if (arg_count > 2) printf(", %-13a", fC);
1167 printf(") = %-13a", Res.dbl);
1168 if (status) printf("\n\texpected %a", Exp.dbl);
1169 putchar('\n');
1170
1171 if (debug) {
1172 print_double("hi", hi);
1173 print_double("lo", lo);
1174 print_double("expected", expected);
1175 print_double("got", Res.dbl);
1176 }
1177 }
1178
1179 return status;
1180 }
1181
1182 int test_float_arithmetic_ops()
1183 {
1184 int status = 0;
1185 flt_op_t op;
1186
1187 /*
1188 ** choose FP operands whose result should be rounded to either
1189 ** lo or hi.
1190 */
1191
1192 printf("-------------------------- %s --------------------------\n",
1193 "test rounding of float operators without guard bits");
1194 for (op = FADDS; op <= FDIVS; op++) {
1195 status |= check_single_arithmetic_op(op);
1196 }
1197
1198 printf("-------------------------- %s --------------------------\n",
1199 "test rounding of float operators with guard bits");
1200 for (op = FADDS; op <= FNMSUBS; op++) {
1201 status |= check_single_guarded_arithmetic_op(op);
1202 }
1203
1204 printf("-------------------------- %s --------------------------\n",
1205 "test rounding of double operators with guard bits");
1206 for (op = FADD; op <= FSQRT; op++) {
1207 status |= check_double_guarded_arithmetic_op(op);
1208 }
1209 return status;
1210 }
1211
1212
1213 int
1214 main()
1215 {
1216 int status = 0;
1217
1218 init();
1219
1220 status |= test_dbl_to_float_convert("test denormalized convert", &denorm_small);
1221 status |= test_dbl_to_float_convert("test normalized convert", &norm_small);
1222 status |= test_int_to_float_convert("test (float)int convert");
1223 status |= test_int_to_float_convert("test (float)int convert");
1224
1225 #ifdef __powerpc64__
1226 status |= test_long_to_double_convert("test (double)long convert");
1227 #endif
1228 status |= test_float_arithmetic_ops();
1229 return status;
1230 }
Instrumentation framework for detecting many memory management and threading bugs, and profile your programs in detail.