Back to index

glibc  2.9
test-arith.c
Go to the documentation of this file.
00001 /* Test floating-point arithmetic operations.
00002    Copyright (C) 1997, 1998 Free Software Foundation, Inc.
00003    This file is part of the GNU C Library.
00004 
00005    The GNU C Library is free software; you can redistribute it and/or
00006    modify it under the terms of the GNU Lesser General Public
00007    License as published by the Free Software Foundation; either
00008    version 2.1 of the License, or (at your option) any later version.
00009 
00010    The GNU C Library is distributed in the hope that it will be useful,
00011    but WITHOUT ANY WARRANTY; without even the implied warranty of
00012    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00013    Lesser General Public License for more details.
00014 
00015    You should have received a copy of the GNU Lesser General Public
00016    License along with the GNU C Library; if not, write to the Free
00017    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
00018    02111-1307 USA.  */
00019 #ifndef _GNU_SOURCE
00020 #define _GNU_SOURCE
00021 #endif
00022 #include <math.h>
00023 #include <stdio.h>
00024 #include <stdlib.h>
00025 #include <string.h>
00026 #include <fenv.h>
00027 #include <assert.h>
00028 
00029 #ifndef ESIZE
00030 typedef double tocheck_t;
00031 #define ESIZE 11
00032 #define MSIZE 52
00033 #define FUNC(x) x
00034 #endif
00035 
00036 #define R_NEAREST 1
00037 #define R_ZERO 2
00038 #define R_UP 4
00039 #define R_DOWN 8
00040 #define R_ALL (R_NEAREST|R_ZERO|R_UP|R_DOWN)
00041 static fenv_t rmodes[4];
00042 static const char * const rmnames[4] =
00043 { "nearest","zero","+Inf","-Inf" };
00044 
00045 typedef union {
00046   tocheck_t tc;
00047   unsigned char c[sizeof(tocheck_t)];
00048 } union_t;
00049 
00050 /* Don't try reading these in a font that doesn't distinguish
00051    O and zero.  */
00052 typedef enum {
00053   P_Z    = 0x0,  /* 00000...0 */
00054   P_000O = 0x1,  /* 00011...1 */
00055   P_001Z = 0x2,  /* 00100...0 */
00056   P_00O  = 0x3,  /* 00111...1 */
00057   P_01Z  = 0x4,  /* 01000...0 */
00058   P_010O = 0x5,  /* 01011...1 */
00059   P_011Z = 0x6,  /* 01100...0 */
00060   P_0O   = 0x7,  /* 01111...1 */
00061   P_1Z   = 0x8,  /* 10000...0 */
00062   P_100O = 0x9,  /* 10011...1 */
00063   P_101Z = 0xa,  /* 10100...0 */
00064   P_10O  = 0xb,  /* 10111...1 */
00065   P_11Z  = 0xc,  /* 11000...0 */
00066   P_110O = 0xd,  /* 11011...1 */
00067   P_111Z = 0xe,  /* 11100...0 */
00068   P_O    = 0xf,  /* 11111...1 */
00069   P_Z1   = 0x11, /* 000...001 */
00070   P_Z10  = 0x12, /* 000...010 */
00071   P_Z11  = 0x13, /* 000...011 */
00072   P_0O00 = 0x14, /* 011...100 */
00073   P_0O01 = 0x15, /* 011...101 */
00074   P_0O0  = 0x16, /* 011...110 */
00075   P_1Z1  = 0x19, /* 100...001 */
00076   P_1Z10 = 0x1a, /* 100...010 */
00077   P_1Z11 = 0x1b, /* 100...011 */
00078   P_O00  = 0x1c, /* 111...100 */
00079   P_O01  = 0x1d, /* 111...101 */
00080   P_O0   = 0x1e, /* 111...110 */
00081   P_R    = 0x20, /* rrr...rrr */ /* ('r' means random. ) */
00082   P_Ro   = 0x21, /* rrr...rrr, with odd parity.  */
00083   P_0R   = 0x22, /* 0rr...rrr */
00084   P_1R   = 0x23, /* 1rr...rrr */
00085   P_Rno  = 0x24, /* rrr...rrr, but not all ones.  */
00086 } pattern_t;
00087 
00088 static void
00089 pattern_fill(pattern_t ptn, unsigned char *start, int bitoffset, int count)
00090 {
00091 #define bitset(count, value) \
00092       start[(count)/8] = (start[(count)/8] & ~(1 << 7-(count)%8)  \
00093                           |  (value) << 7-(count)%8)
00094   int i;
00095 
00096   if (ptn >= 0 && ptn <= 0xf)
00097     {
00098       /* Patterns between 0 and 0xF have the following format:
00099         The LSBit is used to fill the last n-3 bits of the pattern;
00100         The next 3 bits are the first 3 bits of the pattern. */
00101       for (i = 0; i < count; i++)
00102        if (i < 3)
00103          bitset((bitoffset+i), ptn >> (3-i) & 1);
00104        else
00105          bitset((bitoffset+i), ptn >> 0 & 1);
00106     }
00107   else if (ptn <= 0x1f)
00108     {
00109       /* Patterns between 0x10 and 0x1F have the following format:
00110         The two LSBits are the last two bits of the pattern;
00111         The 0x8 bit is the first bit of the pattern;
00112         The 0x4 bit is used to fill the remainder. */
00113       for (i = 0; i < count; i++)
00114        if (i == 0)
00115          bitset((bitoffset+i), ptn >> 3 & 1);
00116        else if (i >= count-2)
00117          bitset((bitoffset+i), ptn >> (count-1-i) & 1);
00118        else
00119          bitset((bitoffset+i), ptn >> 2 & 1);
00120     }
00121   else switch (ptn)
00122     {
00123     case P_0R: case P_1R:
00124       assert(count > 0);
00125       bitset(bitoffset, ptn & 1);
00126       count--;
00127       bitoffset++;
00128     case P_R:
00129       for (; count > 0; count--, bitoffset++)
00130        bitset(bitoffset, rand() & 1);
00131       break;
00132     case P_Ro:
00133       {
00134        int op = 1;
00135        assert(count > 0);
00136        for (; count > 1; count--, bitoffset++)
00137          bitset(bitoffset, op ^= (rand() & 1));
00138        bitset(bitoffset, op);
00139        break;
00140       }
00141     case P_Rno:
00142       {
00143        int op = 1;
00144        assert(count > 0);
00145        for (; count > 1; count--, bitoffset++)
00146        {
00147          int r = rand() & 1;
00148          op &= r;
00149          bitset(bitoffset, r);
00150        }
00151        bitset(bitoffset, rand() & (op ^ 1));
00152        break;
00153       }
00154 
00155     default:
00156       assert(0);
00157     }
00158 #undef bitset
00159 }
00160 
00161 static tocheck_t
00162 pattern(int negative, pattern_t exp, pattern_t mant)
00163 {
00164   union_t result;
00165 #if 0
00166   int i;
00167 #endif
00168 
00169   pattern_fill(negative ? P_O : P_Z, result.c, 0, 1);
00170   pattern_fill(exp, result.c, 1, ESIZE);
00171   pattern_fill(mant, result.c, ESIZE+1, MSIZE);
00172 #if 0
00173   printf("neg=%d exp=%02x mant=%02x: ", negative, exp, mant);
00174   for (i = 0; i < sizeof(tocheck_t); i++)
00175     printf("%02x", result.c[i]);
00176   printf("\n");
00177 #endif
00178   return result.tc;
00179 }
00180 
00181 /* Return the closest different tocheck_t to 'x' in the direction of
00182    'direction', or 'x' if there is no such value.  Assumes 'x' is not
00183    a NaN.  */
00184 static tocheck_t
00185 delta(tocheck_t x, int direction)
00186 {
00187   union_t xx;
00188   int i;
00189 
00190   xx.tc = x;
00191   if (xx.c[0] & 0x80)
00192     direction = -direction;
00193   if (direction == +1)
00194     {
00195       union_t tx;
00196       tx.tc = pattern(xx.c[0] >> 7, P_O, P_Z);
00197       if (memcmp(tx.c, xx.c, sizeof(tocheck_t)) == 0)
00198        return x;
00199     }
00200   for (i = sizeof(tocheck_t)-1; i > 0; i--)
00201     {
00202       xx.c[i] += direction;
00203       if (xx.c[i] != (direction > 0 ? 0 : 0xff))
00204        return xx.tc;
00205     }
00206   if (direction < 0 && (xx.c[0] & 0x7f) == 0)
00207     return pattern(~(xx.c[0] >> 7) & 1, P_Z, P_Z1);
00208   else
00209     {
00210       xx.c[0] += direction;
00211       return xx.tc;
00212     }
00213 }
00214 
00215 static int nerrors = 0;
00216 
00217 #ifdef FE_ALL_INVALID
00218 static const int all_exceptions = FE_ALL_INVALID | FE_ALL_EXCEPT;
00219 #else
00220 static const int all_exceptions = FE_ALL_EXCEPT;
00221 #endif
00222 
00223 static void
00224 check_result(int line, const char *rm, tocheck_t expected, tocheck_t actual)
00225 {
00226   if (memcmp(&expected, &actual, sizeof(tocheck_t)) != 0)
00227     {
00228       unsigned char *ex, *ac;
00229       size_t i;
00230 
00231       printf("%s:%d:round %s:result failed\n"
00232             " expected result 0x", __FILE__, line, rm);
00233       ex = (unsigned char *)&expected;
00234       ac = (unsigned char *)&actual;
00235       for (i = 0; i < sizeof(tocheck_t); i++)
00236        printf("%02x", ex[i]);
00237       printf(" got 0x");
00238       for (i = 0; i < sizeof(tocheck_t); i++)
00239        printf("%02x", ac[i]);
00240       printf("\n");
00241       nerrors++;
00242     }
00243 }
00244 
00245 static const struct {
00246   int except;
00247   const char *name;
00248 } excepts[] = {
00249 #define except_entry(ex) { ex, #ex } ,
00250 #ifdef FE_INEXACT
00251   except_entry(FE_INEXACT)
00252 #else
00253 # define FE_INEXACT 0
00254 #endif
00255 #ifdef FE_DIVBYZERO
00256   except_entry(FE_DIVBYZERO)
00257 #else
00258 # define FE_DIVBYZERO 0
00259 #endif
00260 #ifdef FE_UNDERFLOW
00261   except_entry(FE_UNDERFLOW)
00262 #else
00263 # define FE_UNDERFLOW 0
00264 #endif
00265 #ifdef FE_OVERFLOW
00266   except_entry(FE_OVERFLOW)
00267 #else
00268 # define FE_OVERFLOW 0
00269 #endif
00270 #ifdef FE_INVALID
00271   except_entry(FE_INVALID)
00272 #else
00273 # define FE_INVALID 0
00274 #endif
00275 #ifdef FE_INVALID_SNAN
00276   except_entry(FE_INVALID_SNAN)
00277 #else
00278 # define FE_INVALID_SNAN FE_INVALID
00279 #endif
00280 #ifdef FE_INVALID_ISI
00281   except_entry(FE_INVALID_ISI)
00282 #else
00283 # define FE_INVALID_ISI FE_INVALID
00284 #endif
00285 #ifdef FE_INVALID_IDI
00286   except_entry(FE_INVALID_IDI)
00287 #else
00288 # define FE_INVALID_IDI FE_INVALID
00289 #endif
00290 #ifdef FE_INVALID_ZDZ
00291   except_entry(FE_INVALID_ZDZ)
00292 #else
00293 # define FE_INVALID_ZDZ FE_INVALID
00294 #endif
00295 #ifdef FE_INVALID_COMPARE
00296   except_entry(FE_INVALID_COMPARE)
00297 #else
00298 # define FE_INVALID_COMPARE FE_INVALID
00299 #endif
00300 #ifdef FE_INVALID_SOFTWARE
00301   except_entry(FE_INVALID_SOFTWARE)
00302 #else
00303 # define FE_INVALID_SOFTWARE FE_INVALID
00304 #endif
00305 #ifdef FE_INVALID_SQRT
00306   except_entry(FE_INVALID_SQRT)
00307 #else
00308 # define FE_INVALID_SQRT FE_INVALID
00309 #endif
00310 #ifdef FE_INVALID_INTEGER_CONVERSION
00311   except_entry(FE_INVALID_INTEGER_CONVERSION)
00312 #else
00313 # define FE_INVALID_INTEGER_CONVERSION FE_INVALID
00314 #endif
00315 };
00316 
00317 static int excepts_missing = 0;
00318 
00319 static void
00320 check_excepts(int line, const char *rm, int expected, int actual)
00321 {
00322   if (expected & excepts_missing)
00323     expected = expected & ~excepts_missing | FE_INVALID_SNAN;
00324   if ((expected & all_exceptions) != actual)
00325     {
00326       size_t i;
00327       printf("%s:%d:round %s:exceptions failed\n"
00328             " expected exceptions ", __FILE__, line,rm);
00329       for (i = 0; i < sizeof(excepts)/sizeof(excepts[0]); i++)
00330        if (expected & excepts[i].except)
00331          printf("%s ",excepts[i].name);
00332       if ((expected & all_exceptions) == 0)
00333        printf("- ");
00334       printf("got");
00335       for (i = 0; i < sizeof(excepts)/sizeof(excepts[0]); i++)
00336        if (actual & excepts[i].except)
00337          printf(" %s",excepts[i].name);
00338       if ((actual & all_exceptions) == 0)
00339        printf("- ");
00340       printf(".\n");
00341       nerrors++;
00342     }
00343 }
00344 
00345 typedef enum {
00346   B_ADD, B_SUB, B_MUL, B_DIV, B_NEG, B_ABS, B_SQRT
00347 } op_t;
00348 typedef struct {
00349   int line;
00350   op_t op;
00351   int a_sgn;
00352   pattern_t a_exp, a_mant;
00353   int b_sgn;
00354   pattern_t b_exp, b_mant;
00355   int rmode;
00356   int excepts;
00357   int x_sgn;
00358   pattern_t x_exp, x_mant;
00359 } optest_t;
00360 static const optest_t optests[] = {
00361   /* Additions of zero.  */
00362   {__LINE__,B_ADD, 0,P_Z,P_Z, 0,P_Z,P_Z, R_ALL,0, 0,P_Z,P_Z },
00363   {__LINE__,B_ADD, 1,P_Z,P_Z, 0,P_Z,P_Z, R_ALL & ~R_DOWN,0, 0,P_Z,P_Z },
00364   {__LINE__,B_ADD, 1,P_Z,P_Z, 0,P_Z,P_Z, R_DOWN,0, 1,P_Z,P_Z },
00365   {__LINE__,B_ADD, 1,P_Z,P_Z, 1,P_Z,P_Z, R_ALL,0, 1,P_Z,P_Z },
00366 
00367   /* Additions with NaN.  */
00368   {__LINE__,B_ADD, 0,P_O,P_101Z, 0,P_Z,P_Z, R_ALL,0, 0,P_O,P_101Z },
00369   {__LINE__,B_ADD, 0,P_O,P_01Z, 0,P_Z,P_Z, R_ALL,
00370    FE_INVALID | FE_INVALID_SNAN, 0,P_O,P_11Z },
00371   {__LINE__,B_ADD, 0,P_O,P_Z, 0,P_O,P_0O, R_ALL,
00372    FE_INVALID | FE_INVALID_SNAN, 0,P_O,P_O },
00373   {__LINE__,B_ADD, 0,P_Z,P_Z, 0,P_O,P_11Z, R_ALL,0, 0,P_O,P_11Z },
00374   {__LINE__,B_ADD, 0,P_O,P_001Z, 0,P_O,P_001Z, R_ALL,
00375    FE_INVALID | FE_INVALID_SNAN, 0,P_O,P_101Z },
00376   {__LINE__,B_ADD, 0,P_O,P_1Z, 0,P_Z,P_Z, R_ALL,0, 0,P_O,P_1Z },
00377   {__LINE__,B_ADD, 0,P_0O,P_Z, 0,P_O,P_10O, R_ALL,0, 0,P_O,P_10O },
00378 
00379   /* Additions with infinity.  */
00380   {__LINE__,B_ADD, 0,P_O,P_Z, 0,P_Z,P_Z, R_ALL,0, 0,P_O,P_Z },
00381   {__LINE__,B_ADD, 0,P_O,P_Z, 1,P_Z,P_Z, R_ALL,0, 0,P_O,P_Z },
00382   {__LINE__,B_ADD, 1,P_O,P_Z, 0,P_Z,P_Z, R_ALL,0, 1,P_O,P_Z },
00383   {__LINE__,B_ADD, 1,P_O,P_Z, 1,P_Z,P_Z, R_ALL,0, 1,P_O,P_Z },
00384   {__LINE__,B_ADD, 0,P_O,P_Z, 0,P_O,P_Z, R_ALL,0, 0,P_O,P_Z },
00385   {__LINE__,B_ADD, 1,P_O,P_Z, 1,P_O,P_Z, R_ALL,0, 1,P_O,P_Z },
00386   {__LINE__,B_ADD, 0,P_O,P_Z, 1,P_O,P_Z, R_ALL,
00387    FE_INVALID | FE_INVALID_ISI, 0,P_O,P_1Z },
00388   {__LINE__,B_ADD, 1,P_O,P_Z, 0,P_O,P_Z, R_ALL,
00389    FE_INVALID | FE_INVALID_ISI, 0,P_O,P_1Z },
00390   {__LINE__,B_ADD, 0,P_O,P_Z, 0,P_0O,P_Z, R_ALL,0, 0,P_O,P_Z },
00391   {__LINE__,B_ADD, 1,P_O,P_Z, 0,P_0O,P_Z, R_ALL,0, 1,P_O,P_Z },
00392   {__LINE__,B_ADD, 0,P_O,P_Z, 1,P_0O,P_Z, R_ALL,0, 0,P_O,P_Z },
00393   {__LINE__,B_ADD, 1,P_O,P_Z, 1,P_0O,P_Z, R_ALL,0, 1,P_O,P_Z },
00394 
00395   /* Overflow (and zero).  */
00396   {__LINE__,B_ADD, 0,P_O0,P_Z, 0,P_O0,P_Z, R_NEAREST | R_UP,
00397    FE_INEXACT | FE_OVERFLOW, 0,P_O,P_Z },
00398   {__LINE__,B_ADD, 0,P_O0,P_Z, 0,P_O0,P_Z, R_ZERO | R_DOWN,
00399    FE_INEXACT | FE_OVERFLOW, 0,P_O0,P_O },
00400   {__LINE__,B_ADD, 1,P_O0,P_Z, 1,P_O0,P_Z, R_NEAREST | R_DOWN,
00401    FE_INEXACT | FE_OVERFLOW, 1,P_O,P_Z },
00402   {__LINE__,B_ADD, 1,P_O0,P_Z, 1,P_O0,P_Z, R_ZERO | R_UP,
00403    FE_INEXACT | FE_OVERFLOW, 1,P_O0,P_O },
00404   {__LINE__,B_ADD, 0,P_O0,P_Z, 1,P_O0,P_Z, R_ALL & ~R_DOWN,
00405    0, 0,P_Z,P_Z },
00406   {__LINE__,B_ADD, 0,P_O0,P_Z, 1,P_O0,P_Z, R_DOWN,
00407    0, 1,P_Z,P_Z },
00408 
00409   /* Negation.  */
00410   {__LINE__,B_NEG, 0,P_Z,P_Z,   0,0,0, R_ALL, 0, 1,P_Z,P_Z },
00411   {__LINE__,B_NEG, 1,P_Z,P_Z,   0,0,0, R_ALL, 0, 0,P_Z,P_Z },
00412   {__LINE__,B_NEG, 0,P_O,P_Z,   0,0,0, R_ALL, 0, 1,P_O,P_Z },
00413   {__LINE__,B_NEG, 1,P_O,P_Z,   0,0,0, R_ALL, 0, 0,P_O,P_Z },
00414   {__LINE__,B_NEG, 0,P_O,P_1Z,  0,0,0, R_ALL, 0, 1,P_O,P_1Z },
00415   {__LINE__,B_NEG, 1,P_O,P_1Z,  0,0,0, R_ALL, 0, 0,P_O,P_1Z },
00416   {__LINE__,B_NEG, 0,P_O,P_01Z, 0,0,0, R_ALL, 0, 1,P_O,P_01Z },
00417   {__LINE__,B_NEG, 1,P_O,P_01Z, 0,0,0, R_ALL, 0, 0,P_O,P_01Z },
00418   {__LINE__,B_NEG, 0,P_1Z,P_1Z1, 0,0,0, R_ALL, 0, 1,P_1Z,P_1Z1 },
00419   {__LINE__,B_NEG, 1,P_1Z,P_1Z1, 0,0,0, R_ALL, 0, 0,P_1Z,P_1Z1 },
00420   {__LINE__,B_NEG, 0,P_Z,P_Z1,  0,0,0, R_ALL, 0, 1,P_Z,P_Z1 },
00421   {__LINE__,B_NEG, 1,P_Z,P_Z1,  0,0,0, R_ALL, 0, 0,P_Z,P_Z1 },
00422 
00423   /* Absolute value.  */
00424   {__LINE__,B_ABS, 0,P_Z,P_Z,   0,0,0, R_ALL, 0, 0,P_Z,P_Z },
00425   {__LINE__,B_ABS, 1,P_Z,P_Z,   0,0,0, R_ALL, 0, 0,P_Z,P_Z },
00426   {__LINE__,B_ABS, 0,P_O,P_Z,   0,0,0, R_ALL, 0, 0,P_O,P_Z },
00427   {__LINE__,B_ABS, 1,P_O,P_Z,   0,0,0, R_ALL, 0, 0,P_O,P_Z },
00428   {__LINE__,B_ABS, 0,P_O,P_1Z,  0,0,0, R_ALL, 0, 0,P_O,P_1Z },
00429   {__LINE__,B_ABS, 1,P_O,P_1Z,  0,0,0, R_ALL, 0, 0,P_O,P_1Z },
00430   {__LINE__,B_ABS, 0,P_O,P_01Z, 0,0,0, R_ALL, 0, 0,P_O,P_01Z },
00431   {__LINE__,B_ABS, 1,P_O,P_01Z, 0,0,0, R_ALL, 0, 0,P_O,P_01Z },
00432   {__LINE__,B_ABS, 0,P_1Z,P_1Z1, 0,0,0, R_ALL, 0, 0,P_1Z,P_1Z1 },
00433   {__LINE__,B_ABS, 1,P_1Z,P_1Z1, 0,0,0, R_ALL, 0, 0,P_1Z,P_1Z1 },
00434   {__LINE__,B_ABS, 0,P_Z,P_Z1,  0,0,0, R_ALL, 0, 0,P_Z,P_Z1 },
00435   {__LINE__,B_ABS, 1,P_Z,P_Z1,  0,0,0, R_ALL, 0, 0,P_Z,P_Z1 },
00436 
00437   /* Square root.  */
00438   {__LINE__,B_SQRT, 0,P_Z,P_Z,   0,0,0, R_ALL, 0, 0,P_Z,P_Z },
00439   {__LINE__,B_SQRT, 1,P_Z,P_Z,   0,0,0, R_ALL, 0, 1,P_Z,P_Z },
00440   {__LINE__,B_SQRT, 0,P_O,P_1Z,  0,0,0, R_ALL, 0, 0,P_O,P_1Z },
00441   {__LINE__,B_SQRT, 1,P_O,P_1Z,  0,0,0, R_ALL, 0, 1,P_O,P_1Z },
00442   {__LINE__,B_SQRT, 0,P_O,P_01Z, 0,0,0, R_ALL,
00443    FE_INVALID | FE_INVALID_SNAN, 0,P_O,P_11Z },
00444   {__LINE__,B_SQRT, 1,P_O,P_01Z, 0,0,0, R_ALL,
00445    FE_INVALID | FE_INVALID_SNAN, 1,P_O,P_11Z },
00446 
00447   {__LINE__,B_SQRT, 0,P_O,P_Z,   0,0,0, R_ALL, 0, 0,P_O,P_Z },
00448   {__LINE__,B_SQRT, 0,P_0O,P_Z,  0,0,0, R_ALL, 0, 0,P_0O,P_Z },
00449 
00450   {__LINE__,B_SQRT, 1,P_O,P_Z,   0,0,0, R_ALL,
00451    FE_INVALID | FE_INVALID_SQRT, 0,P_O,P_1Z },
00452   {__LINE__,B_SQRT, 1,P_1Z,P_1Z1, 0,0,0, R_ALL,
00453    FE_INVALID | FE_INVALID_SQRT, 0,P_O,P_1Z },
00454   {__LINE__,B_SQRT, 1,P_Z,P_Z1,  0,0,0, R_ALL,
00455    FE_INVALID | FE_INVALID_SQRT, 0,P_O,P_1Z },
00456 
00457 };
00458 
00459 static void
00460 check_op(void)
00461 {
00462   size_t i;
00463   int j;
00464   tocheck_t r, a, b, x;
00465   int raised;
00466 
00467   for (i = 0; i < sizeof(optests)/sizeof(optests[0]); i++)
00468     {
00469       a = pattern(optests[i].a_sgn, optests[i].a_exp,
00470                 optests[i].a_mant);
00471       b = pattern(optests[i].b_sgn, optests[i].b_exp,
00472                 optests[i].b_mant);
00473       x = pattern(optests[i].x_sgn, optests[i].x_exp,
00474                 optests[i].x_mant);
00475       for (j = 0; j < 4; j++)
00476        if (optests[i].rmode & 1<<j)
00477          {
00478            fesetenv(rmodes+j);
00479            switch (optests[i].op)
00480              {
00481              case B_ADD: r = a + b; break;
00482              case B_SUB: r = a - b; break;
00483              case B_MUL: r = a * b; break;
00484              case B_DIV: r = a / b; break;
00485              case B_NEG: r = -a; break;
00486              case B_ABS: r = FUNC(fabs)(a); break;
00487              case B_SQRT: r = FUNC(sqrt)(a); break;
00488              }
00489            raised = fetestexcept(all_exceptions);
00490            check_result(optests[i].line,rmnames[j],x,r);
00491            check_excepts(optests[i].line,rmnames[j],
00492                        optests[i].excepts,raised);
00493          }
00494     }
00495 }
00496 
00497 static void
00498 fail_xr(int line, const char *rm, tocheck_t x, tocheck_t r, tocheck_t xx,
00499        int xflag)
00500 {
00501   size_t i;
00502   unsigned char *cx, *cr, *cxx;
00503 
00504   printf("%s:%d:round %s:fail\n with x=0x", __FILE__, line,rm);
00505   cx = (unsigned char *)&x;
00506   cr = (unsigned char *)&r;
00507   cxx = (unsigned char *)&xx;
00508   for (i = 0; i < sizeof(tocheck_t); i++)
00509     printf("%02x", cx[i]);
00510   printf(" r=0x");
00511   for (i = 0; i < sizeof(tocheck_t); i++)
00512     printf("%02x", cr[i]);
00513   printf(" xx=0x");
00514   for (i = 0; i < sizeof(tocheck_t); i++)
00515     printf("%02x", cxx[i]);
00516   printf(" inexact=%d\n", xflag != 0);
00517   nerrors++;
00518 }
00519 
00520 static void
00521 check_sqrt(tocheck_t a)
00522 {
00523   int j;
00524   tocheck_t r0, r1, r2, x0, x1, x2;
00525   int raised = 0;
00526   int ok;
00527 
00528   for (j = 0; j < 4; j++)
00529     {
00530       int excepts;
00531 
00532       fesetenv(rmodes+j);
00533       r1 = FUNC(sqrt)(a);
00534       excepts = fetestexcept(all_exceptions);
00535       fesetenv(FE_DFL_ENV);
00536       raised |= excepts & ~FE_INEXACT;
00537       x1 = r1 * r1 - a;
00538       if (excepts & FE_INEXACT)
00539        {
00540          r0 = delta(r1,-1); r2 = delta(r1,1);
00541          switch (1 << j)
00542            {
00543            case R_NEAREST:
00544              x0 = r0 * r0 - a; x2 = r2 * r2 - a;
00545              ok = fabs(x0) >= fabs(x1) && fabs(x1) <= fabs(x2);
00546              break;
00547            case R_ZERO:  case R_DOWN:
00548              x2 = r2 * r2 - a;
00549              ok = x1 <= 0 && x2 >= 0;
00550              break;
00551            case R_UP:
00552              x0 = r0 * r0 - a;
00553              ok = x1 >= 0 && x0 <= 0;
00554              break;
00555            default:
00556              assert(0);
00557            }
00558        }
00559       else
00560        ok = x1 == 0;
00561       if (!ok)
00562        fail_xr(__LINE__,rmnames[j],a,r1,x1,excepts&FE_INEXACT);
00563     }
00564   check_excepts(__LINE__,"all",0,raised);
00565 }
00566 
00567 int main(int argc, char **argv)
00568 {
00569   int i;
00570 
00571   _LIB_VERSION = _IEEE_;
00572 
00573   /* Set up environments for rounding modes.  */
00574   fesetenv(FE_DFL_ENV);
00575   fesetround(FE_TONEAREST);
00576   fegetenv(rmodes+0);
00577   fesetround(FE_TOWARDZERO);
00578   fegetenv(rmodes+1);
00579   fesetround(FE_UPWARD);
00580   fegetenv(rmodes+2);
00581   fesetround(FE_DOWNWARD);
00582   fegetenv(rmodes+3);
00583 
00584 #if defined(FE_INVALID_SOFTWARE) || defined(FE_INVALID_SQRT)
00585   /* There's this really stupid feature of the 601... */
00586   fesetenv(FE_DFL_ENV);
00587   feraiseexcept(FE_INVALID_SOFTWARE);
00588   if (!fetestexcept(FE_INVALID_SOFTWARE))
00589     excepts_missing |= FE_INVALID_SOFTWARE;
00590   fesetenv(FE_DFL_ENV);
00591   feraiseexcept(FE_INVALID_SQRT);
00592   if (!fetestexcept(FE_INVALID_SQRT))
00593     excepts_missing |= FE_INVALID_SQRT;
00594 #endif
00595 
00596   check_op();
00597   for (i = 0; i < 100000; i++)
00598     check_sqrt(pattern(0, P_Rno, P_R));
00599   for (i = 0; i < 100; i++)
00600     check_sqrt(pattern(0, P_Z, P_R));
00601   check_sqrt(pattern(0,P_Z,P_Z1));
00602 
00603   printf("%d errors.\n", nerrors);
00604   return nerrors == 0 ? 0 : 1;
00605 }