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
fpequal(long double x,long double y)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
test_zeroes(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
test_infinities(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
test_nans(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
test_small_z(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
test_big_z(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
test_accuracy(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
test_double_rounding(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
main(int argc,char * argv[])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