fma_test.c revision 177876
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 177876 2008-04-03 06:15:58Z 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)) || (isnan(x) && isnan(y))); 89} 90 91static void 92test_zeroes(void) 93{ 94 const int rd = (fegetround() == FE_DOWNWARD); 95 96 testall(0.0, 0.0, 0.0, 0.0, ALL_STD_EXCEPT, 0); 97 testall(1.0, 0.0, 0.0, 0.0, ALL_STD_EXCEPT, 0); 98 testall(0.0, 1.0, 0.0, 0.0, ALL_STD_EXCEPT, 0); 99 testall(0.0, 0.0, 1.0, 1.0, ALL_STD_EXCEPT, 0); 100 101 testall(-0.0, 0.0, 0.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); 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, 0.0, ALL_STD_EXCEPT, 0); 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 107 testall(-0.0, 0.0, -0.0, -0.0, ALL_STD_EXCEPT, 0); 108 testall(0.0, -0.0, -0.0, -0.0, ALL_STD_EXCEPT, 0); 109 110 testall(-1.0, 1.0, 1.0, rd ? -0.0 : 0.0, ALL_STD_EXCEPT, 0); 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 114 switch (fegetround()) { 115 case FE_TONEAREST: 116 case FE_TOWARDZERO: 117 test(fmaf, -FLT_MIN, FLT_MIN, 0.0, -0.0, 118 ALL_STD_EXCEPT, FE_INEXACT | FE_UNDERFLOW); 119 test(fma, -DBL_MIN, DBL_MIN, 0.0, -0.0, 120 ALL_STD_EXCEPT, FE_INEXACT | FE_UNDERFLOW); 121 test(fmal, -LDBL_MIN, LDBL_MIN, 0.0, -0.0, 122 ALL_STD_EXCEPT, FE_INEXACT | FE_UNDERFLOW); 123 } 124} 125 126static void 127test_infinities(void) 128{ 129 130 testall(INFINITY, 1.0, -1.0, INFINITY, ALL_STD_EXCEPT, 0); 131 testall(-1.0, INFINITY, 0.0, -INFINITY, ALL_STD_EXCEPT, 0); 132 testall(0.0, 0.0, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); 133 testall(1.0, 1.0, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); 134 testall(1.0, 1.0, -INFINITY, -INFINITY, ALL_STD_EXCEPT, 0); 135 136 testall(INFINITY, -INFINITY, 1.0, -INFINITY, ALL_STD_EXCEPT, 0); 137 testall(INFINITY, INFINITY, 1.0, INFINITY, ALL_STD_EXCEPT, 0); 138 testall(-INFINITY, -INFINITY, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); 139 140 testall(0.0, INFINITY, 1.0, NAN, ALL_STD_EXCEPT, FE_INVALID); 141 testall(INFINITY, 0.0, -0.0, NAN, ALL_STD_EXCEPT, FE_INVALID); 142 143 /* The invalid exception is optional in this case. */ 144 testall(INFINITY, 0.0, NAN, NAN, ALL_STD_EXCEPT & ~FE_INVALID, 0); 145 146 testall(INFINITY, INFINITY, -INFINITY, NAN, 147 ALL_STD_EXCEPT, FE_INVALID); 148 testall(-INFINITY, INFINITY, INFINITY, NAN, 149 ALL_STD_EXCEPT, FE_INVALID); 150 testall(INFINITY, -1.0, INFINITY, NAN, 151 ALL_STD_EXCEPT, FE_INVALID); 152 153 test(fmaf, FLT_MAX, FLT_MAX, -INFINITY, -INFINITY, ALL_STD_EXCEPT, 0); 154 test(fma, DBL_MAX, DBL_MAX, -INFINITY, -INFINITY, ALL_STD_EXCEPT, 0); 155 test(fmal, LDBL_MAX, LDBL_MAX, -INFINITY, -INFINITY, 156 ALL_STD_EXCEPT, 0); 157 test(fmaf, FLT_MAX, -FLT_MAX, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); 158 test(fma, DBL_MAX, -DBL_MAX, INFINITY, INFINITY, ALL_STD_EXCEPT, 0); 159 test(fmal, LDBL_MAX, -LDBL_MAX, INFINITY, INFINITY, 160 ALL_STD_EXCEPT, 0); 161} 162 163static void 164test_nans(void) 165{ 166 167 testall(NAN, 0.0, 0.0, NAN, ALL_STD_EXCEPT, 0); 168 testall(1.0, NAN, 1.0, NAN, ALL_STD_EXCEPT, 0); 169 testall(1.0, -1.0, NAN, NAN, ALL_STD_EXCEPT, 0); 170 testall(0.0, 0.0, NAN, NAN, ALL_STD_EXCEPT, 0); 171 testall(NAN, NAN, NAN, NAN, ALL_STD_EXCEPT, 0); 172 173 /* x*y should not raise an inexact/overflow/underflow if z is NaN. */ 174 testall(M_PI, M_PI, NAN, NAN, ALL_STD_EXCEPT, 0); 175 test(fmaf, FLT_MIN, FLT_MIN, NAN, NAN, ALL_STD_EXCEPT, 0); 176 test(fma, DBL_MIN, DBL_MIN, NAN, NAN, ALL_STD_EXCEPT, 0); 177 test(fmal, LDBL_MIN, LDBL_MIN, NAN, NAN, ALL_STD_EXCEPT, 0); 178 test(fmaf, FLT_MAX, FLT_MAX, NAN, NAN, ALL_STD_EXCEPT, 0); 179 test(fma, DBL_MAX, DBL_MAX, NAN, NAN, ALL_STD_EXCEPT, 0); 180 test(fmal, LDBL_MAX, LDBL_MAX, NAN, NAN, ALL_STD_EXCEPT, 0); 181} 182 183/* 184 * Tests for cases where z is very small compared to x*y. 185 */ 186static void 187test_small_z(void) 188{ 189 190 /* x*y positive, z positive */ 191 if (fegetround() == FE_UPWARD) { 192 test(fmaf, 1.0, 1.0, 0x1.0p-100, 1.0 + FLT_EPSILON, 193 ALL_STD_EXCEPT, FE_INEXACT); 194 test(fma, 1.0, 1.0, 0x1.0p-200, 1.0 + DBL_EPSILON, 195 ALL_STD_EXCEPT, FE_INEXACT); 196 test(fmal, 1.0, 1.0, 0x1.0p-200, 1.0 + LDBL_EPSILON, 197 ALL_STD_EXCEPT, FE_INEXACT); 198 } else { 199 testall(0x1.0p100, 1.0, 0x1.0p-100, 0x1.0p100, 200 ALL_STD_EXCEPT, FE_INEXACT); 201 } 202 203 /* x*y negative, z negative */ 204 if (fegetround() == FE_DOWNWARD) { 205 test(fmaf, -1.0, 1.0, -0x1.0p-100, -(1.0 + FLT_EPSILON), 206 ALL_STD_EXCEPT, FE_INEXACT); 207 test(fma, -1.0, 1.0, -0x1.0p-200, -(1.0 + DBL_EPSILON), 208 ALL_STD_EXCEPT, FE_INEXACT); 209 test(fmal, -1.0, 1.0, -0x1.0p-200, -(1.0 + LDBL_EPSILON), 210 ALL_STD_EXCEPT, FE_INEXACT); 211 } else { 212 testall(0x1.0p100, -1.0, -0x1.0p-100, -0x1.0p100, 213 ALL_STD_EXCEPT, FE_INEXACT); 214 } 215 216 /* x*y positive, z negative */ 217 if (fegetround() == FE_DOWNWARD || fegetround() == FE_TOWARDZERO) { 218 test(fmaf, 1.0, 1.0, -0x1.0p-100, 1.0 - FLT_EPSILON / 2, 219 ALL_STD_EXCEPT, FE_INEXACT); 220 test(fma, 1.0, 1.0, -0x1.0p-200, 1.0 - DBL_EPSILON / 2, 221 ALL_STD_EXCEPT, FE_INEXACT); 222 test(fmal, 1.0, 1.0, -0x1.0p-200, 1.0 - LDBL_EPSILON / 2, 223 ALL_STD_EXCEPT, FE_INEXACT); 224 } else { 225 testall(0x1.0p100, 1.0, -0x1.0p-100, 0x1.0p100, 226 ALL_STD_EXCEPT, FE_INEXACT); 227 } 228 229 /* x*y negative, z positive */ 230 if (fegetround() == FE_UPWARD || fegetround() == FE_TOWARDZERO) { 231 test(fmaf, -1.0, 1.0, 0x1.0p-100, -1.0 + FLT_EPSILON / 2, 232 ALL_STD_EXCEPT, FE_INEXACT); 233 test(fma, -1.0, 1.0, 0x1.0p-200, -1.0 + DBL_EPSILON / 2, 234 ALL_STD_EXCEPT, FE_INEXACT); 235 test(fmal, -1.0, 1.0, 0x1.0p-200, -1.0 + LDBL_EPSILON / 2, 236 ALL_STD_EXCEPT, FE_INEXACT); 237 } else { 238 testall(-0x1.0p100, 1.0, 0x1.0p-100, -0x1.0p100, 239 ALL_STD_EXCEPT, FE_INEXACT); 240 } 241} 242 243/* 244 * Tests for cases where z is very large compared to x*y. 245 */ 246static void 247test_big_z(void) 248{ 249 250 /* z positive, x*y positive */ 251 if (fegetround() == FE_UPWARD) { 252 test(fmaf, 0x1.0p-50, 0x1.0p-50, 1.0, 1.0 + FLT_EPSILON, 253 ALL_STD_EXCEPT, FE_INEXACT); 254 test(fma, 0x1.0p-100, 0x1.0p-100, 1.0, 1.0 + DBL_EPSILON, 255 ALL_STD_EXCEPT, FE_INEXACT); 256 test(fmal, 0x1.0p-100, 0x1.0p-100, 1.0, 1.0 + LDBL_EPSILON, 257 ALL_STD_EXCEPT, FE_INEXACT); 258 } else { 259 testall(-0x1.0p-50, -0x1.0p-50, 0x1.0p100, 0x1.0p100, 260 ALL_STD_EXCEPT, FE_INEXACT); 261 } 262 263 /* z negative, x*y negative */ 264 if (fegetround() == FE_DOWNWARD) { 265 test(fmaf, -0x1.0p-50, 0x1.0p-50, -1.0, -(1.0 + FLT_EPSILON), 266 ALL_STD_EXCEPT, FE_INEXACT); 267 test(fma, -0x1.0p-100, 0x1.0p-100, -1.0, -(1.0 + DBL_EPSILON), 268 ALL_STD_EXCEPT, FE_INEXACT); 269 test(fmal, -0x1.0p-100, 0x1.0p-100, -1.0, -(1.0 + LDBL_EPSILON), 270 ALL_STD_EXCEPT, FE_INEXACT); 271 } else { 272 testall(0x1.0p-50, -0x1.0p-50, -0x1.0p100, -0x1.0p100, 273 ALL_STD_EXCEPT, FE_INEXACT); 274 } 275 276 /* z negative, x*y positive */ 277 if (fegetround() == FE_UPWARD || fegetround() == FE_TOWARDZERO) { 278 test(fmaf, -0x1.0p-50, -0x1.0p-50, -1.0, 279 -1.0 + FLT_EPSILON / 2, ALL_STD_EXCEPT, FE_INEXACT); 280 test(fma, -0x1.0p-100, -0x1.0p-100, -1.0, 281 -1.0 + DBL_EPSILON / 2, ALL_STD_EXCEPT, FE_INEXACT); 282 test(fmal, -0x1.0p-100, -0x1.0p-100, -1.0, 283 -1.0 + LDBL_EPSILON / 2, ALL_STD_EXCEPT, FE_INEXACT); 284 } else { 285 testall(0x1.0p-50, 0x1.0p-50, -0x1.0p100, -0x1.0p100, 286 ALL_STD_EXCEPT, FE_INEXACT); 287 } 288 289 /* z positive, x*y negative */ 290 if (fegetround() == FE_DOWNWARD || fegetround() == FE_TOWARDZERO) { 291 test(fmaf, 0x1.0p-50, -0x1.0p-50, 1.0, 1.0 - FLT_EPSILON / 2, 292 ALL_STD_EXCEPT, FE_INEXACT); 293 test(fma, 0x1.0p-100, -0x1.0p-100, 1.0, 1.0 - DBL_EPSILON / 2, 294 ALL_STD_EXCEPT, FE_INEXACT); 295 test(fmal, 0x1.0p-100, -0x1.0p-100, 1.0, 1.0 - LDBL_EPSILON / 2, 296 ALL_STD_EXCEPT, FE_INEXACT); 297 } else { 298 testall(-0x1.0p-50, 0x1.0p-50, 0x1.0p100, 0x1.0p100, 299 ALL_STD_EXCEPT, FE_INEXACT); 300 } 301} 302 303static void 304test_accuracy(void) 305{ 306 307 /* ilogb(x*y) - ilogb(z) = 20 */ 308 testrnd(fmaf, -0x1.c139d8p-51, -0x1.600e7ap32, 0x1.26558cp-38, 309 0x1.34e48ap-18, 0x1.34e48cp-18, 0x1.34e48ap-18, 0x1.34e48ap-18, 310 ALL_STD_EXCEPT, FE_INEXACT); 311 testrnd(fma, -0x1.c139d7b84f1a3p-51, -0x1.600e7a2a16484p32, 312 0x1.26558cac31580p-38, 0x1.34e48a78aae97p-18, 313 0x1.34e48a78aae97p-18, 0x1.34e48a78aae96p-18, 314 0x1.34e48a78aae96p-18, ALL_STD_EXCEPT, FE_INEXACT); 315#if LDBL_MANT_DIG == 113 316 testrnd(fmal, -0x1.c139d7b84f1a3079263afcc5bae3p-51L, 317 -0x1.600e7a2a164840edbe2e7d301a72p32L, 318 0x1.26558cac315807eb07e448042101p-38L, 319 0x1.34e48a78aae96c76ed36077dd387p-18L, 320 0x1.34e48a78aae96c76ed36077dd388p-18L, 321 0x1.34e48a78aae96c76ed36077dd387p-18L, 322 0x1.34e48a78aae96c76ed36077dd387p-18L, 323 ALL_STD_EXCEPT, FE_INEXACT); 324#elif LDBL_MANT_DIG == 64 325 testrnd(fmal, -0x1.c139d7b84f1a307ap-51L, -0x1.600e7a2a164840eep32L, 326 0x1.26558cac315807ecp-38L, 0x1.34e48a78aae96c78p-18L, 327 0x1.34e48a78aae96c78p-18L, 0x1.34e48a78aae96c76p-18L, 328 0x1.34e48a78aae96c76p-18L, ALL_STD_EXCEPT, FE_INEXACT); 329#elif LDBL_MANT_DIG == 53 330 testrnd(fmal, -0x1.c139d7b84f1a3p-51L, -0x1.600e7a2a16484p32L, 331 0x1.26558cac31580p-38L, 0x1.34e48a78aae97p-18L, 332 0x1.34e48a78aae97p-18L, 0x1.34e48a78aae96p-18L, 333 0x1.34e48a78aae96p-18L, ALL_STD_EXCEPT, FE_INEXACT); 334#endif 335 336 /* ilogb(x*y) - ilogb(z) = -40 */ 337 testrnd(fmaf, 0x1.98210ap53, 0x1.9556acp-24, 0x1.d87da4p70, 338 0x1.d87da4p70, 0x1.d87da6p70, 0x1.d87da4p70, 0x1.d87da4p70, 339 ALL_STD_EXCEPT, FE_INEXACT); 340 testrnd(fma, 0x1.98210ac83fe2bp53, 0x1.9556ac1475f0fp-24, 341 0x1.d87da3aafc60ep70, 0x1.d87da3aafda40p70, 342 0x1.d87da3aafda40p70, 0x1.d87da3aafda3fp70, 343 0x1.d87da3aafda3fp70, ALL_STD_EXCEPT, FE_INEXACT); 344#if LDBL_MANT_DIG == 113 345 testrnd(fmal, 0x1.98210ac83fe2a8f65b6278b74cebp53L, 346 0x1.9556ac1475f0f28968b61d0de65ap-24L, 347 0x1.d87da3aafc60d830aa4c6d73b749p70L, 348 0x1.d87da3aafda3f36a69eb86488224p70L, 349 0x1.d87da3aafda3f36a69eb86488225p70L, 350 0x1.d87da3aafda3f36a69eb86488224p70L, 351 0x1.d87da3aafda3f36a69eb86488224p70L, 352 ALL_STD_EXCEPT, FE_INEXACT); 353#elif LDBL_MANT_DIG == 64 354 testrnd(fmal, 0x1.98210ac83fe2a8f6p53L, 0x1.9556ac1475f0f28ap-24L, 355 0x1.d87da3aafc60d83p70L, 0x1.d87da3aafda3f36ap70L, 356 0x1.d87da3aafda3f36ap70L, 0x1.d87da3aafda3f368p70L, 357 0x1.d87da3aafda3f368p70L, ALL_STD_EXCEPT, FE_INEXACT); 358#elif LDBL_MANT_DIG == 53 359 testrnd(fmal, 0x1.98210ac83fe2bp53L, 0x1.9556ac1475f0fp-24L, 360 0x1.d87da3aafc60ep70L, 0x1.d87da3aafda40p70L, 361 0x1.d87da3aafda40p70L, 0x1.d87da3aafda3fp70L, 362 0x1.d87da3aafda3fp70L, ALL_STD_EXCEPT, FE_INEXACT); 363#endif 364} 365 366int 367main(int argc, char *argv[]) 368{ 369 int rmodes[] = { FE_TONEAREST, FE_UPWARD, FE_DOWNWARD, FE_TOWARDZERO }; 370 int i; 371 372 printf("1..18\n"); 373 374 for (i = 0; i < 4; i++) { 375 fesetround(rmodes[i]); 376 test_zeroes(); 377 printf("ok %d - fma zeroes\n", i + 1); 378 } 379 380 for (i = 0; i < 4; i++) { 381 fesetround(rmodes[i]); 382 test_infinities(); 383 printf("ok %d - fma infinities\n", i + 5); 384 } 385 386 fesetround(FE_TONEAREST); 387 test_nans(); 388 printf("ok 9 - fma NaNs\n"); 389 390 for (i = 0; i < 4; i++) { 391 fesetround(rmodes[i]); 392 test_small_z(); 393 printf("ok %d - fma small z\n", i + 10); 394 } 395 396 for (i = 0; i < 4; i++) { 397 fesetround(rmodes[i]); 398 test_big_z(); 399 printf("ok %d - fma big z\n", i + 14); 400 } 401 402 fesetround(FE_TONEAREST); 403 test_accuracy(); 404 printf("ok 18 - fma accuracy\n"); 405 406 /* 407 * TODO: 408 * - Tests for subnormals 409 * - Cancellation tests (e.g., z = (double)x*y, but x*y is inexact) 410 */ 411 412 return (0); 413} 414