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