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 $");
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 /* ilogb(x*y) - ilogb(z) = 0 */
367 testrnd(fmaf, 0x1.31ad02p+100, 0x1.2fbf7ap-42, -0x1.c3e106p+58,
368 -0x1.64c27cp+56, -0x1.64c27ap+56, -0x1.64c27cp+56,
369 -0x1.64c27ap+56, ALL_STD_EXCEPT, FE_INEXACT);
370 testrnd(fma, 0x1.31ad012ede8aap+100, 0x1.2fbf79c839067p-42,
371 -0x1.c3e106929056ep+58, -0x1.64c282b970a5fp+56,
372 -0x1.64c282b970a5ep+56, -0x1.64c282b970a5fp+56,
373 -0x1.64c282b970a5ep+56, ALL_STD_EXCEPT, FE_INEXACT);
374#if LDBL_MANT_DIG == 113
375 testrnd(fmal, 0x1.31ad012ede8aa282fa1c19376d16p+100L,
376 0x1.2fbf79c839066f0f5c68f6d2e814p-42L,
377 -0x1.c3e106929056ec19de72bfe64215p+58L,
378 -0x1.64c282b970a612598fc025ca8cddp+56L,
379 -0x1.64c282b970a612598fc025ca8cddp+56L,
380 -0x1.64c282b970a612598fc025ca8cdep+56L,
381 -0x1.64c282b970a612598fc025ca8cddp+56L,
382 ALL_STD_EXCEPT, FE_INEXACT);
383#elif LDBL_MANT_DIG == 64
384 testrnd(fmal, 0x1.31ad012ede8aa4eap+100L, 0x1.2fbf79c839066aeap-42L,
385 -0x1.c3e106929056e61p+58L, -0x1.64c282b970a60298p+56L,
386 -0x1.64c282b970a60298p+56L, -0x1.64c282b970a6029ap+56L,
387 -0x1.64c282b970a60298p+56L, ALL_STD_EXCEPT, FE_INEXACT);
388#elif LDBL_MANT_DIG == 53
389 testrnd(fmal, 0x1.31ad012ede8aap+100L, 0x1.2fbf79c839067p-42L,
390 -0x1.c3e106929056ep+58L, -0x1.64c282b970a5fp+56L,
391 -0x1.64c282b970a5ep+56L, -0x1.64c282b970a5fp+56L,
392 -0x1.64c282b970a5ep+56L, ALL_STD_EXCEPT, FE_INEXACT);
393#endif
394
395 /* x*y (rounded) ~= -z */
396 /* XXX spurious inexact exceptions */
397 testrnd(fmaf, 0x1.bbffeep-30, -0x1.1d164cp-74, 0x1.ee7296p-104,
398 -0x1.c46ea8p-128, -0x1.c46ea8p-128, -0x1.c46ea8p-128,
399 -0x1.c46ea8p-128, ALL_STD_EXCEPT & ~FE_INEXACT, 0);
400 testrnd(fma, 0x1.bbffeea6fc7d6p-30, 0x1.1d164c6cbf078p-74,
401 -0x1.ee72993aff948p-104, -0x1.71f72ac7d9d8p-159,
402 -0x1.71f72ac7d9d8p-159, -0x1.71f72ac7d9d8p-159,
403 -0x1.71f72ac7d9d8p-159, ALL_STD_EXCEPT & ~FE_INEXACT, 0);
404#if LDBL_MANT_DIG == 113
405 testrnd(fmal, 0x1.bbffeea6fc7d65927d147f437675p-30L,
406 0x1.1d164c6cbf078b7a22607d1cd6a2p-74L,
407 -0x1.ee72993aff94973876031bec0944p-104L,
408 0x1.64e086175b3a2adc36e607058814p-217L,
409 0x1.64e086175b3a2adc36e607058814p-217L,
410 0x1.64e086175b3a2adc36e607058814p-217L,
411 0x1.64e086175b3a2adc36e607058814p-217L,
412 ALL_STD_EXCEPT & ~FE_INEXACT, 0);
413#elif LDBL_MANT_DIG == 64
414 testrnd(fmal, 0x1.bbffeea6fc7d6592p-30L, 0x1.1d164c6cbf078b7ap-74L,
415 -0x1.ee72993aff949736p-104L, 0x1.af190e7a1ee6ad94p-168L,
416 0x1.af190e7a1ee6ad94p-168L, 0x1.af190e7a1ee6ad94p-168L,
417 0x1.af190e7a1ee6ad94p-168L, ALL_STD_EXCEPT & ~FE_INEXACT, 0);
418#elif LDBL_MANT_DIG == 53
419 testrnd(fmal, 0x1.bbffeea6fc7d6p-30L, 0x1.1d164c6cbf078p-74L,
420 -0x1.ee72993aff948p-104L, -0x1.71f72ac7d9d8p-159L,
421 -0x1.71f72ac7d9d8p-159L, -0x1.71f72ac7d9d8p-159L,
422 -0x1.71f72ac7d9d8p-159L, ALL_STD_EXCEPT & ~FE_INEXACT, 0);
423#endif
424}
425
426static void
427test_double_rounding(void)
428{
429
430 /*
431 * a = 0x1.8000000000001p0
432 * b = 0x1.8000000000001p0
433 * c = -0x0.0000000000000000000000000080...1p+1
434 * a * b = 0x1.2000000000001800000000000080p+1
435 *
436 * The correct behavior is to round DOWN to 0x1.2000000000001p+1 in
437 * round-to-nearest mode. An implementation that computes a*b+c in
438 * double+double precision, however, will get 0x1.20000000000018p+1,
439 * and then round UP.
440 */
441 fesetround(FE_TONEAREST);
442 test(fma, 0x1.8000000000001p0, 0x1.8000000000001p0,
443 -0x1.0000000000001p-104, 0x1.2000000000001p+1,
444 ALL_STD_EXCEPT, FE_INEXACT);
445 fesetround(FE_DOWNWARD);
446 test(fma, 0x1.8000000000001p0, 0x1.8000000000001p0,
447 -0x1.0000000000001p-104, 0x1.2000000000001p+1,
448 ALL_STD_EXCEPT, FE_INEXACT);
449 fesetround(FE_UPWARD);
450 test(fma, 0x1.8000000000001p0, 0x1.8000000000001p0,
451 -0x1.0000000000001p-104, 0x1.2000000000002p+1,
452 ALL_STD_EXCEPT, FE_INEXACT);
453
454 fesetround(FE_TONEAREST);
455 test(fmaf, 0x1.800002p+0, 0x1.800002p+0, -0x1.000002p-46, 0x1.200002p+1,
456 ALL_STD_EXCEPT, FE_INEXACT);
457 fesetround(FE_DOWNWARD);
458 test(fmaf, 0x1.800002p+0, 0x1.800002p+0, -0x1.000002p-46, 0x1.200002p+1,
459 ALL_STD_EXCEPT, FE_INEXACT);
460 fesetround(FE_UPWARD);
461 test(fmaf, 0x1.800002p+0, 0x1.800002p+0, -0x1.000002p-46, 0x1.200004p+1,
462 ALL_STD_EXCEPT, FE_INEXACT);
463
464 fesetround(FE_TONEAREST);
465#if LDBL_MANT_DIG == 64
466 test(fmal, 0x1.4p+0L, 0x1.0000000000000004p+0L, 0x1p-128L,
467 0x1.4000000000000006p+0L, ALL_STD_EXCEPT, FE_INEXACT);
468#elif LDBL_MANT_DIG == 113
469 test(fmal, 0x1.8000000000000000000000000001p+0L,
470 0x1.8000000000000000000000000001p+0L,
471 -0x1.0000000000000000000000000001p-224L,
472 0x1.2000000000000000000000000001p+1L, ALL_STD_EXCEPT, FE_INEXACT);
473#endif
474
475}
476
477int
478main(int argc, char *argv[])
479{
480 int rmodes[] = { FE_TONEAREST, FE_UPWARD, FE_DOWNWARD, FE_TOWARDZERO };
481 int i;
482
483 printf("1..19\n");
484
485 for (i = 0; i < 4; i++) {
486 fesetround(rmodes[i]);
487 test_zeroes();
488 printf("ok %d - fma zeroes\n", i + 1);
489 }
490
491 for (i = 0; i < 4; i++) {
492 fesetround(rmodes[i]);
493 test_infinities();
494 printf("ok %d - fma infinities\n", i + 5);
495 }
496
497 fesetround(FE_TONEAREST);
498 test_nans();
499 printf("ok 9 - fma NaNs\n");
500
501 for (i = 0; i < 4; i++) {
502 fesetround(rmodes[i]);
503 test_small_z();
504 printf("ok %d - fma small z\n", i + 10);
505 }
506
507 for (i = 0; i < 4; i++) {
508 fesetround(rmodes[i]);
509 test_big_z();
510 printf("ok %d - fma big z\n", i + 14);
511 }
512
513 fesetround(FE_TONEAREST);
514 test_accuracy();
515 printf("ok 18 - fma accuracy\n");
516
517 test_double_rounding();
518 printf("ok 19 - fma double rounding\n");
519
520 /*
521 * TODO:
522 * - Tests for subnormals
523 * - Cancellation tests (e.g., z = (double)x*y, but x*y is inexact)
524 */
525
526 return (0);
527}
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 $");
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 /* ilogb(x*y) - ilogb(z) = 0 */
367 testrnd(fmaf, 0x1.31ad02p+100, 0x1.2fbf7ap-42, -0x1.c3e106p+58,
368 -0x1.64c27cp+56, -0x1.64c27ap+56, -0x1.64c27cp+56,
369 -0x1.64c27ap+56, ALL_STD_EXCEPT, FE_INEXACT);
370 testrnd(fma, 0x1.31ad012ede8aap+100, 0x1.2fbf79c839067p-42,
371 -0x1.c3e106929056ep+58, -0x1.64c282b970a5fp+56,
372 -0x1.64c282b970a5ep+56, -0x1.64c282b970a5fp+56,
373 -0x1.64c282b970a5ep+56, ALL_STD_EXCEPT, FE_INEXACT);
374#if LDBL_MANT_DIG == 113
375 testrnd(fmal, 0x1.31ad012ede8aa282fa1c19376d16p+100L,
376 0x1.2fbf79c839066f0f5c68f6d2e814p-42L,
377 -0x1.c3e106929056ec19de72bfe64215p+58L,
378 -0x1.64c282b970a612598fc025ca8cddp+56L,
379 -0x1.64c282b970a612598fc025ca8cddp+56L,
380 -0x1.64c282b970a612598fc025ca8cdep+56L,
381 -0x1.64c282b970a612598fc025ca8cddp+56L,
382 ALL_STD_EXCEPT, FE_INEXACT);
383#elif LDBL_MANT_DIG == 64
384 testrnd(fmal, 0x1.31ad012ede8aa4eap+100L, 0x1.2fbf79c839066aeap-42L,
385 -0x1.c3e106929056e61p+58L, -0x1.64c282b970a60298p+56L,
386 -0x1.64c282b970a60298p+56L, -0x1.64c282b970a6029ap+56L,
387 -0x1.64c282b970a60298p+56L, ALL_STD_EXCEPT, FE_INEXACT);
388#elif LDBL_MANT_DIG == 53
389 testrnd(fmal, 0x1.31ad012ede8aap+100L, 0x1.2fbf79c839067p-42L,
390 -0x1.c3e106929056ep+58L, -0x1.64c282b970a5fp+56L,
391 -0x1.64c282b970a5ep+56L, -0x1.64c282b970a5fp+56L,
392 -0x1.64c282b970a5ep+56L, ALL_STD_EXCEPT, FE_INEXACT);
393#endif
394
395 /* x*y (rounded) ~= -z */
396 /* XXX spurious inexact exceptions */
397 testrnd(fmaf, 0x1.bbffeep-30, -0x1.1d164cp-74, 0x1.ee7296p-104,
398 -0x1.c46ea8p-128, -0x1.c46ea8p-128, -0x1.c46ea8p-128,
399 -0x1.c46ea8p-128, ALL_STD_EXCEPT & ~FE_INEXACT, 0);
400 testrnd(fma, 0x1.bbffeea6fc7d6p-30, 0x1.1d164c6cbf078p-74,
401 -0x1.ee72993aff948p-104, -0x1.71f72ac7d9d8p-159,
402 -0x1.71f72ac7d9d8p-159, -0x1.71f72ac7d9d8p-159,
403 -0x1.71f72ac7d9d8p-159, ALL_STD_EXCEPT & ~FE_INEXACT, 0);
404#if LDBL_MANT_DIG == 113
405 testrnd(fmal, 0x1.bbffeea6fc7d65927d147f437675p-30L,
406 0x1.1d164c6cbf078b7a22607d1cd6a2p-74L,
407 -0x1.ee72993aff94973876031bec0944p-104L,
408 0x1.64e086175b3a2adc36e607058814p-217L,
409 0x1.64e086175b3a2adc36e607058814p-217L,
410 0x1.64e086175b3a2adc36e607058814p-217L,
411 0x1.64e086175b3a2adc36e607058814p-217L,
412 ALL_STD_EXCEPT & ~FE_INEXACT, 0);
413#elif LDBL_MANT_DIG == 64
414 testrnd(fmal, 0x1.bbffeea6fc7d6592p-30L, 0x1.1d164c6cbf078b7ap-74L,
415 -0x1.ee72993aff949736p-104L, 0x1.af190e7a1ee6ad94p-168L,
416 0x1.af190e7a1ee6ad94p-168L, 0x1.af190e7a1ee6ad94p-168L,
417 0x1.af190e7a1ee6ad94p-168L, ALL_STD_EXCEPT & ~FE_INEXACT, 0);
418#elif LDBL_MANT_DIG == 53
419 testrnd(fmal, 0x1.bbffeea6fc7d6p-30L, 0x1.1d164c6cbf078p-74L,
420 -0x1.ee72993aff948p-104L, -0x1.71f72ac7d9d8p-159L,
421 -0x1.71f72ac7d9d8p-159L, -0x1.71f72ac7d9d8p-159L,
422 -0x1.71f72ac7d9d8p-159L, ALL_STD_EXCEPT & ~FE_INEXACT, 0);
423#endif
424}
425
426static void
427test_double_rounding(void)
428{
429
430 /*
431 * a = 0x1.8000000000001p0
432 * b = 0x1.8000000000001p0
433 * c = -0x0.0000000000000000000000000080...1p+1
434 * a * b = 0x1.2000000000001800000000000080p+1
435 *
436 * The correct behavior is to round DOWN to 0x1.2000000000001p+1 in
437 * round-to-nearest mode. An implementation that computes a*b+c in
438 * double+double precision, however, will get 0x1.20000000000018p+1,
439 * and then round UP.
440 */
441 fesetround(FE_TONEAREST);
442 test(fma, 0x1.8000000000001p0, 0x1.8000000000001p0,
443 -0x1.0000000000001p-104, 0x1.2000000000001p+1,
444 ALL_STD_EXCEPT, FE_INEXACT);
445 fesetround(FE_DOWNWARD);
446 test(fma, 0x1.8000000000001p0, 0x1.8000000000001p0,
447 -0x1.0000000000001p-104, 0x1.2000000000001p+1,
448 ALL_STD_EXCEPT, FE_INEXACT);
449 fesetround(FE_UPWARD);
450 test(fma, 0x1.8000000000001p0, 0x1.8000000000001p0,
451 -0x1.0000000000001p-104, 0x1.2000000000002p+1,
452 ALL_STD_EXCEPT, FE_INEXACT);
453
454 fesetround(FE_TONEAREST);
455 test(fmaf, 0x1.800002p+0, 0x1.800002p+0, -0x1.000002p-46, 0x1.200002p+1,
456 ALL_STD_EXCEPT, FE_INEXACT);
457 fesetround(FE_DOWNWARD);
458 test(fmaf, 0x1.800002p+0, 0x1.800002p+0, -0x1.000002p-46, 0x1.200002p+1,
459 ALL_STD_EXCEPT, FE_INEXACT);
460 fesetround(FE_UPWARD);
461 test(fmaf, 0x1.800002p+0, 0x1.800002p+0, -0x1.000002p-46, 0x1.200004p+1,
462 ALL_STD_EXCEPT, FE_INEXACT);
463
464 fesetround(FE_TONEAREST);
465#if LDBL_MANT_DIG == 64
466 test(fmal, 0x1.4p+0L, 0x1.0000000000000004p+0L, 0x1p-128L,
467 0x1.4000000000000006p+0L, ALL_STD_EXCEPT, FE_INEXACT);
468#elif LDBL_MANT_DIG == 113
469 test(fmal, 0x1.8000000000000000000000000001p+0L,
470 0x1.8000000000000000000000000001p+0L,
471 -0x1.0000000000000000000000000001p-224L,
472 0x1.2000000000000000000000000001p+1L, ALL_STD_EXCEPT, FE_INEXACT);
473#endif
474
475}
476
477int
478main(int argc, char *argv[])
479{
480 int rmodes[] = { FE_TONEAREST, FE_UPWARD, FE_DOWNWARD, FE_TOWARDZERO };
481 int i;
482
483 printf("1..19\n");
484
485 for (i = 0; i < 4; i++) {
486 fesetround(rmodes[i]);
487 test_zeroes();
488 printf("ok %d - fma zeroes\n", i + 1);
489 }
490
491 for (i = 0; i < 4; i++) {
492 fesetround(rmodes[i]);
493 test_infinities();
494 printf("ok %d - fma infinities\n", i + 5);
495 }
496
497 fesetround(FE_TONEAREST);
498 test_nans();
499 printf("ok 9 - fma NaNs\n");
500
501 for (i = 0; i < 4; i++) {
502 fesetround(rmodes[i]);
503 test_small_z();
504 printf("ok %d - fma small z\n", i + 10);
505 }
506
507 for (i = 0; i < 4; i++) {
508 fesetround(rmodes[i]);
509 test_big_z();
510 printf("ok %d - fma big z\n", i + 14);
511 }
512
513 fesetround(FE_TONEAREST);
514 test_accuracy();
515 printf("ok 18 - fma accuracy\n");
516
517 test_double_rounding();
518 printf("ok 19 - fma double rounding\n");
519
520 /*
521 * TODO:
522 * - Tests for subnormals
523 * - Cancellation tests (e.g., z = (double)x*y, but x*y is inexact)
524 */
525
526 return (0);
527}