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 * $FreeBSD: src/tools/regression/lib/msun/test-fma.c,v 1.6 2013/05/28 00:27:56 svnexp Exp $ 27 */ 28 29 /* 30 * Tests for fma{,f,l}(). 31 */ 32 33 #include <assert.h> 34 #include <fenv.h> 35 #include <float.h> 36 #include <math.h> 37 #include <stdio.h> 38 39 #define ALL_STD_EXCEPT (FE_DIVBYZERO | FE_INEXACT | FE_INVALID | \ 40 FE_OVERFLOW | FE_UNDERFLOW) 41 42 #pragma STDC FENV_ACCESS ON 43 44 /* 45 * Test that a function returns the correct value and sets the 46 * exception flags correctly. The exceptmask specifies which 47 * exceptions we should check. We need to be lenient for several 48 * reasons, but mainly because on some architectures it's impossible 49 * to raise FE_OVERFLOW without raising FE_INEXACT. 50 * 51 * These are macros instead of functions so that assert provides more 52 * meaningful error messages. 53 */ 54 #define test(func, x, y, z, result, exceptmask, excepts) do { \ 55 assert(feclearexcept(FE_ALL_EXCEPT) == 0); \ 56 assert(fpequal((func)((x), (y), (z)), (result))); \ 57 assert(((func), fetestexcept(exceptmask) == (excepts))); \ 58 } while (0) 59 60 #define testall(x, y, z, result, exceptmask, excepts) do { \ 61 test(fma, (x), (y), (z), (double)(result), (exceptmask), (excepts)); \ 62 test(fmaf, (x), (y), (z), (float)(result), (exceptmask), (excepts)); \ 63 test(fmal, (x), (y), (z), (result), (exceptmask), (excepts)); \ 64 } while (0) 65 66 /* Test in all rounding modes. */ 67 #define testrnd(func, x, y, z, rn, ru, rd, rz, exceptmask, excepts) do { \ 68 fesetround(FE_TONEAREST); \ 69 test((func), (x), (y), (z), (rn), (exceptmask), (excepts)); \ 70 fesetround(FE_UPWARD); \ 71 test((func), (x), (y), (z), (ru), (exceptmask), (excepts)); \ 72 fesetround(FE_DOWNWARD); \ 73 test((func), (x), (y), (z), (rd), (exceptmask), (excepts)); \ 74 fesetround(FE_TOWARDZERO); \ 75 test((func), (x), (y), (z), (rz), (exceptmask), (excepts)); \ 76 } while (0) 77 78 /* 79 * This is needed because clang constant-folds fma in ways that are incorrect 80 * in rounding modes other than FE_TONEAREST. 81 */ 82 volatile double one = 1.0; 83 84 /* 85 * Determine whether x and y are equal, with two special rules: 86 * +0.0 != -0.0 87 * NaN == NaN 88 */ 89 int 90 fpequal(long double x, long double y) 91 { 92 93 return ((x == y && !signbit(x) == !signbit(y)) 94 || (isnan(x) && isnan(y))); 95 } 96 97 static void 98 test_zeroes(void) 99 { 100 const int rd = (fegetround() == FE_DOWNWARD); 101 102 testall(0.0, 0.0, 0.0, 0.0, ALL_STD_EXCEPT, 0); 103 testall(1.0, 0.0, 0.0, 0.0, ALL_STD_EXCEPT, 0); 104 testall(0.0, 1.0, 0.0, 0.0, ALL_STD_EXCEPT, 0); 105 testall(0.0, 0.0, 1.0, 1.0, ALL_STD_EXCEPT, 0); 106 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 testall(-0.0, -0.0, 0.0, 0.0, ALL_STD_EXCEPT, 0); 110 testall(0.0, 0.0, -0.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); 111 testall(-0.0, -0.0, -0.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); 112 113 testall(-0.0, 0.0, -0.0, -0.0, ALL_STD_EXCEPT, 0); 114 testall(0.0, -0.0, -0.0, -0.0, ALL_STD_EXCEPT, 0); 115 116 testall(-one, one, one, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); 117 testall(one, -one, one, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); 118 testall(-one, -one, -one, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); 119 120 switch (fegetround()) { 121 case FE_TONEAREST: 122 case FE_TOWARDZERO: 123 test(fmaf, -FLT_MIN, FLT_MIN, 0.0, -0.0, 124 ALL_STD_EXCEPT, FE_INEXACT | FE_UNDERFLOW); 125 test(fma, -DBL_MIN, DBL_MIN, 0.0, -0.0, 126 ALL_STD_EXCEPT, FE_INEXACT | FE_UNDERFLOW); 127 test(fmal, -LDBL_MIN, LDBL_MIN, 0.0, -0.0, 128 ALL_STD_EXCEPT, FE_INEXACT | FE_UNDERFLOW); 129 } 130 } 131 132 static void 133 test_infinities(void) 134 { 135 136 testall(INFINITY, 1.0, -1.0, INFINITY, ALL_STD_EXCEPT, 0); 137 testall(-1.0, INFINITY, 0.0, -INFINITY, ALL_STD_EXCEPT, 0); 138 testall(0.0, 0.0, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); 139 testall(1.0, 1.0, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); 140 testall(1.0, 1.0, -INFINITY, -INFINITY, ALL_STD_EXCEPT, 0); 141 142 testall(INFINITY, -INFINITY, 1.0, -INFINITY, ALL_STD_EXCEPT, 0); 143 testall(INFINITY, INFINITY, 1.0, INFINITY, ALL_STD_EXCEPT, 0); 144 testall(-INFINITY, -INFINITY, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); 145 146 testall(0.0, INFINITY, 1.0, NAN, ALL_STD_EXCEPT, FE_INVALID); 147 testall(INFINITY, 0.0, -0.0, NAN, ALL_STD_EXCEPT, FE_INVALID); 148 149 /* The invalid exception is optional in this case. */ 150 testall(INFINITY, 0.0, NAN, NAN, ALL_STD_EXCEPT & ~FE_INVALID, 0); 151 152 testall(INFINITY, INFINITY, -INFINITY, NAN, 153 ALL_STD_EXCEPT, FE_INVALID); 154 testall(-INFINITY, INFINITY, INFINITY, NAN, 155 ALL_STD_EXCEPT, FE_INVALID); 156 testall(INFINITY, -1.0, INFINITY, NAN, 157 ALL_STD_EXCEPT, FE_INVALID); 158 159 test(fmaf, FLT_MAX, FLT_MAX, -INFINITY, -INFINITY, ALL_STD_EXCEPT, 0); 160 test(fma, DBL_MAX, DBL_MAX, -INFINITY, -INFINITY, ALL_STD_EXCEPT, 0); 161 test(fmal, LDBL_MAX, LDBL_MAX, -INFINITY, -INFINITY, 162 ALL_STD_EXCEPT, 0); 163 test(fmaf, FLT_MAX, -FLT_MAX, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); 164 test(fma, DBL_MAX, -DBL_MAX, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); 165 test(fmal, LDBL_MAX, -LDBL_MAX, INFINITY, INFINITY, 166 ALL_STD_EXCEPT, 0); 167 } 168 169 static void 170 test_nans(void) 171 { 172 173 testall(NAN, 0.0, 0.0, NAN, ALL_STD_EXCEPT, 0); 174 testall(1.0, NAN, 1.0, NAN, ALL_STD_EXCEPT, 0); 175 testall(1.0, -1.0, NAN, NAN, ALL_STD_EXCEPT, 0); 176 testall(0.0, 0.0, NAN, NAN, ALL_STD_EXCEPT, 0); 177 testall(NAN, NAN, NAN, NAN, ALL_STD_EXCEPT, 0); 178 179 /* x*y should not raise an inexact/overflow/underflow if z is NaN. */ 180 testall(M_PI, M_PI, NAN, NAN, ALL_STD_EXCEPT, 0); 181 test(fmaf, FLT_MIN, FLT_MIN, NAN, NAN, ALL_STD_EXCEPT, 0); 182 test(fma, DBL_MIN, DBL_MIN, NAN, NAN, ALL_STD_EXCEPT, 0); 183 test(fmal, LDBL_MIN, LDBL_MIN, NAN, NAN, ALL_STD_EXCEPT, 0); 184 test(fmaf, FLT_MAX, FLT_MAX, NAN, NAN, ALL_STD_EXCEPT, 0); 185 test(fma, DBL_MAX, DBL_MAX, NAN, NAN, ALL_STD_EXCEPT, 0); 186 test(fmal, LDBL_MAX, LDBL_MAX, NAN, NAN, ALL_STD_EXCEPT, 0); 187 } 188 189 /* 190 * Tests for cases where z is very small compared to x*y. 191 */ 192 static void 193 test_small_z(void) 194 { 195 196 /* x*y positive, z positive */ 197 if (fegetround() == FE_UPWARD) { 198 test(fmaf, one, one, 0x1.0p-100, 1.0 + FLT_EPSILON, 199 ALL_STD_EXCEPT, FE_INEXACT); 200 test(fma, one, one, 0x1.0p-200, 1.0 + DBL_EPSILON, 201 ALL_STD_EXCEPT, FE_INEXACT); 202 test(fmal, one, one, 0x1.0p-200, 1.0 + LDBL_EPSILON, 203 ALL_STD_EXCEPT, FE_INEXACT); 204 } else { 205 testall(0x1.0p100, one, 0x1.0p-100, 0x1.0p100, 206 ALL_STD_EXCEPT, FE_INEXACT); 207 } 208 209 /* x*y negative, z negative */ 210 if (fegetround() == FE_DOWNWARD) { 211 test(fmaf, -one, one, -0x1.0p-100, -(1.0 + FLT_EPSILON), 212 ALL_STD_EXCEPT, FE_INEXACT); 213 test(fma, -one, one, -0x1.0p-200, -(1.0 + DBL_EPSILON), 214 ALL_STD_EXCEPT, FE_INEXACT); 215 test(fmal, -one, one, -0x1.0p-200, -(1.0 + LDBL_EPSILON), 216 ALL_STD_EXCEPT, FE_INEXACT); 217 } else { 218 testall(0x1.0p100, -one, -0x1.0p-100, -0x1.0p100, 219 ALL_STD_EXCEPT, FE_INEXACT); 220 } 221 222 /* x*y positive, z negative */ 223 if (fegetround() == FE_DOWNWARD || fegetround() == FE_TOWARDZERO) { 224 test(fmaf, one, one, -0x1.0p-100, 1.0 - FLT_EPSILON / 2, 225 ALL_STD_EXCEPT, FE_INEXACT); 226 test(fma, one, one, -0x1.0p-200, 1.0 - DBL_EPSILON / 2, 227 ALL_STD_EXCEPT, FE_INEXACT); 228 test(fmal, one, one, -0x1.0p-200, 1.0 - LDBL_EPSILON / 2, 229 ALL_STD_EXCEPT, FE_INEXACT); 230 } else { 231 testall(0x1.0p100, one, -0x1.0p-100, 0x1.0p100, 232 ALL_STD_EXCEPT, FE_INEXACT); 233 } 234 235 /* x*y negative, z positive */ 236 if (fegetround() == FE_UPWARD || fegetround() == FE_TOWARDZERO) { 237 test(fmaf, -one, one, 0x1.0p-100, -1.0 + FLT_EPSILON / 2, 238 ALL_STD_EXCEPT, FE_INEXACT); 239 test(fma, -one, one, 0x1.0p-200, -1.0 + DBL_EPSILON / 2, 240 ALL_STD_EXCEPT, FE_INEXACT); 241 test(fmal, -one, one, 0x1.0p-200, -1.0 + LDBL_EPSILON / 2, 242 ALL_STD_EXCEPT, FE_INEXACT); 243 } else { 244 testall(-0x1.0p100, one, 0x1.0p-100, -0x1.0p100, 245 ALL_STD_EXCEPT, FE_INEXACT); 246 } 247 } 248 249 /* 250 * Tests for cases where z is very large compared to x*y. 251 */ 252 static void 253 test_big_z(void) 254 { 255 256 /* z positive, x*y positive */ 257 if (fegetround() == FE_UPWARD) { 258 test(fmaf, 0x1.0p-50, 0x1.0p-50, 1.0, 1.0 + FLT_EPSILON, 259 ALL_STD_EXCEPT, FE_INEXACT); 260 test(fma, 0x1.0p-100, 0x1.0p-100, 1.0, 1.0 + DBL_EPSILON, 261 ALL_STD_EXCEPT, FE_INEXACT); 262 test(fmal, 0x1.0p-100, 0x1.0p-100, 1.0, 1.0 + LDBL_EPSILON, 263 ALL_STD_EXCEPT, FE_INEXACT); 264 } else { 265 testall(-0x1.0p-50, -0x1.0p-50, 0x1.0p100, 0x1.0p100, 266 ALL_STD_EXCEPT, FE_INEXACT); 267 } 268 269 /* z negative, x*y negative */ 270 if (fegetround() == FE_DOWNWARD) { 271 test(fmaf, -0x1.0p-50, 0x1.0p-50, -1.0, -(1.0 + FLT_EPSILON), 272 ALL_STD_EXCEPT, FE_INEXACT); 273 test(fma, -0x1.0p-100, 0x1.0p-100, -1.0, -(1.0 + DBL_EPSILON), 274 ALL_STD_EXCEPT, FE_INEXACT); 275 test(fmal, -0x1.0p-100, 0x1.0p-100, -1.0, -(1.0 + LDBL_EPSILON), 276 ALL_STD_EXCEPT, FE_INEXACT); 277 } else { 278 testall(0x1.0p-50, -0x1.0p-50, -0x1.0p100, -0x1.0p100, 279 ALL_STD_EXCEPT, FE_INEXACT); 280 } 281 282 /* z negative, x*y positive */ 283 if (fegetround() == FE_UPWARD || fegetround() == FE_TOWARDZERO) { 284 test(fmaf, -0x1.0p-50, -0x1.0p-50, -1.0, 285 -1.0 + FLT_EPSILON / 2, ALL_STD_EXCEPT, FE_INEXACT); 286 test(fma, -0x1.0p-100, -0x1.0p-100, -1.0, 287 -1.0 + DBL_EPSILON / 2, ALL_STD_EXCEPT, FE_INEXACT); 288 test(fmal, -0x1.0p-100, -0x1.0p-100, -1.0, 289 -1.0 + LDBL_EPSILON / 2, ALL_STD_EXCEPT, FE_INEXACT); 290 } else { 291 testall(0x1.0p-50, 0x1.0p-50, -0x1.0p100, -0x1.0p100, 292 ALL_STD_EXCEPT, FE_INEXACT); 293 } 294 295 /* z positive, x*y negative */ 296 if (fegetround() == FE_DOWNWARD || fegetround() == FE_TOWARDZERO) { 297 test(fmaf, 0x1.0p-50, -0x1.0p-50, 1.0, 1.0 - FLT_EPSILON / 2, 298 ALL_STD_EXCEPT, FE_INEXACT); 299 test(fma, 0x1.0p-100, -0x1.0p-100, 1.0, 1.0 - DBL_EPSILON / 2, 300 ALL_STD_EXCEPT, FE_INEXACT); 301 test(fmal, 0x1.0p-100, -0x1.0p-100, 1.0, 1.0 - LDBL_EPSILON / 2, 302 ALL_STD_EXCEPT, FE_INEXACT); 303 } else { 304 testall(-0x1.0p-50, 0x1.0p-50, 0x1.0p100, 0x1.0p100, 305 ALL_STD_EXCEPT, FE_INEXACT); 306 } 307 } 308 309 static void 310 test_accuracy(void) 311 { 312 313 /* ilogb(x*y) - ilogb(z) = 20 */ 314 testrnd(fmaf, -0x1.c139d8p-51, -0x1.600e7ap32, 0x1.26558cp-38, 315 0x1.34e48ap-18, 0x1.34e48cp-18, 0x1.34e48ap-18, 0x1.34e48ap-18, 316 ALL_STD_EXCEPT, FE_INEXACT); 317 testrnd(fma, -0x1.c139d7b84f1a3p-51, -0x1.600e7a2a16484p32, 318 0x1.26558cac31580p-38, 0x1.34e48a78aae97p-18, 319 0x1.34e48a78aae97p-18, 0x1.34e48a78aae96p-18, 320 0x1.34e48a78aae96p-18, ALL_STD_EXCEPT, FE_INEXACT); 321 #if LDBL_MANT_DIG == 113 322 testrnd(fmal, -0x1.c139d7b84f1a3079263afcc5bae3p-51L, 323 -0x1.600e7a2a164840edbe2e7d301a72p32L, 324 0x1.26558cac315807eb07e448042101p-38L, 325 0x1.34e48a78aae96c76ed36077dd387p-18L, 326 0x1.34e48a78aae96c76ed36077dd388p-18L, 327 0x1.34e48a78aae96c76ed36077dd387p-18L, 328 0x1.34e48a78aae96c76ed36077dd387p-18L, 329 ALL_STD_EXCEPT, FE_INEXACT); 330 #elif LDBL_MANT_DIG == 64 331 testrnd(fmal, -0x1.c139d7b84f1a307ap-51L, -0x1.600e7a2a164840eep32L, 332 0x1.26558cac315807ecp-38L, 0x1.34e48a78aae96c78p-18L, 333 0x1.34e48a78aae96c78p-18L, 0x1.34e48a78aae96c76p-18L, 334 0x1.34e48a78aae96c76p-18L, ALL_STD_EXCEPT, FE_INEXACT); 335 #elif LDBL_MANT_DIG == 53 336 testrnd(fmal, -0x1.c139d7b84f1a3p-51L, -0x1.600e7a2a16484p32L, 337 0x1.26558cac31580p-38L, 0x1.34e48a78aae97p-18L, 338 0x1.34e48a78aae97p-18L, 0x1.34e48a78aae96p-18L, 339 0x1.34e48a78aae96p-18L, ALL_STD_EXCEPT, FE_INEXACT); 340 #endif 341 342 /* ilogb(x*y) - ilogb(z) = -40 */ 343 testrnd(fmaf, 0x1.98210ap53, 0x1.9556acp-24, 0x1.d87da4p70, 344 0x1.d87da4p70, 0x1.d87da6p70, 0x1.d87da4p70, 0x1.d87da4p70, 345 ALL_STD_EXCEPT, FE_INEXACT); 346 testrnd(fma, 0x1.98210ac83fe2bp53, 0x1.9556ac1475f0fp-24, 347 0x1.d87da3aafc60ep70, 0x1.d87da3aafda40p70, 348 0x1.d87da3aafda40p70, 0x1.d87da3aafda3fp70, 349 0x1.d87da3aafda3fp70, ALL_STD_EXCEPT, FE_INEXACT); 350 #if LDBL_MANT_DIG == 113 351 testrnd(fmal, 0x1.98210ac83fe2a8f65b6278b74cebp53L, 352 0x1.9556ac1475f0f28968b61d0de65ap-24L, 353 0x1.d87da3aafc60d830aa4c6d73b749p70L, 354 0x1.d87da3aafda3f36a69eb86488224p70L, 355 0x1.d87da3aafda3f36a69eb86488225p70L, 356 0x1.d87da3aafda3f36a69eb86488224p70L, 357 0x1.d87da3aafda3f36a69eb86488224p70L, 358 ALL_STD_EXCEPT, FE_INEXACT); 359 #elif LDBL_MANT_DIG == 64 360 testrnd(fmal, 0x1.98210ac83fe2a8f6p53L, 0x1.9556ac1475f0f28ap-24L, 361 0x1.d87da3aafc60d83p70L, 0x1.d87da3aafda3f36ap70L, 362 0x1.d87da3aafda3f36ap70L, 0x1.d87da3aafda3f368p70L, 363 0x1.d87da3aafda3f368p70L, ALL_STD_EXCEPT, FE_INEXACT); 364 #elif LDBL_MANT_DIG == 53 365 testrnd(fmal, 0x1.98210ac83fe2bp53L, 0x1.9556ac1475f0fp-24L, 366 0x1.d87da3aafc60ep70L, 0x1.d87da3aafda40p70L, 367 0x1.d87da3aafda40p70L, 0x1.d87da3aafda3fp70L, 368 0x1.d87da3aafda3fp70L, ALL_STD_EXCEPT, FE_INEXACT); 369 #endif 370 371 /* ilogb(x*y) - ilogb(z) = 0 */ 372 testrnd(fmaf, 0x1.31ad02p+100, 0x1.2fbf7ap-42, -0x1.c3e106p+58, 373 -0x1.64c27cp+56, -0x1.64c27ap+56, -0x1.64c27cp+56, 374 -0x1.64c27ap+56, ALL_STD_EXCEPT, FE_INEXACT); 375 testrnd(fma, 0x1.31ad012ede8aap+100, 0x1.2fbf79c839067p-42, 376 -0x1.c3e106929056ep+58, -0x1.64c282b970a5fp+56, 377 -0x1.64c282b970a5ep+56, -0x1.64c282b970a5fp+56, 378 -0x1.64c282b970a5ep+56, ALL_STD_EXCEPT, FE_INEXACT); 379 #if LDBL_MANT_DIG == 113 380 testrnd(fmal, 0x1.31ad012ede8aa282fa1c19376d16p+100L, 381 0x1.2fbf79c839066f0f5c68f6d2e814p-42L, 382 -0x1.c3e106929056ec19de72bfe64215p+58L, 383 -0x1.64c282b970a612598fc025ca8cddp+56L, 384 -0x1.64c282b970a612598fc025ca8cddp+56L, 385 -0x1.64c282b970a612598fc025ca8cdep+56L, 386 -0x1.64c282b970a612598fc025ca8cddp+56L, 387 ALL_STD_EXCEPT, FE_INEXACT); 388 #elif LDBL_MANT_DIG == 64 389 testrnd(fmal, 0x1.31ad012ede8aa4eap+100L, 0x1.2fbf79c839066aeap-42L, 390 -0x1.c3e106929056e61p+58L, -0x1.64c282b970a60298p+56L, 391 -0x1.64c282b970a60298p+56L, -0x1.64c282b970a6029ap+56L, 392 -0x1.64c282b970a60298p+56L, ALL_STD_EXCEPT, FE_INEXACT); 393 #elif LDBL_MANT_DIG == 53 394 testrnd(fmal, 0x1.31ad012ede8aap+100L, 0x1.2fbf79c839067p-42L, 395 -0x1.c3e106929056ep+58L, -0x1.64c282b970a5fp+56L, 396 -0x1.64c282b970a5ep+56L, -0x1.64c282b970a5fp+56L, 397 -0x1.64c282b970a5ep+56L, ALL_STD_EXCEPT, FE_INEXACT); 398 #endif 399 400 /* x*y (rounded) ~= -z */ 401 /* XXX spurious inexact exceptions */ 402 testrnd(fmaf, 0x1.bbffeep-30, -0x1.1d164cp-74, 0x1.ee7296p-104, 403 -0x1.c46ea8p-128, -0x1.c46ea8p-128, -0x1.c46ea8p-128, 404 -0x1.c46ea8p-128, ALL_STD_EXCEPT & ~FE_INEXACT, 0); 405 testrnd(fma, 0x1.bbffeea6fc7d6p-30, 0x1.1d164c6cbf078p-74, 406 -0x1.ee72993aff948p-104, -0x1.71f72ac7d9d8p-159, 407 -0x1.71f72ac7d9d8p-159, -0x1.71f72ac7d9d8p-159, 408 -0x1.71f72ac7d9d8p-159, ALL_STD_EXCEPT & ~FE_INEXACT, 0); 409 #if LDBL_MANT_DIG == 113 410 testrnd(fmal, 0x1.bbffeea6fc7d65927d147f437675p-30L, 411 0x1.1d164c6cbf078b7a22607d1cd6a2p-74L, 412 -0x1.ee72993aff94973876031bec0944p-104L, 413 0x1.64e086175b3a2adc36e607058814p-217L, 414 0x1.64e086175b3a2adc36e607058814p-217L, 415 0x1.64e086175b3a2adc36e607058814p-217L, 416 0x1.64e086175b3a2adc36e607058814p-217L, 417 ALL_STD_EXCEPT & ~FE_INEXACT, 0); 418 #elif LDBL_MANT_DIG == 64 419 testrnd(fmal, 0x1.bbffeea6fc7d6592p-30L, 0x1.1d164c6cbf078b7ap-74L, 420 -0x1.ee72993aff949736p-104L, 0x1.af190e7a1ee6ad94p-168L, 421 0x1.af190e7a1ee6ad94p-168L, 0x1.af190e7a1ee6ad94p-168L, 422 0x1.af190e7a1ee6ad94p-168L, ALL_STD_EXCEPT & ~FE_INEXACT, 0); 423 #elif LDBL_MANT_DIG == 53 424 testrnd(fmal, 0x1.bbffeea6fc7d6p-30L, 0x1.1d164c6cbf078p-74L, 425 -0x1.ee72993aff948p-104L, -0x1.71f72ac7d9d8p-159L, 426 -0x1.71f72ac7d9d8p-159L, -0x1.71f72ac7d9d8p-159L, 427 -0x1.71f72ac7d9d8p-159L, ALL_STD_EXCEPT & ~FE_INEXACT, 0); 428 #endif 429 } 430 431 static void 432 test_double_rounding(void) 433 { 434 435 /* 436 * a = 0x1.8000000000001p0 437 * b = 0x1.8000000000001p0 438 * c = -0x0.0000000000000000000000000080...1p+1 439 * a * b = 0x1.2000000000001800000000000080p+1 440 * 441 * The correct behavior is to round DOWN to 0x1.2000000000001p+1 in 442 * round-to-nearest mode. An implementation that computes a*b+c in 443 * double+double precision, however, will get 0x1.20000000000018p+1, 444 * and then round UP. 445 */ 446 fesetround(FE_TONEAREST); 447 test(fma, 0x1.8000000000001p0, 0x1.8000000000001p0, 448 -0x1.0000000000001p-104, 0x1.2000000000001p+1, 449 ALL_STD_EXCEPT, FE_INEXACT); 450 fesetround(FE_DOWNWARD); 451 test(fma, 0x1.8000000000001p0, 0x1.8000000000001p0, 452 -0x1.0000000000001p-104, 0x1.2000000000001p+1, 453 ALL_STD_EXCEPT, FE_INEXACT); 454 fesetround(FE_UPWARD); 455 test(fma, 0x1.8000000000001p0, 0x1.8000000000001p0, 456 -0x1.0000000000001p-104, 0x1.2000000000002p+1, 457 ALL_STD_EXCEPT, FE_INEXACT); 458 459 fesetround(FE_TONEAREST); 460 test(fmaf, 0x1.800002p+0, 0x1.800002p+0, -0x1.000002p-46, 0x1.200002p+1, 461 ALL_STD_EXCEPT, FE_INEXACT); 462 fesetround(FE_DOWNWARD); 463 test(fmaf, 0x1.800002p+0, 0x1.800002p+0, -0x1.000002p-46, 0x1.200002p+1, 464 ALL_STD_EXCEPT, FE_INEXACT); 465 fesetround(FE_UPWARD); 466 test(fmaf, 0x1.800002p+0, 0x1.800002p+0, -0x1.000002p-46, 0x1.200004p+1, 467 ALL_STD_EXCEPT, FE_INEXACT); 468 469 fesetround(FE_TONEAREST); 470 #if LDBL_MANT_DIG == 64 471 test(fmal, 0x1.4p+0L, 0x1.0000000000000004p+0L, 0x1p-128L, 472 0x1.4000000000000006p+0L, ALL_STD_EXCEPT, FE_INEXACT); 473 #elif LDBL_MANT_DIG == 113 474 test(fmal, 0x1.8000000000000000000000000001p+0L, 475 0x1.8000000000000000000000000001p+0L, 476 -0x1.0000000000000000000000000001p-224L, 477 0x1.2000000000000000000000000001p+1L, ALL_STD_EXCEPT, FE_INEXACT); 478 #endif 479 480 } 481 482 int 483 main(int argc, char *argv[]) 484 { 485 int rmodes[] = { FE_TONEAREST, FE_UPWARD, FE_DOWNWARD, FE_TOWARDZERO }; 486 int i; 487 488 printf("1..19\n"); 489 490 for (i = 0; i < 4; i++) { 491 fesetround(rmodes[i]); 492 test_zeroes(); 493 printf("ok %d - fma zeroes\n", i + 1); 494 } 495 496 for (i = 0; i < 4; i++) { 497 fesetround(rmodes[i]); 498 test_infinities(); 499 printf("ok %d - fma infinities\n", i + 5); 500 } 501 502 fesetround(FE_TONEAREST); 503 test_nans(); 504 printf("ok 9 - fma NaNs\n"); 505 506 for (i = 0; i < 4; i++) { 507 fesetround(rmodes[i]); 508 test_small_z(); 509 printf("ok %d - fma small z\n", i + 10); 510 } 511 512 for (i = 0; i < 4; i++) { 513 fesetround(rmodes[i]); 514 test_big_z(); 515 printf("ok %d - fma big z\n", i + 14); 516 } 517 518 fesetround(FE_TONEAREST); 519 test_accuracy(); 520 printf("ok 18 - fma accuracy\n"); 521 522 test_double_rounding(); 523 printf("ok 19 - fma double rounding\n"); 524 525 /* 526 * TODO: 527 * - Tests for subnormals 528 * - Cancellation tests (e.g., z = (double)x*y, but x*y is inexact) 529 */ 530 531 return (0); 532 } 533