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>
| 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/*
| 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/*
|
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
| 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);
|
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) {
| 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,
|
194 ALL_STD_EXCEPT, FE_INEXACT);
| 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,
|
196 ALL_STD_EXCEPT, FE_INEXACT);
| 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,
|
198 ALL_STD_EXCEPT, FE_INEXACT);
199 } else {
| 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,
|
201 ALL_STD_EXCEPT, FE_INEXACT);
202 }
203
204 /* x*y negative, z negative */
205 if (fegetround() == FE_DOWNWARD) {
| 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),
|
207 ALL_STD_EXCEPT, FE_INEXACT);
| 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),
|
209 ALL_STD_EXCEPT, FE_INEXACT);
| 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),
|
211 ALL_STD_EXCEPT, FE_INEXACT);
212 } else {
| 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,
|
214 ALL_STD_EXCEPT, FE_INEXACT);
215 }
216
217 /* x*y positive, z negative */
218 if (fegetround() == FE_DOWNWARD || fegetround() == FE_TOWARDZERO) {
| 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,
|
220 ALL_STD_EXCEPT, FE_INEXACT);
| 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,
|
222 ALL_STD_EXCEPT, FE_INEXACT);
| 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,
|
224 ALL_STD_EXCEPT, FE_INEXACT);
225 } else {
| 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,
|
227 ALL_STD_EXCEPT, FE_INEXACT);
228 }
229
230 /* x*y negative, z positive */
231 if (fegetround() == FE_UPWARD || fegetround() == FE_TOWARDZERO) {
| 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,
|
233 ALL_STD_EXCEPT, FE_INEXACT);
| 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,
|
235 ALL_STD_EXCEPT, FE_INEXACT);
| 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,
|
237 ALL_STD_EXCEPT, FE_INEXACT);
238 } else {
| 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,
|
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}
| 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}
|