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