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