tools/regression/lib/msun/test-fma.c

   1 /*-
   2  * Copyright (c) 2008 David Schultz <das@FreeBSD.org>
   3  * All rights reserved.
   4  *
   5  * Redistribution and use in source and binary forms, with or without
   6  * modification, are permitted provided that the following conditions
   7  * are met:
   8  * 1. Redistributions of source code must retain the above copyright
   9  *    notice, this list of conditions and the following disclaimer.
  10  * 2. Redistributions in binary form must reproduce the above copyright
  11  *    notice, this list of conditions and the following disclaimer in the
  12  *    documentation and/or other materials provided with the distribution.
  13  *
  14  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
  15  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  16  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  17  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  18  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  19  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  20  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  21  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  22  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  23  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  24  * SUCH DAMAGE.
  25  */
  26
  27 /*
  28  * Tests for fma{,f,l}().
  29  */
  30
  31 #include <sys/cdefs.h>
  32 __FBSDID("$FreeBSD$");
  33
  34 #include <sys/param.h>
  35 #include <assert.h>
  36 #include <fenv.h>
  37 #include <float.h>
  38 #include <math.h>
  39 #include <stdio.h>
  40 #include <stdlib.h>
  41
  42 #define ALL_STD_EXCEPT  (FE_DIVBYZERO | FE_INEXACT | FE_INVALID | \
  43                          FE_OVERFLOW | FE_UNDERFLOW)
  44
  45 #pragma STDC FENV_ACCESS ON
  46
  47 /*
  48  * Test that a function returns the correct value and sets the
  49  * exception flags correctly. The exceptmask specifies which
  50  * exceptions we should check. We need to be lenient for several
  51  * reasons, but mainly because on some architectures it's impossible
  52  * to raise FE_OVERFLOW without raising FE_INEXACT.
  53  *
  54  * These are macros instead of functions so that assert provides more
  55  * meaningful error messages.
  56  */
  57 #define test(func, x, y, z, result, exceptmask, excepts) do {           \
  58         assert(feclearexcept(FE_ALL_EXCEPT) == 0);                      \
  59         assert(fpequal((func)((x), (y), (z)), (result)));               \
  60         assert(((func), fetestexcept(exceptmask) == (excepts)));        \
  61 } while (0)
  62
  63 #define testall(x, y, z, result, exceptmask, excepts)   do {            \
  64         test(fma, (x), (y), (z), (double)(result), (exceptmask), (excepts)); \
  65         test(fmaf, (x), (y), (z), (float)(result), (exceptmask), (excepts)); \
  66         test(fmal, (x), (y), (z), (result), (exceptmask), (excepts));   \
  67 } while (0)
  68
  69 /* Test in all rounding modes. */
  70 #define testrnd(func, x, y, z, rn, ru, rd, rz, exceptmask, excepts)     do { \
  71         fesetround(FE_TONEAREST);                                       \
  72         test((func), (x), (y), (z), (rn), (exceptmask), (excepts));     \
  73         fesetround(FE_UPWARD);                                          \
  74         test((func), (x), (y), (z), (ru), (exceptmask), (excepts));     \
  75         fesetround(FE_DOWNWARD);                                        \
  76         test((func), (x), (y), (z), (rd), (exceptmask), (excepts));     \
  77         fesetround(FE_TOWARDZERO);                                      \
  78         test((func), (x), (y), (z), (rz), (exceptmask), (excepts));     \
  79 } while (0)
  80
  81 /*
  82  * Determine whether x and y are equal, with two special rules:
  83  *      +0.0 != -0.0
  84  *       NaN == NaN
  85  */
  86 int
  87 fpequal(long double x, long double y)
  88 {
  89
  90         return ((x == y && !signbit(x) == !signbit(y))
  91                 || (isnan(x) && isnan(y)));
  92 }
  93
  94 static void
  95 test_zeroes(void)
  96 {
  97         const int rd = (fegetround() == FE_DOWNWARD);
  98
  99         testall(0.0, 0.0, 0.0, 0.0, ALL_STD_EXCEPT, 0);
 100         testall(1.0, 0.0, 0.0, 0.0, ALL_STD_EXCEPT, 0);
 101         testall(0.0, 1.0, 0.0, 0.0, ALL_STD_EXCEPT, 0);
 102         testall(0.0, 0.0, 1.0, 1.0, ALL_STD_EXCEPT, 0);
 103
 104         testall(-0.0, 0.0, 0.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0);
 105         testall(0.0, -0.0, 0.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0);
 106         testall(-0.0, -0.0, 0.0, 0.0, ALL_STD_EXCEPT, 0);
 107         testall(0.0, 0.0, -0.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0);
 108         testall(-0.0, -0.0, -0.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0);
 109
 110         testall(-0.0, 0.0, -0.0, -0.0, ALL_STD_EXCEPT, 0);
 111         testall(0.0, -0.0, -0.0, -0.0, ALL_STD_EXCEPT, 0);
 112
 113         testall(-1.0, 1.0, 1.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0);
 114         testall(1.0, -1.0, 1.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0);
 115         testall(-1.0, -1.0, -1.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0);
 116
 117         switch (fegetround()) {
 118         case FE_TONEAREST:
 119         case FE_TOWARDZERO:
 120                 test(fmaf, -FLT_MIN, FLT_MIN, 0.0, -0.0,
 121                      ALL_STD_EXCEPT, FE_INEXACT | FE_UNDERFLOW);
 122                 test(fma, -DBL_MIN, DBL_MIN, 0.0, -0.0,
 123                      ALL_STD_EXCEPT, FE_INEXACT | FE_UNDERFLOW);
 124                 test(fmal, -LDBL_MIN, LDBL_MIN, 0.0, -0.0,
 125                      ALL_STD_EXCEPT, FE_INEXACT | FE_UNDERFLOW);
 126         }
 127 }
 128
 129 static void
 130 test_infinities(void)
 131 {
 132
 133         testall(INFINITY, 1.0, -1.0, INFINITY, ALL_STD_EXCEPT, 0);
 134         testall(-1.0, INFINITY, 0.0, -INFINITY, ALL_STD_EXCEPT, 0);
 135         testall(0.0, 0.0, INFINITY, INFINITY, ALL_STD_EXCEPT, 0);
 136         testall(1.0, 1.0, INFINITY, INFINITY, ALL_STD_EXCEPT, 0);
 137         testall(1.0, 1.0, -INFINITY, -INFINITY, ALL_STD_EXCEPT, 0);
 138
 139         testall(INFINITY, -INFINITY, 1.0, -INFINITY, ALL_STD_EXCEPT, 0);
 140         testall(INFINITY, INFINITY, 1.0, INFINITY, ALL_STD_EXCEPT, 0);
 141         testall(-INFINITY, -INFINITY, INFINITY, INFINITY, ALL_STD_EXCEPT, 0);
 142
 143         testall(0.0, INFINITY, 1.0, NAN, ALL_STD_EXCEPT, FE_INVALID);
 144         testall(INFINITY, 0.0, -0.0, NAN, ALL_STD_EXCEPT, FE_INVALID);
 145
 146         /* The invalid exception is optional in this case. */
 147         testall(INFINITY, 0.0, NAN, NAN, ALL_STD_EXCEPT & ~FE_INVALID, 0);
 148
 149         testall(INFINITY, INFINITY, -INFINITY, NAN,
 150                 ALL_STD_EXCEPT, FE_INVALID);
 151         testall(-INFINITY, INFINITY, INFINITY, NAN,
 152                 ALL_STD_EXCEPT, FE_INVALID);
 153         testall(INFINITY, -1.0, INFINITY, NAN,
 154                 ALL_STD_EXCEPT, FE_INVALID);
 155
 156         test(fmaf, FLT_MAX, FLT_MAX, -INFINITY, -INFINITY, ALL_STD_EXCEPT, 0);
 157         test(fma, DBL_MAX, DBL_MAX, -INFINITY, -INFINITY, ALL_STD_EXCEPT, 0);
 158         test(fmal, LDBL_MAX, LDBL_MAX, -INFINITY, -INFINITY,
 159              ALL_STD_EXCEPT, 0);
 160         test(fmaf, FLT_MAX, -FLT_MAX, INFINITY, INFINITY, ALL_STD_EXCEPT, 0);
 161         test(fma, DBL_MAX, -DBL_MAX, INFINITY, INFINITY, ALL_STD_EXCEPT, 0);
 162         test(fmal, LDBL_MAX, -LDBL_MAX, INFINITY, INFINITY,
 163              ALL_STD_EXCEPT, 0);
 164 }
 165
 166 static void
 167 test_nans(void)
 168 {
 169
 170         testall(NAN, 0.0, 0.0, NAN, ALL_STD_EXCEPT, 0);
 171         testall(1.0, NAN, 1.0, NAN, ALL_STD_EXCEPT, 0);
 172         testall(1.0, -1.0, NAN, NAN, ALL_STD_EXCEPT, 0);
 173         testall(0.0, 0.0, NAN, NAN, ALL_STD_EXCEPT, 0);
 174         testall(NAN, NAN, NAN, NAN, ALL_STD_EXCEPT, 0);
 175
 176         /* x*y should not raise an inexact/overflow/underflow if z is NaN. */
 177         testall(M_PI, M_PI, NAN, NAN, ALL_STD_EXCEPT, 0);
 178         test(fmaf, FLT_MIN, FLT_MIN, NAN, NAN, ALL_STD_EXCEPT, 0);
 179         test(fma, DBL_MIN, DBL_MIN, NAN, NAN, ALL_STD_EXCEPT, 0);
 180         test(fmal, LDBL_MIN, LDBL_MIN, NAN, NAN, ALL_STD_EXCEPT, 0);
 181         test(fmaf, FLT_MAX, FLT_MAX, NAN, NAN, ALL_STD_EXCEPT, 0);
 182         test(fma, DBL_MAX, DBL_MAX, NAN, NAN, ALL_STD_EXCEPT, 0);
 183         test(fmal, LDBL_MAX, LDBL_MAX, NAN, NAN, ALL_STD_EXCEPT, 0);
 184 }
 185
 186 /*
 187  * Tests for cases where z is very small compared to x*y.
 188  */
 189 static void
 190 test_small_z(void)
 191 {
 192
 193         /* x*y positive, z positive */
 194         if (fegetround() == FE_UPWARD) {
 195                 test(fmaf, 1.0, 1.0, 0x1.0p-100, 1.0 + FLT_EPSILON,
 196                      ALL_STD_EXCEPT, FE_INEXACT);
 197                 test(fma, 1.0, 1.0, 0x1.0p-200, 1.0 + DBL_EPSILON,
 198                      ALL_STD_EXCEPT, FE_INEXACT);
 199                 test(fmal, 1.0, 1.0, 0x1.0p-200, 1.0 + LDBL_EPSILON,
 200                      ALL_STD_EXCEPT, FE_INEXACT);
 201         } else {
 202                 testall(0x1.0p100, 1.0, 0x1.0p-100, 0x1.0p100,
 203                         ALL_STD_EXCEPT, FE_INEXACT);
 204         }
 205
 206         /* x*y negative, z negative */
 207         if (fegetround() == FE_DOWNWARD) {
 208                 test(fmaf, -1.0, 1.0, -0x1.0p-100, -(1.0 + FLT_EPSILON),
 209                      ALL_STD_EXCEPT, FE_INEXACT);
 210                 test(fma, -1.0, 1.0, -0x1.0p-200, -(1.0 + DBL_EPSILON),
 211                      ALL_STD_EXCEPT, FE_INEXACT);
 212                 test(fmal, -1.0, 1.0, -0x1.0p-200, -(1.0 + LDBL_EPSILON),
 213                      ALL_STD_EXCEPT, FE_INEXACT);
 214         } else {
 215                 testall(0x1.0p100, -1.0, -0x1.0p-100, -0x1.0p100,
 216                         ALL_STD_EXCEPT, FE_INEXACT);
 217         }
 218
 219         /* x*y positive, z negative */
 220         if (fegetround() == FE_DOWNWARD || fegetround() == FE_TOWARDZERO) {
 221                 test(fmaf, 1.0, 1.0, -0x1.0p-100, 1.0 - FLT_EPSILON / 2,
 222                      ALL_STD_EXCEPT, FE_INEXACT);
 223                 test(fma, 1.0, 1.0, -0x1.0p-200, 1.0 - DBL_EPSILON / 2,
 224                      ALL_STD_EXCEPT, FE_INEXACT);
 225                 test(fmal, 1.0, 1.0, -0x1.0p-200, 1.0 - LDBL_EPSILON / 2,
 226                      ALL_STD_EXCEPT, FE_INEXACT);
 227         } else {
 228                 testall(0x1.0p100, 1.0, -0x1.0p-100, 0x1.0p100,
 229                         ALL_STD_EXCEPT, FE_INEXACT);
 230         }
 231
 232         /* x*y negative, z positive */
 233         if (fegetround() == FE_UPWARD || fegetround() == FE_TOWARDZERO) {
 234                 test(fmaf, -1.0, 1.0, 0x1.0p-100, -1.0 + FLT_EPSILON / 2,
 235                      ALL_STD_EXCEPT, FE_INEXACT);
 236                 test(fma, -1.0, 1.0, 0x1.0p-200, -1.0 + DBL_EPSILON / 2,
 237                      ALL_STD_EXCEPT, FE_INEXACT);
 238                 test(fmal, -1.0, 1.0, 0x1.0p-200, -1.0 + LDBL_EPSILON / 2,
 239                      ALL_STD_EXCEPT, FE_INEXACT);
 240         } else {
 241                 testall(-0x1.0p100, 1.0, 0x1.0p-100, -0x1.0p100,
 242                         ALL_STD_EXCEPT, FE_INEXACT);
 243         }
 244 }
 245
 246 /*
 247  * Tests for cases where z is very large compared to x*y.
 248  */
 249 static void
 250 test_big_z(void)
 251 {
 252
 253         /* z positive, x*y positive */
 254         if (fegetround() == FE_UPWARD) {
 255                 test(fmaf, 0x1.0p-50, 0x1.0p-50, 1.0, 1.0 + FLT_EPSILON,
 256                      ALL_STD_EXCEPT, FE_INEXACT);
 257                 test(fma, 0x1.0p-100, 0x1.0p-100, 1.0, 1.0 + DBL_EPSILON,
 258                      ALL_STD_EXCEPT, FE_INEXACT);
 259                 test(fmal, 0x1.0p-100, 0x1.0p-100, 1.0, 1.0 + LDBL_EPSILON,
 260                      ALL_STD_EXCEPT, FE_INEXACT);
 261         } else {
 262                 testall(-0x1.0p-50, -0x1.0p-50, 0x1.0p100, 0x1.0p100,
 263                         ALL_STD_EXCEPT, FE_INEXACT);
 264         }
 265
 266         /* z negative, x*y negative */
 267         if (fegetround() == FE_DOWNWARD) {
 268                 test(fmaf, -0x1.0p-50, 0x1.0p-50, -1.0, -(1.0 + FLT_EPSILON),
 269                      ALL_STD_EXCEPT, FE_INEXACT);
 270                 test(fma, -0x1.0p-100, 0x1.0p-100, -1.0, -(1.0 + DBL_EPSILON),
 271                      ALL_STD_EXCEPT, FE_INEXACT);
 272                 test(fmal, -0x1.0p-100, 0x1.0p-100, -1.0, -(1.0 + LDBL_EPSILON),
 273                      ALL_STD_EXCEPT, FE_INEXACT);
 274         } else {
 275                 testall(0x1.0p-50, -0x1.0p-50, -0x1.0p100, -0x1.0p100,
 276                         ALL_STD_EXCEPT, FE_INEXACT);
 277         }
 278
 279         /* z negative, x*y positive */
 280         if (fegetround() == FE_UPWARD || fegetround() == FE_TOWARDZERO) {
 281                 test(fmaf, -0x1.0p-50, -0x1.0p-50, -1.0,
 282                      -1.0 + FLT_EPSILON / 2, ALL_STD_EXCEPT, FE_INEXACT);
 283                 test(fma, -0x1.0p-100, -0x1.0p-100, -1.0,
 284                      -1.0 + DBL_EPSILON / 2, ALL_STD_EXCEPT, FE_INEXACT);
 285                 test(fmal, -0x1.0p-100, -0x1.0p-100, -1.0,
 286                      -1.0 + LDBL_EPSILON / 2, ALL_STD_EXCEPT, FE_INEXACT);
 287         } else {
 288                 testall(0x1.0p-50, 0x1.0p-50, -0x1.0p100, -0x1.0p100,
 289                         ALL_STD_EXCEPT, FE_INEXACT);
 290         }
 291
 292         /* z positive, x*y negative */
 293         if (fegetround() == FE_DOWNWARD || fegetround() == FE_TOWARDZERO) {
 294                 test(fmaf, 0x1.0p-50, -0x1.0p-50, 1.0, 1.0 - FLT_EPSILON / 2,
 295                      ALL_STD_EXCEPT, FE_INEXACT);
 296                 test(fma, 0x1.0p-100, -0x1.0p-100, 1.0, 1.0 - DBL_EPSILON / 2,
 297                      ALL_STD_EXCEPT, FE_INEXACT);
 298                 test(fmal, 0x1.0p-100, -0x1.0p-100, 1.0, 1.0 - LDBL_EPSILON / 2,
 299                      ALL_STD_EXCEPT, FE_INEXACT);
 300         } else {
 301                 testall(-0x1.0p-50, 0x1.0p-50, 0x1.0p100, 0x1.0p100,
 302                         ALL_STD_EXCEPT, FE_INEXACT);
 303         }
 304 }
 305
 306 static void
 307 test_accuracy(void)
 308 {
 309
 310         /* ilogb(x*y) - ilogb(z) = 20 */
 311         testrnd(fmaf, -0x1.c139d8p-51, -0x1.600e7ap32, 0x1.26558cp-38,
 312                 0x1.34e48ap-18, 0x1.34e48cp-18, 0x1.34e48ap-18, 0x1.34e48ap-18,
 313                 ALL_STD_EXCEPT, FE_INEXACT);
 314         testrnd(fma, -0x1.c139d7b84f1a3p-51, -0x1.600e7a2a16484p32,
 315                 0x1.26558cac31580p-38, 0x1.34e48a78aae97p-18,
 316                 0x1.34e48a78aae97p-18, 0x1.34e48a78aae96p-18,
 317                 0x1.34e48a78aae96p-18, ALL_STD_EXCEPT, FE_INEXACT);
 318 #if LDBL_MANT_DIG == 113
 319         testrnd(fmal, -0x1.c139d7b84f1a3079263afcc5bae3p-51L,
 320                 -0x1.600e7a2a164840edbe2e7d301a72p32L,
 321                 0x1.26558cac315807eb07e448042101p-38L,
 322                 0x1.34e48a78aae96c76ed36077dd387p-18L,
 323                 0x1.34e48a78aae96c76ed36077dd388p-18L,
 324                 0x1.34e48a78aae96c76ed36077dd387p-18L,
 325                 0x1.34e48a78aae96c76ed36077dd387p-18L,
 326                 ALL_STD_EXCEPT, FE_INEXACT);
 327 #elif LDBL_MANT_DIG == 64
 328         testrnd(fmal, -0x1.c139d7b84f1a307ap-51L, -0x1.600e7a2a164840eep32L,
 329                 0x1.26558cac315807ecp-38L, 0x1.34e48a78aae96c78p-18L,
 330                 0x1.34e48a78aae96c78p-18L, 0x1.34e48a78aae96c76p-18L,
 331                 0x1.34e48a78aae96c76p-18L, ALL_STD_EXCEPT, FE_INEXACT);
 332 #elif LDBL_MANT_DIG == 53
 333         testrnd(fmal, -0x1.c139d7b84f1a3p-51L, -0x1.600e7a2a16484p32L,
 334                 0x1.26558cac31580p-38L, 0x1.34e48a78aae97p-18L,
 335                 0x1.34e48a78aae97p-18L, 0x1.34e48a78aae96p-18L,
 336                 0x1.34e48a78aae96p-18L, ALL_STD_EXCEPT, FE_INEXACT);
 337 #endif
 338
 339         /* ilogb(x*y) - ilogb(z) = -40 */
 340         testrnd(fmaf, 0x1.98210ap53, 0x1.9556acp-24, 0x1.d87da4p70,
 341                 0x1.d87da4p70, 0x1.d87da6p70, 0x1.d87da4p70, 0x1.d87da4p70,
 342                 ALL_STD_EXCEPT, FE_INEXACT);
 343         testrnd(fma, 0x1.98210ac83fe2bp53, 0x1.9556ac1475f0fp-24,
 344                 0x1.d87da3aafc60ep70, 0x1.d87da3aafda40p70,
 345                 0x1.d87da3aafda40p70, 0x1.d87da3aafda3fp70,
 346                 0x1.d87da3aafda3fp70, ALL_STD_EXCEPT, FE_INEXACT);
 347 #if LDBL_MANT_DIG == 113
 348         testrnd(fmal, 0x1.98210ac83fe2a8f65b6278b74cebp53L,
 349                 0x1.9556ac1475f0f28968b61d0de65ap-24L,
 350                 0x1.d87da3aafc60d830aa4c6d73b749p70L,
 351                 0x1.d87da3aafda3f36a69eb86488224p70L,
 352                 0x1.d87da3aafda3f36a69eb86488225p70L,
 353                 0x1.d87da3aafda3f36a69eb86488224p70L,
 354                 0x1.d87da3aafda3f36a69eb86488224p70L,
 355                 ALL_STD_EXCEPT, FE_INEXACT);
 356 #elif LDBL_MANT_DIG == 64
 357         testrnd(fmal, 0x1.98210ac83fe2a8f6p53L, 0x1.9556ac1475f0f28ap-24L,
 358                 0x1.d87da3aafc60d83p70L, 0x1.d87da3aafda3f36ap70L,
 359                 0x1.d87da3aafda3f36ap70L, 0x1.d87da3aafda3f368p70L,
 360                 0x1.d87da3aafda3f368p70L, ALL_STD_EXCEPT, FE_INEXACT);
 361 #elif LDBL_MANT_DIG == 53
 362         testrnd(fmal, 0x1.98210ac83fe2bp53L, 0x1.9556ac1475f0fp-24L,
 363                 0x1.d87da3aafc60ep70L, 0x1.d87da3aafda40p70L,
 364                 0x1.d87da3aafda40p70L, 0x1.d87da3aafda3fp70L,
 365                 0x1.d87da3aafda3fp70L, ALL_STD_EXCEPT, FE_INEXACT);
 366 #endif
 367 }
 368
 369 static void
 370 test_double_rounding(void)
 371 {
 372
 373         /*
 374          *     a =  0x1.8000000000001p0
 375          *     b =  0x1.8000000000001p0
 376          *     c = -0x0.0000000000000000000000000080...1p+1
 377          * a * b =  0x1.2000000000001800000000000080p+1
 378          *
 379          * The correct behavior is to round DOWN to 0x1.2000000000001p+1 in
 380          * round-to-nearest mode.  An implementation that computes a*b+c in
 381          * double+double precision, however, will get 0x1.20000000000018p+1,
 382          * and then round UP.
 383          */
 384         fesetround(FE_TONEAREST);
 385         test(fma, 0x1.8000000000001p0, 0x1.8000000000001p0,
 386              -0x1.0000000000001p-104, 0x1.2000000000001p+1,
 387              ALL_STD_EXCEPT, FE_INEXACT);
 388         fesetround(FE_DOWNWARD);
 389         test(fma, 0x1.8000000000001p0, 0x1.8000000000001p0,
 390              -0x1.0000000000001p-104, 0x1.2000000000001p+1,
 391              ALL_STD_EXCEPT, FE_INEXACT);
 392         fesetround(FE_UPWARD);
 393         test(fma, 0x1.8000000000001p0, 0x1.8000000000001p0,
 394              -0x1.0000000000001p-104, 0x1.2000000000002p+1,
 395              ALL_STD_EXCEPT, FE_INEXACT);
 396
 397         fesetround(FE_TONEAREST);
 398         test(fmaf, 0x1.800002p+0, 0x1.800002p+0, -0x1.000002p-46, 0x1.200002p+1,
 399              ALL_STD_EXCEPT, FE_INEXACT);
 400         fesetround(FE_DOWNWARD);
 401         test(fmaf, 0x1.800002p+0, 0x1.800002p+0, -0x1.000002p-46, 0x1.200002p+1,
 402              ALL_STD_EXCEPT, FE_INEXACT);
 403         fesetround(FE_UPWARD);
 404         test(fmaf, 0x1.800002p+0, 0x1.800002p+0, -0x1.000002p-46, 0x1.200004p+1,
 405              ALL_STD_EXCEPT, FE_INEXACT);
 406
 407         fesetround(FE_TONEAREST);
 408 #if LDBL_MANT_DIG == 64
 409         test(fmal, 0x1.4p+0L, 0x1.0000000000000004p+0L, 0x1p-128L,
 410              0x1.4000000000000006p+0L, ALL_STD_EXCEPT, FE_INEXACT);
 411 #elif LDBL_MANT_DIG == 113
 412         /* XXX untested test */
 413         test(fmal, 0x1.4p+0L, 0x1.0000000000000000000000000002p+0L, 0x1p-224L,
 414              0x1.4000000000000000000000000003p+0L, ALL_STD_EXCEPT, FE_INEXACT);
 415 #endif
 416
 417 }
 418
 419 int
 420 main(int argc, char *argv[])
 421 {
 422         int rmodes[] = { FE_TONEAREST, FE_UPWARD, FE_DOWNWARD, FE_TOWARDZERO };
 423         int i, j;
 424
 425         j = 1;
 426
 427 #if defined(__i386__)
 428         printf("1..0 # SKIP all testcases fail on i386\n");
 429         exit(0);
 430 #endif
 431         printf("1..19\n");
 432
 433         for (i = 0; i < nitems(rmodes); i++, j++) {
 434                 printf("rmode = %d\n", rmodes[i]);
 435                 fesetround(rmodes[i]);
 436                 test_zeroes();
 437                 printf("ok %d - fma zeroes\n", j);
 438         }
 439
 440         for (i = 0; i < nitems(rmodes); i++, j++) {
 441                 printf("rmode = %d\n", rmodes[i]);
 442                 fesetround(rmodes[i]);
 443                 test_infinities();
 444                 printf("ok %d - fma infinities\n", j);
 445         }
 446
 447         fesetround(FE_TONEAREST);
 448         test_nans();
 449         printf("ok %d - fma NaNs\n", j);
 450         j++;
 451
 452         for (i = 0; i < nitems(rmodes); i++, j++) {
 453                 printf("rmode = %d\n", rmodes[i]);
 454                 fesetround(rmodes[i]);
 455                 test_small_z();
 456                 printf("ok %d - fma small z\n", j);
 457         }
 458
 459         for (i = 0; i < nitems(rmodes); i++, j++) {
 460                 printf("rmode = %d\n", rmodes[i]);
 461                 fesetround(rmodes[i]);
 462                 test_big_z();
 463                 printf("ok %d - fma big z\n", j);
 464         }
 465
 466         fesetround(FE_TONEAREST);
 467         test_accuracy();
 468         printf("ok %d - fma accuracy\n", j);
 469         j++;
 470
 471         test_double_rounding();
 472         printf("ok %d - fma double rounding\n", j);
 473         j++;
 474
 475         /*
 476          * TODO:
 477          * - Tests for subnormals
 478          * - Cancellation tests (e.g., z = (double)x*y, but x*y is inexact)
 479          */
 480
 481         return (0);
 482 }