test-fma.c revision 216222
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: head/tools/regression/lib/msun/test-fma.c 216222 2010-12-06 00:02:49Z das $"); 33 34#include <assert.h> 35#include <fenv.h> 36#include <float.h> 37#include <math.h> 38#include <stdio.h> 39 40#define ALL_STD_EXCEPT (FE_DIVBYZERO | FE_INEXACT | FE_INVALID | \ 41 FE_OVERFLOW | FE_UNDERFLOW) 42 43#pragma STDC FENV_ACCESS ON 44 45/* 46 * Test that a function returns the correct value and sets the 47 * exception flags correctly. The exceptmask specifies which 48 * exceptions we should check. We need to be lenient for several 49 * reasons, but mainly because on some architectures it's impossible 50 * to raise FE_OVERFLOW without raising FE_INEXACT. 51 * 52 * These are macros instead of functions so that assert provides more 53 * meaningful error messages. 54 */ 55#define test(func, x, y, z, result, exceptmask, excepts) do { \ 56 assert(feclearexcept(FE_ALL_EXCEPT) == 0); \ 57 assert(fpequal((func)((x), (y), (z)), (result))); \ 58 assert(((func), fetestexcept(exceptmask) == (excepts))); \ 59} while (0) 60 61#define testall(x, y, z, result, exceptmask, excepts) do { \ 62 test(fma, (x), (y), (z), (double)(result), (exceptmask), (excepts)); \ 63 test(fmaf, (x), (y), (z), (float)(result), (exceptmask), (excepts)); \ 64 test(fmal, (x), (y), (z), (result), (exceptmask), (excepts)); \ 65} while (0) 66 67/* Test in all rounding modes. */ 68#define testrnd(func, x, y, z, rn, ru, rd, rz, exceptmask, excepts) do { \ 69 fesetround(FE_TONEAREST); \ 70 test((func), (x), (y), (z), (rn), (exceptmask), (excepts)); \ 71 fesetround(FE_UPWARD); \ 72 test((func), (x), (y), (z), (ru), (exceptmask), (excepts)); \ 73 fesetround(FE_DOWNWARD); \ 74 test((func), (x), (y), (z), (rd), (exceptmask), (excepts)); \ 75 fesetround(FE_TOWARDZERO); \ 76 test((func), (x), (y), (z), (rz), (exceptmask), (excepts)); \ 77} while (0) 78 79/* 80 * Determine whether x and y are equal, with two special rules: 81 * +0.0 != -0.0 82 * NaN == NaN 83 */ 84int 85fpequal(long double x, long double y) 86{ 87 88 return ((x == y && !signbit(x) == !signbit(y)) 89 || (isnan(x) && isnan(y))); 90} 91 92static void 93test_zeroes(void) 94{ 95 const int rd = (fegetround() == FE_DOWNWARD); 96 97 testall(0.0, 0.0, 0.0, 0.0, ALL_STD_EXCEPT, 0); 98 testall(1.0, 0.0, 0.0, 0.0, ALL_STD_EXCEPT, 0); 99 testall(0.0, 1.0, 0.0, 0.0, ALL_STD_EXCEPT, 0); 100 testall(0.0, 0.0, 1.0, 1.0, ALL_STD_EXCEPT, 0); 101 102 testall(-0.0, 0.0, 0.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); 103 testall(0.0, -0.0, 0.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); 104 testall(-0.0, -0.0, 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, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); 107 108 testall(-0.0, 0.0, -0.0, -0.0, ALL_STD_EXCEPT, 0); 109 testall(0.0, -0.0, -0.0, -0.0, ALL_STD_EXCEPT, 0); 110 111 testall(-1.0, 1.0, 1.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); 112 testall(1.0, -1.0, 1.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); 113 testall(-1.0, -1.0, -1.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); 114 115 switch (fegetround()) { 116 case FE_TONEAREST: 117 case FE_TOWARDZERO: 118 test(fmaf, -FLT_MIN, FLT_MIN, 0.0, -0.0, 119 ALL_STD_EXCEPT, FE_INEXACT | FE_UNDERFLOW); 120 test(fma, -DBL_MIN, DBL_MIN, 0.0, -0.0, 121 ALL_STD_EXCEPT, FE_INEXACT | FE_UNDERFLOW); 122 test(fmal, -LDBL_MIN, LDBL_MIN, 0.0, -0.0, 123 ALL_STD_EXCEPT, FE_INEXACT | FE_UNDERFLOW); 124 } 125} 126 127static void 128test_infinities(void) 129{ 130 131 testall(INFINITY, 1.0, -1.0, INFINITY, ALL_STD_EXCEPT, 0); 132 testall(-1.0, INFINITY, 0.0, -INFINITY, ALL_STD_EXCEPT, 0); 133 testall(0.0, 0.0, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); 134 testall(1.0, 1.0, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); 135 testall(1.0, 1.0, -INFINITY, -INFINITY, ALL_STD_EXCEPT, 0); 136 137 testall(INFINITY, -INFINITY, 1.0, -INFINITY, ALL_STD_EXCEPT, 0); 138 testall(INFINITY, INFINITY, 1.0, INFINITY, ALL_STD_EXCEPT, 0); 139 testall(-INFINITY, -INFINITY, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); 140 141 testall(0.0, INFINITY, 1.0, NAN, ALL_STD_EXCEPT, FE_INVALID); 142 testall(INFINITY, 0.0, -0.0, NAN, ALL_STD_EXCEPT, FE_INVALID); 143 144 /* The invalid exception is optional in this case. */ 145 testall(INFINITY, 0.0, NAN, NAN, ALL_STD_EXCEPT & ~FE_INVALID, 0); 146 147 testall(INFINITY, INFINITY, -INFINITY, NAN, 148 ALL_STD_EXCEPT, FE_INVALID); 149 testall(-INFINITY, INFINITY, INFINITY, NAN, 150 ALL_STD_EXCEPT, FE_INVALID); 151 testall(INFINITY, -1.0, INFINITY, NAN, 152 ALL_STD_EXCEPT, FE_INVALID); 153 154 test(fmaf, FLT_MAX, FLT_MAX, -INFINITY, -INFINITY, ALL_STD_EXCEPT, 0); 155 test(fma, DBL_MAX, DBL_MAX, -INFINITY, -INFINITY, ALL_STD_EXCEPT, 0); 156 test(fmal, LDBL_MAX, LDBL_MAX, -INFINITY, -INFINITY, 157 ALL_STD_EXCEPT, 0); 158 test(fmaf, FLT_MAX, -FLT_MAX, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); 159 test(fma, DBL_MAX, -DBL_MAX, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); 160 test(fmal, LDBL_MAX, -LDBL_MAX, INFINITY, INFINITY, 161 ALL_STD_EXCEPT, 0); 162} 163 164static void 165test_nans(void) 166{ 167 168 testall(NAN, 0.0, 0.0, NAN, ALL_STD_EXCEPT, 0); 169 testall(1.0, NAN, 1.0, NAN, ALL_STD_EXCEPT, 0); 170 testall(1.0, -1.0, NAN, NAN, ALL_STD_EXCEPT, 0); 171 testall(0.0, 0.0, NAN, NAN, ALL_STD_EXCEPT, 0); 172 testall(NAN, NAN, NAN, NAN, ALL_STD_EXCEPT, 0); 173 174 /* x*y should not raise an inexact/overflow/underflow if z is NaN. */ 175 testall(M_PI, M_PI, NAN, NAN, ALL_STD_EXCEPT, 0); 176 test(fmaf, FLT_MIN, FLT_MIN, NAN, NAN, ALL_STD_EXCEPT, 0); 177 test(fma, DBL_MIN, DBL_MIN, NAN, NAN, ALL_STD_EXCEPT, 0); 178 test(fmal, LDBL_MIN, LDBL_MIN, NAN, NAN, ALL_STD_EXCEPT, 0); 179 test(fmaf, FLT_MAX, FLT_MAX, NAN, NAN, ALL_STD_EXCEPT, 0); 180 test(fma, DBL_MAX, DBL_MAX, NAN, NAN, ALL_STD_EXCEPT, 0); 181 test(fmal, LDBL_MAX, LDBL_MAX, NAN, NAN, ALL_STD_EXCEPT, 0); 182} 183 184/* 185 * Tests for cases where z is very small compared to x*y. 186 */ 187static void 188test_small_z(void) 189{ 190 191 /* x*y positive, z positive */ 192 if (fegetround() == FE_UPWARD) { 193 test(fmaf, 1.0, 1.0, 0x1.0p-100, 1.0 + FLT_EPSILON, 194 ALL_STD_EXCEPT, FE_INEXACT); 195 test(fma, 1.0, 1.0, 0x1.0p-200, 1.0 + DBL_EPSILON, 196 ALL_STD_EXCEPT, FE_INEXACT); 197 test(fmal, 1.0, 1.0, 0x1.0p-200, 1.0 + LDBL_EPSILON, 198 ALL_STD_EXCEPT, FE_INEXACT); 199 } else { 200 testall(0x1.0p100, 1.0, 0x1.0p-100, 0x1.0p100, 201 ALL_STD_EXCEPT, FE_INEXACT); 202 } 203 204 /* x*y negative, z negative */ 205 if (fegetround() == FE_DOWNWARD) { 206 test(fmaf, -1.0, 1.0, -0x1.0p-100, -(1.0 + FLT_EPSILON), 207 ALL_STD_EXCEPT, FE_INEXACT); 208 test(fma, -1.0, 1.0, -0x1.0p-200, -(1.0 + DBL_EPSILON), 209 ALL_STD_EXCEPT, FE_INEXACT); 210 test(fmal, -1.0, 1.0, -0x1.0p-200, -(1.0 + LDBL_EPSILON), 211 ALL_STD_EXCEPT, FE_INEXACT); 212 } else { 213 testall(0x1.0p100, -1.0, -0x1.0p-100, -0x1.0p100, 214 ALL_STD_EXCEPT, FE_INEXACT); 215 } 216 217 /* x*y positive, z negative */ 218 if (fegetround() == FE_DOWNWARD || fegetround() == FE_TOWARDZERO) { 219 test(fmaf, 1.0, 1.0, -0x1.0p-100, 1.0 - FLT_EPSILON / 2, 220 ALL_STD_EXCEPT, FE_INEXACT); 221 test(fma, 1.0, 1.0, -0x1.0p-200, 1.0 - DBL_EPSILON / 2, 222 ALL_STD_EXCEPT, FE_INEXACT); 223 test(fmal, 1.0, 1.0, -0x1.0p-200, 1.0 - LDBL_EPSILON / 2, 224 ALL_STD_EXCEPT, FE_INEXACT); 225 } else { 226 testall(0x1.0p100, 1.0, -0x1.0p-100, 0x1.0p100, 227 ALL_STD_EXCEPT, FE_INEXACT); 228 } 229 230 /* x*y negative, z positive */ 231 if (fegetround() == FE_UPWARD || fegetround() == FE_TOWARDZERO) { 232 test(fmaf, -1.0, 1.0, 0x1.0p-100, -1.0 + FLT_EPSILON / 2, 233 ALL_STD_EXCEPT, FE_INEXACT); 234 test(fma, -1.0, 1.0, 0x1.0p-200, -1.0 + DBL_EPSILON / 2, 235 ALL_STD_EXCEPT, FE_INEXACT); 236 test(fmal, -1.0, 1.0, 0x1.0p-200, -1.0 + LDBL_EPSILON / 2, 237 ALL_STD_EXCEPT, FE_INEXACT); 238 } else { 239 testall(-0x1.0p100, 1.0, 0x1.0p-100, -0x1.0p100, 240 ALL_STD_EXCEPT, FE_INEXACT); 241 } 242} 243 244/* 245 * Tests for cases where z is very large compared to x*y. 246 */ 247static void 248test_big_z(void) 249{ 250 251 /* z positive, x*y positive */ 252 if (fegetround() == FE_UPWARD) { 253 test(fmaf, 0x1.0p-50, 0x1.0p-50, 1.0, 1.0 + FLT_EPSILON, 254 ALL_STD_EXCEPT, FE_INEXACT); 255 test(fma, 0x1.0p-100, 0x1.0p-100, 1.0, 1.0 + DBL_EPSILON, 256 ALL_STD_EXCEPT, FE_INEXACT); 257 test(fmal, 0x1.0p-100, 0x1.0p-100, 1.0, 1.0 + LDBL_EPSILON, 258 ALL_STD_EXCEPT, FE_INEXACT); 259 } else { 260 testall(-0x1.0p-50, -0x1.0p-50, 0x1.0p100, 0x1.0p100, 261 ALL_STD_EXCEPT, FE_INEXACT); 262 } 263 264 /* z negative, x*y negative */ 265 if (fegetround() == FE_DOWNWARD) { 266 test(fmaf, -0x1.0p-50, 0x1.0p-50, -1.0, -(1.0 + FLT_EPSILON), 267 ALL_STD_EXCEPT, FE_INEXACT); 268 test(fma, -0x1.0p-100, 0x1.0p-100, -1.0, -(1.0 + DBL_EPSILON), 269 ALL_STD_EXCEPT, FE_INEXACT); 270 test(fmal, -0x1.0p-100, 0x1.0p-100, -1.0, -(1.0 + LDBL_EPSILON), 271 ALL_STD_EXCEPT, FE_INEXACT); 272 } else { 273 testall(0x1.0p-50, -0x1.0p-50, -0x1.0p100, -0x1.0p100, 274 ALL_STD_EXCEPT, FE_INEXACT); 275 } 276 277 /* z negative, x*y positive */ 278 if (fegetround() == FE_UPWARD || fegetround() == FE_TOWARDZERO) { 279 test(fmaf, -0x1.0p-50, -0x1.0p-50, -1.0, 280 -1.0 + FLT_EPSILON / 2, ALL_STD_EXCEPT, FE_INEXACT); 281 test(fma, -0x1.0p-100, -0x1.0p-100, -1.0, 282 -1.0 + DBL_EPSILON / 2, ALL_STD_EXCEPT, FE_INEXACT); 283 test(fmal, -0x1.0p-100, -0x1.0p-100, -1.0, 284 -1.0 + LDBL_EPSILON / 2, ALL_STD_EXCEPT, FE_INEXACT); 285 } else { 286 testall(0x1.0p-50, 0x1.0p-50, -0x1.0p100, -0x1.0p100, 287 ALL_STD_EXCEPT, FE_INEXACT); 288 } 289 290 /* z positive, x*y negative */ 291 if (fegetround() == FE_DOWNWARD || fegetround() == FE_TOWARDZERO) { 292 test(fmaf, 0x1.0p-50, -0x1.0p-50, 1.0, 1.0 - FLT_EPSILON / 2, 293 ALL_STD_EXCEPT, FE_INEXACT); 294 test(fma, 0x1.0p-100, -0x1.0p-100, 1.0, 1.0 - DBL_EPSILON / 2, 295 ALL_STD_EXCEPT, FE_INEXACT); 296 test(fmal, 0x1.0p-100, -0x1.0p-100, 1.0, 1.0 - LDBL_EPSILON / 2, 297 ALL_STD_EXCEPT, FE_INEXACT); 298 } else { 299 testall(-0x1.0p-50, 0x1.0p-50, 0x1.0p100, 0x1.0p100, 300 ALL_STD_EXCEPT, FE_INEXACT); 301 } 302} 303 304static void 305test_accuracy(void) 306{ 307 308 /* ilogb(x*y) - ilogb(z) = 20 */ 309 testrnd(fmaf, -0x1.c139d8p-51, -0x1.600e7ap32, 0x1.26558cp-38, 310 0x1.34e48ap-18, 0x1.34e48cp-18, 0x1.34e48ap-18, 0x1.34e48ap-18, 311 ALL_STD_EXCEPT, FE_INEXACT); 312 testrnd(fma, -0x1.c139d7b84f1a3p-51, -0x1.600e7a2a16484p32, 313 0x1.26558cac31580p-38, 0x1.34e48a78aae97p-18, 314 0x1.34e48a78aae97p-18, 0x1.34e48a78aae96p-18, 315 0x1.34e48a78aae96p-18, ALL_STD_EXCEPT, FE_INEXACT); 316#if LDBL_MANT_DIG == 113 317 testrnd(fmal, -0x1.c139d7b84f1a3079263afcc5bae3p-51L, 318 -0x1.600e7a2a164840edbe2e7d301a72p32L, 319 0x1.26558cac315807eb07e448042101p-38L, 320 0x1.34e48a78aae96c76ed36077dd387p-18L, 321 0x1.34e48a78aae96c76ed36077dd388p-18L, 322 0x1.34e48a78aae96c76ed36077dd387p-18L, 323 0x1.34e48a78aae96c76ed36077dd387p-18L, 324 ALL_STD_EXCEPT, FE_INEXACT); 325#elif LDBL_MANT_DIG == 64 326 testrnd(fmal, -0x1.c139d7b84f1a307ap-51L, -0x1.600e7a2a164840eep32L, 327 0x1.26558cac315807ecp-38L, 0x1.34e48a78aae96c78p-18L, 328 0x1.34e48a78aae96c78p-18L, 0x1.34e48a78aae96c76p-18L, 329 0x1.34e48a78aae96c76p-18L, ALL_STD_EXCEPT, FE_INEXACT); 330#elif LDBL_MANT_DIG == 53 331 testrnd(fmal, -0x1.c139d7b84f1a3p-51L, -0x1.600e7a2a16484p32L, 332 0x1.26558cac31580p-38L, 0x1.34e48a78aae97p-18L, 333 0x1.34e48a78aae97p-18L, 0x1.34e48a78aae96p-18L, 334 0x1.34e48a78aae96p-18L, ALL_STD_EXCEPT, FE_INEXACT); 335#endif 336 337 /* ilogb(x*y) - ilogb(z) = -40 */ 338 testrnd(fmaf, 0x1.98210ap53, 0x1.9556acp-24, 0x1.d87da4p70, 339 0x1.d87da4p70, 0x1.d87da6p70, 0x1.d87da4p70, 0x1.d87da4p70, 340 ALL_STD_EXCEPT, FE_INEXACT); 341 testrnd(fma, 0x1.98210ac83fe2bp53, 0x1.9556ac1475f0fp-24, 342 0x1.d87da3aafc60ep70, 0x1.d87da3aafda40p70, 343 0x1.d87da3aafda40p70, 0x1.d87da3aafda3fp70, 344 0x1.d87da3aafda3fp70, ALL_STD_EXCEPT, FE_INEXACT); 345#if LDBL_MANT_DIG == 113 346 testrnd(fmal, 0x1.98210ac83fe2a8f65b6278b74cebp53L, 347 0x1.9556ac1475f0f28968b61d0de65ap-24L, 348 0x1.d87da3aafc60d830aa4c6d73b749p70L, 349 0x1.d87da3aafda3f36a69eb86488224p70L, 350 0x1.d87da3aafda3f36a69eb86488225p70L, 351 0x1.d87da3aafda3f36a69eb86488224p70L, 352 0x1.d87da3aafda3f36a69eb86488224p70L, 353 ALL_STD_EXCEPT, FE_INEXACT); 354#elif LDBL_MANT_DIG == 64 355 testrnd(fmal, 0x1.98210ac83fe2a8f6p53L, 0x1.9556ac1475f0f28ap-24L, 356 0x1.d87da3aafc60d83p70L, 0x1.d87da3aafda3f36ap70L, 357 0x1.d87da3aafda3f36ap70L, 0x1.d87da3aafda3f368p70L, 358 0x1.d87da3aafda3f368p70L, ALL_STD_EXCEPT, FE_INEXACT); 359#elif LDBL_MANT_DIG == 53 360 testrnd(fmal, 0x1.98210ac83fe2bp53L, 0x1.9556ac1475f0fp-24L, 361 0x1.d87da3aafc60ep70L, 0x1.d87da3aafda40p70L, 362 0x1.d87da3aafda40p70L, 0x1.d87da3aafda3fp70L, 363 0x1.d87da3aafda3fp70L, ALL_STD_EXCEPT, FE_INEXACT); 364#endif 365} 366 367int 368main(int argc, char *argv[]) 369{ 370 int rmodes[] = { FE_TONEAREST, FE_UPWARD, FE_DOWNWARD, FE_TOWARDZERO }; 371 int i; 372 373 printf("1..18\n"); 374 375 for (i = 0; i < 4; i++) { 376 fesetround(rmodes[i]); 377 test_zeroes(); 378 printf("ok %d - fma zeroes\n", i + 1); 379 } 380 381 for (i = 0; i < 4; i++) { 382 fesetround(rmodes[i]); 383 test_infinities(); 384 printf("ok %d - fma infinities\n", i + 5); 385 } 386 387 fesetround(FE_TONEAREST); 388 test_nans(); 389 printf("ok 9 - fma NaNs\n"); 390 391 for (i = 0; i < 4; i++) { 392 fesetround(rmodes[i]); 393 test_small_z(); 394 printf("ok %d - fma small z\n", i + 10); 395 } 396 397 for (i = 0; i < 4; i++) { 398 fesetround(rmodes[i]); 399 test_big_z(); 400 printf("ok %d - fma big z\n", i + 14); 401 } 402 403 fesetround(FE_TONEAREST); 404 test_accuracy(); 405 printf("ok 18 - fma accuracy\n"); 406 407 /* 408 * TODO: 409 * - Tests for subnormals 410 * - Cancellation tests (e.g., z = (double)x*y, but x*y is inexact) 411 */ 412 413 return (0); 414} 415