sysdeps/powerpc/test-arith.c

   1 /* Test floating-point arithmetic operations.
   2    Copyright (C) 1997, 1998 Free Software Foundation, Inc.
   3    This file is part of the GNU C Library.
   4
   5    The GNU C Library is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU Library 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    The GNU C Library is distributed in the hope that it will be useful,
  11    but WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13    Library General Public License for more details.
  14
  15    You should have received a copy of the GNU Library General Public
  16    License along with the GNU C Library; see the file COPYING.LIB.  If not,
  17    write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  18    Boston, MA 02111-1307, USA.  */
  19 #ifndef _GNU_SOURCE
  20 #define _GNU_SOURCE
  21 #endif
  22 #include <math.h>
  23 #include <stdio.h>
  24 #include <stdlib.h>
  25 #include <string.h>
  26 #include <fenv.h>
  27 #include <assert.h>
  28
  29 #ifndef ESIZE
  30 typedef double tocheck_t;
  31 #define ESIZE 11
  32 #define MSIZE 52
  33 #define FUNC(x) x
  34 #endif
  35
  36 #define R_NEAREST 1
  37 #define R_ZERO 2
  38 #define R_UP 4
  39 #define R_DOWN 8
  40 #define R_ALL (R_NEAREST|R_ZERO|R_UP|R_DOWN)
  41 static fenv_t rmodes[4];
  42 static const char * const rmnames[4] =
  43 { "nearest","zero","+Inf","-Inf" };
  44
  45 typedef union {
  46   tocheck_t tc;
  47   unsigned char c[sizeof(tocheck_t)];
  48 } union_t;
  49
  50 /* Don't try reading these in a font that doesn't distinguish
  51    O and zero.  */
  52 typedef enum {
  53   P_Z    = 0x0,  /* 00000...0 */
  54   P_000O = 0x1,  /* 00011...1 */
  55   P_001Z = 0x2,  /* 00100...0 */
  56   P_00O  = 0x3,  /* 00111...1 */
  57   P_01Z  = 0x4,  /* 01000...0 */
  58   P_010O = 0x5,  /* 01011...1 */
  59   P_011Z = 0x6,  /* 01100...0 */
  60   P_0O   = 0x7,  /* 01111...1 */
  61   P_1Z   = 0x8,  /* 10000...0 */
  62   P_100O = 0x9,  /* 10011...1 */
  63   P_101Z = 0xa,  /* 10100...0 */
  64   P_10O  = 0xb,  /* 10111...1 */
  65   P_11Z  = 0xc,  /* 11000...0 */
  66   P_110O = 0xd,  /* 11011...1 */
  67   P_111Z = 0xe,  /* 11100...0 */
  68   P_O    = 0xf,  /* 11111...1 */
  69   P_Z1   = 0x11, /* 000...001 */
  70   P_Z10  = 0x12, /* 000...010 */
  71   P_Z11  = 0x13, /* 000...011 */
  72   P_0O00 = 0x14, /* 011...100 */
  73   P_0O01 = 0x15, /* 011...101 */
  74   P_0O0  = 0x16, /* 011...110 */
  75   P_1Z1  = 0x19, /* 100...001 */
  76   P_1Z10 = 0x1a, /* 100...010 */
  77   P_1Z11 = 0x1b, /* 100...011 */
  78   P_O00  = 0x1c, /* 111...100 */
  79   P_O01  = 0x1d, /* 111...101 */
  80   P_O0   = 0x1e, /* 111...110 */
  81   P_R    = 0x20, /* rrr...rrr */ /* ('r' means random. ) */
  82   P_Ro   = 0x21, /* rrr...rrr, with odd parity.  */
  83   P_0R   = 0x22, /* 0rr...rrr */
  84   P_1R   = 0x23, /* 1rr...rrr */
  85   P_Rno  = 0x24, /* rrr...rrr, but not all ones.  */
  86 } pattern_t;
  87
  88 static void
  89 pattern_fill(pattern_t ptn, unsigned char *start, int bitoffset, int count)
  90 {
  91 #define bitset(count, value) \
  92       start[(count)/8] = (start[(count)/8] & ~(1 << 7-(count)%8)  \
  93                           |  (value) << 7-(count)%8)
  94   int i;
  95
  96   if (ptn >= 0 && ptn <= 0xf)
  97     {
  98       /* Patterns between 0 and 0xF have the following format:
  99          The LSBit is used to fill the last n-3 bits of the pattern;
 100          The next 3 bits are the first 3 bits of the pattern. */
 101       for (i = 0; i < count; i++)
 102         if (i < 3)
 103           bitset((bitoffset+i), ptn >> (3-i) & 1);
 104         else
 105           bitset((bitoffset+i), ptn >> 0 & 1);
 106     }
 107   else if (ptn <= 0x1f)
 108     {
 109       /* Patterns between 0x10 and 0x1F have the following format:
 110          The two LSBits are the last two bits of the pattern;
 111          The 0x8 bit is the first bit of the pattern;
 112          The 0x4 bit is used to fill the remainder. */
 113       for (i = 0; i < count; i++)
 114         if (i == 0)
 115           bitset((bitoffset+i), ptn >> 3 & 1);
 116         else if (i >= count-2)
 117           bitset((bitoffset+i), ptn >> (count-1-i) & 1);
 118         else
 119           bitset((bitoffset+i), ptn >> 2 & 1);
 120     }
 121   else switch (ptn)
 122     {
 123     case P_0R: case P_1R:
 124       assert(count > 0);
 125       bitset(bitoffset, ptn & 1);
 126       count--;
 127       bitoffset++;
 128     case P_R:
 129       for (; count > 0; count--, bitoffset++)
 130         bitset(bitoffset, rand() & 1);
 131       break;
 132     case P_Ro:
 133       {
 134         int op = 1;
 135         assert(count > 0);
 136         for (; count > 1; count--, bitoffset++)
 137           bitset(bitoffset, op ^= (rand() & 1));
 138         bitset(bitoffset, op);
 139         break;
 140       }
 141     case P_Rno:
 142       {
 143         int op = 1;
 144         assert(count > 0);
 145         for (; count > 1; count--, bitoffset++)
 146         {
 147           int r = rand() & 1;
 148           op &= r;
 149           bitset(bitoffset, r);
 150         }
 151         bitset(bitoffset, rand() & (op ^ 1));
 152         break;
 153       }
 154
 155     default:
 156       assert(0);
 157     }
 158 #undef bitset
 159 }
 160
 161 static tocheck_t
 162 pattern(int negative, pattern_t exp, pattern_t mant)
 163 {
 164   union_t result;
 165 #if 0
 166   int i;
 167 #endif
 168
 169   pattern_fill(negative ? P_O : P_Z, result.c, 0, 1);
 170   pattern_fill(exp, result.c, 1, ESIZE);
 171   pattern_fill(mant, result.c, ESIZE+1, MSIZE);
 172 #if 0
 173   printf("neg=%d exp=%02x mant=%02x: ", negative, exp, mant);
 174   for (i = 0; i < sizeof(tocheck_t); i++)
 175     printf("%02x", result.c[i]);
 176   printf("\n");
 177 #endif
 178   return result.tc;
 179 }
 180
 181 /* Return the closest different tocheck_t to 'x' in the direction of
 182    'direction', or 'x' if there is no such value.  Assumes 'x' is not
 183    a NaN.  */
 184 static tocheck_t
 185 delta(tocheck_t x, int direction)
 186 {
 187   union_t xx;
 188   int i;
 189
 190   xx.tc = x;
 191   if (xx.c[0] & 0x80)
 192     direction = -direction;
 193   if (direction == +1)
 194     {
 195       union_t tx;
 196       tx.tc = pattern(xx.c[0] >> 7, P_O, P_Z);
 197       if (memcmp(tx.c, xx.c, sizeof(tocheck_t)) == 0)
 198         return x;
 199     }
 200   for (i = sizeof(tocheck_t)-1; i > 0; i--)
 201     {
 202       xx.c[i] += direction;
 203       if (xx.c[i] != (direction > 0 ? 0 : 0xff))
 204         return xx.tc;
 205     }
 206   if (direction < 0 && (xx.c[0] & 0x7f) == 0)
 207     return pattern(~(xx.c[0] >> 7) & 1, P_Z, P_Z1);
 208   else
 209     {
 210       xx.c[0] += direction;
 211       return xx.tc;
 212     }
 213 }
 214
 215 static int nerrors = 0;
 216
 217 #ifdef FE_ALL_INVALID
 218 static const int all_exceptions = FE_ALL_INVALID | FE_ALL_EXCEPT;
 219 #else
 220 static const int all_exceptions = FE_ALL_EXCEPT;
 221 #endif
 222
 223 static void
 224 check_result(int line, const char *rm, tocheck_t expected, tocheck_t actual)
 225 {
 226   if (memcmp(&expected, &actual, sizeof(tocheck_t)) != 0)
 227     {
 228       unsigned char *ex, *ac;
 229       size_t i;
 230
 231       printf("%s:%d:round %s:result failed\n"
 232              " expected result 0x", __FILE__, line, rm);
 233       ex = (unsigned char *)&expected;
 234       ac = (unsigned char *)&actual;
 235       for (i = 0; i < sizeof(tocheck_t); i++)
 236         printf("%02x", ex[i]);
 237       printf(" got 0x");
 238       for (i = 0; i < sizeof(tocheck_t); i++)
 239         printf("%02x", ac[i]);
 240       printf("\n");
 241       nerrors++;
 242     }
 243 }
 244
 245 static const struct {
 246   int except;
 247   const char *name;
 248 } excepts[] = {
 249 #define except_entry(ex) { ex, #ex } ,
 250 #ifdef FE_INEXACT
 251   except_entry(FE_INEXACT)
 252 #else
 253 # define FE_INEXACT 0
 254 #endif
 255 #ifdef FE_DIVBYZERO
 256   except_entry(FE_DIVBYZERO)
 257 #else
 258 # define FE_DIVBYZERO 0
 259 #endif
 260 #ifdef FE_UNDERFLOW
 261   except_entry(FE_UNDERFLOW)
 262 #else
 263 # define FE_UNDERFLOW 0
 264 #endif
 265 #ifdef FE_OVERFLOW
 266   except_entry(FE_OVERFLOW)
 267 #else
 268 # define FE_OVERFLOW 0
 269 #endif
 270 #ifdef FE_INVALID
 271   except_entry(FE_INVALID)
 272 #else
 273 # define FE_INVALID 0
 274 #endif
 275 #ifdef FE_INVALID_SNAN
 276   except_entry(FE_INVALID_SNAN)
 277 #else
 278 # define FE_INVALID_SNAN FE_INVALID
 279 #endif
 280 #ifdef FE_INVALID_ISI
 281   except_entry(FE_INVALID_ISI)
 282 #else
 283 # define FE_INVALID_ISI FE_INVALID
 284 #endif
 285 #ifdef FE_INVALID_IDI
 286   except_entry(FE_INVALID_IDI)
 287 #else
 288 # define FE_INVALID_IDI FE_INVALID
 289 #endif
 290 #ifdef FE_INVALID_ZDZ
 291   except_entry(FE_INVALID_ZDZ)
 292 #else
 293 # define FE_INVALID_ZDZ FE_INVALID
 294 #endif
 295 #ifdef FE_INVALID_COMPARE
 296   except_entry(FE_INVALID_COMPARE)
 297 #else
 298 # define FE_INVALID_COMPARE FE_INVALID
 299 #endif
 300 #ifdef FE_INVALID_SOFTWARE
 301   except_entry(FE_INVALID_SOFTWARE)
 302 #else
 303 # define FE_INVALID_SOFTWARE FE_INVALID
 304 #endif
 305 #ifdef FE_INVALID_SQRT
 306   except_entry(FE_INVALID_SQRT)
 307 #else
 308 # define FE_INVALID_SQRT FE_INVALID
 309 #endif
 310 #ifdef FE_INVALID_INTEGER_CONVERSION
 311   except_entry(FE_INVALID_INTEGER_CONVERSION)
 312 #else
 313 # define FE_INVALID_INTEGER_CONVERSION FE_INVALID
 314 #endif
 315 };
 316
 317 static int excepts_missing = 0;
 318
 319 static void
 320 check_excepts(int line, const char *rm, int expected, int actual)
 321 {
 322   if (expected & excepts_missing)
 323     expected = expected & ~excepts_missing | FE_INVALID_SNAN;
 324   if ((expected & all_exceptions) != actual)
 325     {
 326       size_t i;
 327       printf("%s:%d:round %s:exceptions failed\n"
 328              " expected exceptions ", __FILE__, line,rm);
 329       for (i = 0; i < sizeof(excepts)/sizeof(excepts[0]); i++)
 330         if (expected & excepts[i].except)
 331           printf("%s ",excepts[i].name);
 332       if ((expected & all_exceptions) == 0)
 333         printf("- ");
 334       printf("got");
 335       for (i = 0; i < sizeof(excepts)/sizeof(excepts[0]); i++)
 336         if (actual & excepts[i].except)
 337           printf(" %s",excepts[i].name);
 338       if ((actual & all_exceptions) == 0)
 339         printf("- ");
 340       printf(".\n");
 341       nerrors++;
 342     }
 343 }
 344
 345 typedef enum {
 346   B_ADD, B_SUB, B_MUL, B_DIV, B_NEG, B_ABS, B_SQRT
 347 } op_t;
 348 typedef struct {
 349   int line;
 350   op_t op;
 351   int a_sgn;
 352   pattern_t a_exp, a_mant;
 353   int b_sgn;
 354   pattern_t b_exp, b_mant;
 355   int rmode;
 356   int excepts;
 357   int x_sgn;
 358   pattern_t x_exp, x_mant;
 359 } optest_t;
 360 static const optest_t optests[] = {
 361   /* Additions of zero.  */
 362   {__LINE__,B_ADD, 0,P_Z,P_Z, 0,P_Z,P_Z, R_ALL,0, 0,P_Z,P_Z },
 363   {__LINE__,B_ADD, 1,P_Z,P_Z, 0,P_Z,P_Z, R_ALL & ~R_DOWN,0, 0,P_Z,P_Z },
 364   {__LINE__,B_ADD, 1,P_Z,P_Z, 0,P_Z,P_Z, R_DOWN,0, 1,P_Z,P_Z },
 365   {__LINE__,B_ADD, 1,P_Z,P_Z, 1,P_Z,P_Z, R_ALL,0, 1,P_Z,P_Z },
 366
 367   /* Additions with NaN.  */
 368   {__LINE__,B_ADD, 0,P_O,P_101Z, 0,P_Z,P_Z, R_ALL,0, 0,P_O,P_101Z },
 369   {__LINE__,B_ADD, 0,P_O,P_01Z, 0,P_Z,P_Z, R_ALL,
 370    FE_INVALID | FE_INVALID_SNAN, 0,P_O,P_11Z },
 371   {__LINE__,B_ADD, 0,P_O,P_Z, 0,P_O,P_0O, R_ALL,
 372    FE_INVALID | FE_INVALID_SNAN, 0,P_O,P_O },
 373   {__LINE__,B_ADD, 0,P_Z,P_Z, 0,P_O,P_11Z, R_ALL,0, 0,P_O,P_11Z },
 374   {__LINE__,B_ADD, 0,P_O,P_001Z, 0,P_O,P_001Z, R_ALL,
 375    FE_INVALID | FE_INVALID_SNAN, 0,P_O,P_101Z },
 376   {__LINE__,B_ADD, 0,P_O,P_1Z, 0,P_Z,P_Z, R_ALL,0, 0,P_O,P_1Z },
 377   {__LINE__,B_ADD, 0,P_0O,P_Z, 0,P_O,P_10O, R_ALL,0, 0,P_O,P_10O },
 378
 379   /* Additions with infinity.  */
 380   {__LINE__,B_ADD, 0,P_O,P_Z, 0,P_Z,P_Z, R_ALL,0, 0,P_O,P_Z },
 381   {__LINE__,B_ADD, 0,P_O,P_Z, 1,P_Z,P_Z, R_ALL,0, 0,P_O,P_Z },
 382   {__LINE__,B_ADD, 1,P_O,P_Z, 0,P_Z,P_Z, R_ALL,0, 1,P_O,P_Z },
 383   {__LINE__,B_ADD, 1,P_O,P_Z, 1,P_Z,P_Z, R_ALL,0, 1,P_O,P_Z },
 384   {__LINE__,B_ADD, 0,P_O,P_Z, 0,P_O,P_Z, R_ALL,0, 0,P_O,P_Z },
 385   {__LINE__,B_ADD, 1,P_O,P_Z, 1,P_O,P_Z, R_ALL,0, 1,P_O,P_Z },
 386   {__LINE__,B_ADD, 0,P_O,P_Z, 1,P_O,P_Z, R_ALL,
 387    FE_INVALID | FE_INVALID_ISI, 0,P_O,P_1Z },
 388   {__LINE__,B_ADD, 1,P_O,P_Z, 0,P_O,P_Z, R_ALL,
 389    FE_INVALID | FE_INVALID_ISI, 0,P_O,P_1Z },
 390   {__LINE__,B_ADD, 0,P_O,P_Z, 0,P_0O,P_Z, R_ALL,0, 0,P_O,P_Z },
 391   {__LINE__,B_ADD, 1,P_O,P_Z, 0,P_0O,P_Z, R_ALL,0, 1,P_O,P_Z },
 392   {__LINE__,B_ADD, 0,P_O,P_Z, 1,P_0O,P_Z, R_ALL,0, 0,P_O,P_Z },
 393   {__LINE__,B_ADD, 1,P_O,P_Z, 1,P_0O,P_Z, R_ALL,0, 1,P_O,P_Z },
 394
 395   /* Overflow (and zero).  */
 396   {__LINE__,B_ADD, 0,P_O0,P_Z, 0,P_O0,P_Z, R_NEAREST | R_UP,
 397    FE_INEXACT | FE_OVERFLOW, 0,P_O,P_Z },
 398   {__LINE__,B_ADD, 0,P_O0,P_Z, 0,P_O0,P_Z, R_ZERO | R_DOWN,
 399    FE_INEXACT | FE_OVERFLOW, 0,P_O0,P_O },
 400   {__LINE__,B_ADD, 1,P_O0,P_Z, 1,P_O0,P_Z, R_NEAREST | R_DOWN,
 401    FE_INEXACT | FE_OVERFLOW, 1,P_O,P_Z },
 402   {__LINE__,B_ADD, 1,P_O0,P_Z, 1,P_O0,P_Z, R_ZERO | R_UP,
 403    FE_INEXACT | FE_OVERFLOW, 1,P_O0,P_O },
 404   {__LINE__,B_ADD, 0,P_O0,P_Z, 1,P_O0,P_Z, R_ALL & ~R_DOWN,
 405    0, 0,P_Z,P_Z },
 406   {__LINE__,B_ADD, 0,P_O0,P_Z, 1,P_O0,P_Z, R_DOWN,
 407    0, 1,P_Z,P_Z },
 408
 409   /* Negation.  */
 410   {__LINE__,B_NEG, 0,P_Z,P_Z,   0,0,0, R_ALL, 0, 1,P_Z,P_Z },
 411   {__LINE__,B_NEG, 1,P_Z,P_Z,   0,0,0, R_ALL, 0, 0,P_Z,P_Z },
 412   {__LINE__,B_NEG, 0,P_O,P_Z,   0,0,0, R_ALL, 0, 1,P_O,P_Z },
 413   {__LINE__,B_NEG, 1,P_O,P_Z,   0,0,0, R_ALL, 0, 0,P_O,P_Z },
 414   {__LINE__,B_NEG, 0,P_O,P_1Z,  0,0,0, R_ALL, 0, 1,P_O,P_1Z },
 415   {__LINE__,B_NEG, 1,P_O,P_1Z,  0,0,0, R_ALL, 0, 0,P_O,P_1Z },
 416   {__LINE__,B_NEG, 0,P_O,P_01Z, 0,0,0, R_ALL, 0, 1,P_O,P_01Z },
 417   {__LINE__,B_NEG, 1,P_O,P_01Z, 0,0,0, R_ALL, 0, 0,P_O,P_01Z },
 418   {__LINE__,B_NEG, 0,P_1Z,P_1Z1, 0,0,0, R_ALL, 0, 1,P_1Z,P_1Z1 },
 419   {__LINE__,B_NEG, 1,P_1Z,P_1Z1, 0,0,0, R_ALL, 0, 0,P_1Z,P_1Z1 },
 420   {__LINE__,B_NEG, 0,P_Z,P_Z1,  0,0,0, R_ALL, 0, 1,P_Z,P_Z1 },
 421   {__LINE__,B_NEG, 1,P_Z,P_Z1,  0,0,0, R_ALL, 0, 0,P_Z,P_Z1 },
 422
 423   /* Absolute value.  */
 424   {__LINE__,B_ABS, 0,P_Z,P_Z,   0,0,0, R_ALL, 0, 0,P_Z,P_Z },
 425   {__LINE__,B_ABS, 1,P_Z,P_Z,   0,0,0, R_ALL, 0, 0,P_Z,P_Z },
 426   {__LINE__,B_ABS, 0,P_O,P_Z,   0,0,0, R_ALL, 0, 0,P_O,P_Z },
 427   {__LINE__,B_ABS, 1,P_O,P_Z,   0,0,0, R_ALL, 0, 0,P_O,P_Z },
 428   {__LINE__,B_ABS, 0,P_O,P_1Z,  0,0,0, R_ALL, 0, 0,P_O,P_1Z },
 429   {__LINE__,B_ABS, 1,P_O,P_1Z,  0,0,0, R_ALL, 0, 0,P_O,P_1Z },
 430   {__LINE__,B_ABS, 0,P_O,P_01Z, 0,0,0, R_ALL, 0, 0,P_O,P_01Z },
 431   {__LINE__,B_ABS, 1,P_O,P_01Z, 0,0,0, R_ALL, 0, 0,P_O,P_01Z },
 432   {__LINE__,B_ABS, 0,P_1Z,P_1Z1, 0,0,0, R_ALL, 0, 0,P_1Z,P_1Z1 },
 433   {__LINE__,B_ABS, 1,P_1Z,P_1Z1, 0,0,0, R_ALL, 0, 0,P_1Z,P_1Z1 },
 434   {__LINE__,B_ABS, 0,P_Z,P_Z1,  0,0,0, R_ALL, 0, 0,P_Z,P_Z1 },
 435   {__LINE__,B_ABS, 1,P_Z,P_Z1,  0,0,0, R_ALL, 0, 0,P_Z,P_Z1 },
 436
 437   /* Square root.  */
 438   {__LINE__,B_SQRT, 0,P_Z,P_Z,   0,0,0, R_ALL, 0, 0,P_Z,P_Z },
 439   {__LINE__,B_SQRT, 1,P_Z,P_Z,   0,0,0, R_ALL, 0, 1,P_Z,P_Z },
 440   {__LINE__,B_SQRT, 0,P_O,P_1Z,  0,0,0, R_ALL, 0, 0,P_O,P_1Z },
 441   {__LINE__,B_SQRT, 1,P_O,P_1Z,  0,0,0, R_ALL, 0, 1,P_O,P_1Z },
 442   {__LINE__,B_SQRT, 0,P_O,P_01Z, 0,0,0, R_ALL,
 443    FE_INVALID | FE_INVALID_SNAN, 0,P_O,P_11Z },
 444   {__LINE__,B_SQRT, 1,P_O,P_01Z, 0,0,0, R_ALL,
 445    FE_INVALID | FE_INVALID_SNAN, 1,P_O,P_11Z },
 446
 447   {__LINE__,B_SQRT, 0,P_O,P_Z,   0,0,0, R_ALL, 0, 0,P_O,P_Z },
 448   {__LINE__,B_SQRT, 0,P_0O,P_Z,  0,0,0, R_ALL, 0, 0,P_0O,P_Z },
 449
 450   {__LINE__,B_SQRT, 1,P_O,P_Z,   0,0,0, R_ALL,
 451    FE_INVALID | FE_INVALID_SQRT, 0,P_O,P_1Z },
 452   {__LINE__,B_SQRT, 1,P_1Z,P_1Z1, 0,0,0, R_ALL,
 453    FE_INVALID | FE_INVALID_SQRT, 0,P_O,P_1Z },
 454   {__LINE__,B_SQRT, 1,P_Z,P_Z1,  0,0,0, R_ALL,
 455    FE_INVALID | FE_INVALID_SQRT, 0,P_O,P_1Z },
 456
 457 };
 458
 459 static void
 460 check_op(void)
 461 {
 462   size_t i;
 463   int j;
 464   tocheck_t r, a, b, x;
 465   int raised;
 466
 467   for (i = 0; i < sizeof(optests)/sizeof(optests[0]); i++)
 468     {
 469       a = pattern(optests[i].a_sgn, optests[i].a_exp,
 470                   optests[i].a_mant);
 471       b = pattern(optests[i].b_sgn, optests[i].b_exp,
 472                   optests[i].b_mant);
 473       x = pattern(optests[i].x_sgn, optests[i].x_exp,
 474                   optests[i].x_mant);
 475       for (j = 0; j < 4; j++)
 476         if (optests[i].rmode & 1<<j)
 477           {
 478             fesetenv(rmodes+j);
 479             switch (optests[i].op)
 480               {
 481               case B_ADD: r = a + b; break;
 482               case B_SUB: r = a - b; break;
 483               case B_MUL: r = a * b; break;
 484               case B_DIV: r = a / b; break;
 485               case B_NEG: r = -a; break;
 486               case B_ABS: r = FUNC(fabs)(a); break;
 487               case B_SQRT: r = FUNC(sqrt)(a); break;
 488               }
 489             raised = fetestexcept(all_exceptions);
 490             check_result(optests[i].line,rmnames[j],x,r);
 491             check_excepts(optests[i].line,rmnames[j],
 492                           optests[i].excepts,raised);
 493           }
 494     }
 495 }
 496
 497 static void
 498 fail_xr(int line, const char *rm, tocheck_t x, tocheck_t r, tocheck_t xx,
 499         int xflag)
 500 {
 501   size_t i;
 502   unsigned char *cx, *cr, *cxx;
 503
 504   printf("%s:%d:round %s:fail\n with x=0x", __FILE__, line,rm);
 505   cx = (unsigned char *)&x;
 506   cr = (unsigned char *)&r;
 507   cxx = (unsigned char *)&xx;
 508   for (i = 0; i < sizeof(tocheck_t); i++)
 509     printf("%02x", cx[i]);
 510   printf(" r=0x");
 511   for (i = 0; i < sizeof(tocheck_t); i++)
 512     printf("%02x", cr[i]);
 513   printf(" xx=0x");
 514   for (i = 0; i < sizeof(tocheck_t); i++)
 515     printf("%02x", cxx[i]);
 516   printf(" inexact=%d\n", xflag != 0);
 517   nerrors++;
 518 }
 519
 520 static void
 521 check_sqrt(tocheck_t a)
 522 {
 523   int j;
 524   tocheck_t r0, r1, r2, x0, x1, x2;
 525   int raised = 0;
 526   int ok;
 527
 528   for (j = 0; j < 4; j++)
 529     {
 530       int excepts;
 531
 532       fesetenv(rmodes+j);
 533       r1 = FUNC(sqrt)(a);
 534       excepts = fetestexcept(all_exceptions);
 535       fesetenv(FE_DFL_ENV);
 536       raised |= excepts & ~FE_INEXACT;
 537       x1 = r1 * r1 - a;
 538       if (excepts & FE_INEXACT)
 539         {
 540           r0 = delta(r1,-1); r2 = delta(r1,1);
 541           switch (1 << j)
 542             {
 543             case R_NEAREST:
 544               x0 = r0 * r0 - a; x2 = r2 * r2 - a;
 545               ok = fabs(x0) >= fabs(x1) && fabs(x1) <= fabs(x2);
 546               break;
 547             case R_ZERO:  case R_DOWN:
 548               x2 = r2 * r2 - a;
 549               ok = x1 <= 0 && x2 >= 0;
 550               break;
 551             case R_UP:
 552               x0 = r0 * r0 - a;
 553               ok = x1 >= 0 && x0 <= 0;
 554               break;
 555             default:
 556               assert(0);
 557             }
 558         }
 559       else
 560         ok = x1 == 0;
 561       if (!ok)
 562         fail_xr(__LINE__,rmnames[j],a,r1,x1,excepts&FE_INEXACT);
 563     }
 564   check_excepts(__LINE__,"all",0,raised);
 565 }
 566
 567 int main(int argc, char **argv)
 568 {
 569   int i;
 570
 571   _LIB_VERSION = _IEEE_;
 572
 573   /* Set up environments for rounding modes.  */
 574   fesetenv(FE_DFL_ENV);
 575   fesetround(FE_TONEAREST);
 576   fegetenv(rmodes+0);
 577   fesetround(FE_TOWARDZERO);
 578   fegetenv(rmodes+1);
 579   fesetround(FE_UPWARD);
 580   fegetenv(rmodes+2);
 581   fesetround(FE_DOWNWARD);
 582   fegetenv(rmodes+3);
 583
 584 #if defined(FE_INVALID_SOFTWARE) || defined(FE_INVALID_SQRT)
 585   /* There's this really stupid feature of the 601... */
 586   fesetenv(FE_DFL_ENV);
 587   feraiseexcept(FE_INVALID_SOFTWARE);
 588   if (!fetestexcept(FE_INVALID_SOFTWARE))
 589     excepts_missing |= FE_INVALID_SOFTWARE;
 590   fesetenv(FE_DFL_ENV);
 591   feraiseexcept(FE_INVALID_SQRT);
 592   if (!fetestexcept(FE_INVALID_SQRT))
 593     excepts_missing |= FE_INVALID_SQRT;
 594 #endif
 595
 596   check_op();
 597   for (i = 0; i < 100000; i++)
 598     check_sqrt(pattern(0, P_Rno, P_R));
 599   for (i = 0; i < 100; i++)
 600     check_sqrt(pattern(0, P_Z, P_R));
 601   check_sqrt(pattern(0,P_Z,P_Z1));
 602
 603   printf("%d errors.\n", nerrors);
 604   return nerrors == 0 ? 0 : 1;
 605 }