1 /*-
2 * Copyright (c) 2007 David Schultz <das@FreeBSD.org>
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
7 * are met:
8 * 1. Redistributions of source code must retain the above copyright
9 * notice, this list of conditions and the following disclaimer.
10 * 2. Redistributions in binary form must reproduce the above copyright
11 * notice, this list of conditions and the following disclaimer in the
12 * documentation and/or other materials provided with the distribution.
13 *
14 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
15 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
16 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
17 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
18 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
19 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
20 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
21 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
22 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
23 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
24 * SUCH DAMAGE.
25 *
26 * $FreeBSD: src/tools/regression/lib/msun/test-csqrt.c,v 1.2 2008/03/30 20:09:51 das Exp $
27 */
28
29 /*
30 * Tests for csqrt{,f}()
31 */
32
33 #include <assert.h>
34 #include <complex.h>
35 #include <float.h>
36 #include <math.h>
37 #include <stdio.h>
38
39 #define N(i) (sizeof(i) / sizeof((i)[0]))
40
41 /*
42 * This is a test hook that can point to csqrtl(), _csqrt(), or to _csqrtf().
43 * The latter two convert to float or double, respectively, and test csqrtf()
44 * and csqrt() with the same arguments.
45 */
46 long double complex (*t_csqrt)(long double complex);
47
48 static long double complex
_csqrtf(long double complex d)49 _csqrtf(long double complex d)
50 {
51
52 return (csqrtf((float complex)d));
53 }
54
55 static long double complex
_csqrt(long double complex d)56 _csqrt(long double complex d)
57 {
58
59 return (csqrt((double complex)d));
60 }
61
62 #pragma STDC CX_LIMITED_RANGE off
63
64 /*
65 * XXX gcc implements complex multiplication incorrectly. In
66 * particular, it implements it as if the CX_LIMITED_RANGE pragma
67 * were ON. Consequently, we need this function to form numbers
68 * such as x + INFINITY * I, since gcc evalutes INFINITY * I as
69 * NaN + INFINITY * I.
70 */
71 static inline long double complex
cpackl(long double x,long double y)72 cpackl(long double x, long double y)
73 {
74 long double complex z;
75
76 __real__ z = x;
77 __imag__ z = y;
78 return (z);
79 }
80
81 /*
82 * Compare d1 and d2 using special rules: NaN == NaN and +0 != -0.
83 * Fail an assertion if they differ.
84 */
85 static void
assert_equal(long double complex d1,long double complex d2)86 assert_equal(long double complex d1, long double complex d2)
87 {
88
89 if (isnan(creall(d1))) {
90 assert(isnan(creall(d2)));
91 } else {
92 assert(creall(d1) == creall(d2));
93 assert(copysignl(1.0, creall(d1)) ==
94 copysignl(1.0, creall(d2)));
95 }
96 if (isnan(cimagl(d1))) {
97 assert(isnan(cimagl(d2)));
98 } else {
99 assert(cimagl(d1) == cimagl(d2));
100 assert(copysignl(1.0, cimagl(d1)) ==
101 copysignl(1.0, cimagl(d2)));
102 }
103 }
104
105 /*
106 * Test csqrt for some finite arguments where the answer is exact.
107 * (We do not test if it produces correctly rounded answers when the
108 * result is inexact, nor do we check whether it throws spurious
109 * exceptions.)
110 */
111 static void
test_finite()112 test_finite()
113 {
114 static const double tests[] = {
115 /* csqrt(a + bI) = x + yI */
116 /* a b x y */
117 0, 8, 2, 2,
118 0, -8, 2, -2,
119 4, 0, 2, 0,
120 -4, 0, 0, 2,
121 3, 4, 2, 1,
122 3, -4, 2, -1,
123 -3, 4, 1, 2,
124 -3, -4, 1, -2,
125 5, 12, 3, 2,
126 7, 24, 4, 3,
127 9, 40, 5, 4,
128 11, 60, 6, 5,
129 13, 84, 7, 6,
130 33, 56, 7, 4,
131 39, 80, 8, 5,
132 65, 72, 9, 4,
133 987, 9916, 74, 67,
134 5289, 6640, 83, 40,
135 460766389075.0, 16762287900.0, 678910, 12345
136 };
137 /*
138 * We also test some multiples of the above arguments. This
139 * array defines which multiples we use. Note that these have
140 * to be small enough to not cause overflow for float precision
141 * with all of the constants in the above table.
142 */
143 static const double mults[] = {
144 1,
145 2,
146 3,
147 13,
148 16,
149 0x1.p30,
150 0x1.p-30,
151 };
152
153 double a, b;
154 double x, y;
155 int i, j;
156
157 for (i = 0; i < N(tests); i += 4) {
158 for (j = 0; j < N(mults); j++) {
159 a = tests[i] * mults[j] * mults[j];
160 b = tests[i + 1] * mults[j] * mults[j];
161 x = tests[i + 2] * mults[j];
162 y = tests[i + 3] * mults[j];
163 assert(t_csqrt(cpackl(a, b)) == cpackl(x, y));
164 }
165 }
166
167 }
168
169 /*
170 * Test the handling of +/- 0.
171 */
172 static void
test_zeros()173 test_zeros()
174 {
175
176 assert_equal(t_csqrt(cpackl(0.0, 0.0)), cpackl(0.0, 0.0));
177 assert_equal(t_csqrt(cpackl(-0.0, 0.0)), cpackl(0.0, 0.0));
178 assert_equal(t_csqrt(cpackl(0.0, -0.0)), cpackl(0.0, -0.0));
179 assert_equal(t_csqrt(cpackl(-0.0, -0.0)), cpackl(0.0, -0.0));
180 }
181
182 /*
183 * Test the handling of infinities when the other argument is not NaN.
184 */
185 static void
test_infinities()186 test_infinities()
187 {
188 static const double vals[] = {
189 0.0,
190 -0.0,
191 42.0,
192 -42.0,
193 INFINITY,
194 -INFINITY,
195 };
196
197 int i;
198
199 for (i = 0; i < N(vals); i++) {
200 if (isfinite(vals[i])) {
201 assert_equal(t_csqrt(cpackl(-INFINITY, vals[i])),
202 cpackl(0.0, copysignl(INFINITY, vals[i])));
203 assert_equal(t_csqrt(cpackl(INFINITY, vals[i])),
204 cpackl(INFINITY, copysignl(0.0, vals[i])));
205 }
206 assert_equal(t_csqrt(cpackl(vals[i], INFINITY)),
207 cpackl(INFINITY, INFINITY));
208 assert_equal(t_csqrt(cpackl(vals[i], -INFINITY)),
209 cpackl(INFINITY, -INFINITY));
210 }
211 }
212
213 /*
214 * Test the handling of NaNs.
215 */
216 static void
test_nans()217 test_nans()
218 {
219
220 assert(creall(t_csqrt(cpackl(INFINITY, NAN))) == INFINITY);
221 assert(isnan(cimagl(t_csqrt(cpackl(INFINITY, NAN)))));
222
223 assert(isnan(creall(t_csqrt(cpackl(-INFINITY, NAN)))));
224 assert(isinf(cimagl(t_csqrt(cpackl(-INFINITY, NAN)))));
225
226 assert_equal(t_csqrt(cpackl(NAN, INFINITY)),
227 cpackl(INFINITY, INFINITY));
228 assert_equal(t_csqrt(cpackl(NAN, -INFINITY)),
229 cpackl(INFINITY, -INFINITY));
230
231 assert_equal(t_csqrt(cpackl(0.0, NAN)), cpackl(NAN, NAN));
232 assert_equal(t_csqrt(cpackl(-0.0, NAN)), cpackl(NAN, NAN));
233 assert_equal(t_csqrt(cpackl(42.0, NAN)), cpackl(NAN, NAN));
234 assert_equal(t_csqrt(cpackl(-42.0, NAN)), cpackl(NAN, NAN));
235 assert_equal(t_csqrt(cpackl(NAN, 0.0)), cpackl(NAN, NAN));
236 assert_equal(t_csqrt(cpackl(NAN, -0.0)), cpackl(NAN, NAN));
237 assert_equal(t_csqrt(cpackl(NAN, 42.0)), cpackl(NAN, NAN));
238 assert_equal(t_csqrt(cpackl(NAN, -42.0)), cpackl(NAN, NAN));
239 assert_equal(t_csqrt(cpackl(NAN, NAN)), cpackl(NAN, NAN));
240 }
241
242 /*
243 * Test whether csqrt(a + bi) works for inputs that are large enough to
244 * cause overflow in hypot(a, b) + a. In this case we are using
245 * csqrt(115 + 252*I) == 14 + 9*I
246 * scaled up to near MAX_EXP.
247 */
248 static void
test_overflow(int maxexp)249 test_overflow(int maxexp)
250 {
251 long double a, b;
252 long double complex result;
253
254 a = ldexpl(115 * 0x1p-8, maxexp);
255 b = ldexpl(252 * 0x1p-8, maxexp);
256 result = t_csqrt(cpackl(a, b));
257 assert(creall(result) == ldexpl(14 * 0x1p-4, maxexp / 2));
258 assert(cimagl(result) == ldexpl(9 * 0x1p-4, maxexp / 2));
259 }
260
261 int
main(int argc,char * argv[])262 main(int argc, char *argv[])
263 {
264
265 printf("1..15\n");
266
267 /* Test csqrt() */
268 t_csqrt = _csqrt;
269
270 test_finite();
271 printf("ok 1 - csqrt\n");
272
273 test_zeros();
274 printf("ok 2 - csqrt\n");
275
276 test_infinities();
277 printf("ok 3 - csqrt\n");
278
279 test_nans();
280 printf("ok 4 - csqrt\n");
281
282 test_overflow(DBL_MAX_EXP);
283 printf("ok 5 - csqrt\n");
284
285 /* Now test csqrtf() */
286 t_csqrt = _csqrtf;
287
288 test_finite();
289 printf("ok 6 - csqrt\n");
290
291 test_zeros();
292 printf("ok 7 - csqrt\n");
293
294 test_infinities();
295 printf("ok 8 - csqrt\n");
296
297 test_nans();
298 printf("ok 9 - csqrt\n");
299
300 test_overflow(FLT_MAX_EXP);
301 printf("ok 10 - csqrt\n");
302
303 /* Now test csqrtl() */
304 t_csqrt = csqrtl;
305
306 test_finite();
307 printf("ok 11 - csqrt\n");
308
309 test_zeros();
310 printf("ok 12 - csqrt\n");
311
312 test_infinities();
313 printf("ok 13 - csqrt\n");
314
315 test_nans();
316 printf("ok 14 - csqrt\n");
317
318 test_overflow(LDBL_MAX_EXP);
319 printf("ok 15 - csqrt\n");
320
321 return (0);
322 }
323