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