xref: /freebsd/lib/msun/tests/csqrt_test.c (revision 9768746b)
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$");
33 
34 #include <sys/param.h>
35 
36 #include <complex.h>
37 #include <float.h>
38 #include <math.h>
39 #include <stdio.h>
40 
41 #include "test-utils.h"
42 
43 /*
44  * This is a test hook that can point to csqrtl(), _csqrt(), or to _csqrtf().
45  * The latter two convert to float or double, respectively, and test csqrtf()
46  * and csqrt() with the same arguments.
47  */
48 static long double complex (*t_csqrt)(long double complex);
49 
50 static long double complex
51 _csqrtf(long double complex d)
52 {
53 
54 	return (csqrtf((float complex)d));
55 }
56 
57 static long double complex
58 _csqrt(long double complex d)
59 {
60 
61 	return (csqrt((double complex)d));
62 }
63 
64 #pragma	STDC CX_LIMITED_RANGE	OFF
65 
66 /*
67  * Compare d1 and d2 using special rules: NaN == NaN and +0 != -0.
68  * Fail an assertion if they differ.
69  */
70 #define assert_equal(d1, d2) CHECK_CFPEQUAL_CS(d1, d2, CS_BOTH)
71 
72 /*
73  * Test csqrt for some finite arguments where the answer is exact.
74  * (We do not test if it produces correctly rounded answers when the
75  * result is inexact, nor do we check whether it throws spurious
76  * exceptions.)
77  */
78 static void
79 test_finite(void)
80 {
81 	static const double tests[] = {
82 	     /* csqrt(a + bI) = x + yI */
83 	     /* a	b	x	y */
84 		0,	8,	2,	2,
85 		0,	-8,	2,	-2,
86 		4,	0,	2,	0,
87 		-4,	0,	0,	2,
88 		3,	4,	2,	1,
89 		3,	-4,	2,	-1,
90 		-3,	4,	1,	2,
91 		-3,	-4,	1,	-2,
92 		5,	12,	3,	2,
93 		7,	24,	4,	3,
94 		9,	40,	5,	4,
95 		11,	60,	6,	5,
96 		13,	84,	7,	6,
97 		33,	56,	7,	4,
98 		39,	80,	8,	5,
99 		65,	72,	9,	4,
100 		987,	9916,	74,	67,
101 		5289,	6640,	83,	40,
102 		460766389075.0, 16762287900.0, 678910, 12345
103 	};
104 	/*
105 	 * We also test some multiples of the above arguments. This
106 	 * array defines which multiples we use. Note that these have
107 	 * to be small enough to not cause overflow for float precision
108 	 * with all of the constants in the above table.
109 	 */
110 	static const double mults[] = {
111 		1,
112 		2,
113 		3,
114 		13,
115 		16,
116 		0x1.p30,
117 		0x1.p-30,
118 	};
119 
120 	double a, b;
121 	double x, y;
122 	unsigned i, j;
123 
124 	for (i = 0; i < nitems(tests); i += 4) {
125 		for (j = 0; j < nitems(mults); j++) {
126 			a = tests[i] * mults[j] * mults[j];
127 			b = tests[i + 1] * mults[j] * mults[j];
128 			x = tests[i + 2] * mults[j];
129 			y = tests[i + 3] * mults[j];
130 			ATF_CHECK(t_csqrt(CMPLXL(a, b)) == CMPLXL(x, y));
131 		}
132 	}
133 
134 }
135 
136 /*
137  * Test the handling of +/- 0.
138  */
139 static void
140 test_zeros(void)
141 {
142 
143 	assert_equal(t_csqrt(CMPLXL(0.0, 0.0)), CMPLXL(0.0, 0.0));
144 	assert_equal(t_csqrt(CMPLXL(-0.0, 0.0)), CMPLXL(0.0, 0.0));
145 	assert_equal(t_csqrt(CMPLXL(0.0, -0.0)), CMPLXL(0.0, -0.0));
146 	assert_equal(t_csqrt(CMPLXL(-0.0, -0.0)), CMPLXL(0.0, -0.0));
147 }
148 
149 /*
150  * Test the handling of infinities when the other argument is not NaN.
151  */
152 static void
153 test_infinities(void)
154 {
155 	static const double vals[] = {
156 		0.0,
157 		-0.0,
158 		42.0,
159 		-42.0,
160 		INFINITY,
161 		-INFINITY,
162 	};
163 
164 	unsigned i;
165 
166 	for (i = 0; i < nitems(vals); i++) {
167 		if (isfinite(vals[i])) {
168 			assert_equal(t_csqrt(CMPLXL(-INFINITY, vals[i])),
169 			    CMPLXL(0.0, copysignl(INFINITY, vals[i])));
170 			assert_equal(t_csqrt(CMPLXL(INFINITY, vals[i])),
171 			    CMPLXL(INFINITY, copysignl(0.0, vals[i])));
172 		}
173 		assert_equal(t_csqrt(CMPLXL(vals[i], INFINITY)),
174 		    CMPLXL(INFINITY, INFINITY));
175 		assert_equal(t_csqrt(CMPLXL(vals[i], -INFINITY)),
176 		    CMPLXL(INFINITY, -INFINITY));
177 	}
178 }
179 
180 /*
181  * Test the handling of NaNs.
182  */
183 static void
184 test_nans(void)
185 {
186 
187 	ATF_CHECK(creall(t_csqrt(CMPLXL(INFINITY, NAN))) == INFINITY);
188 	ATF_CHECK(isnan(cimagl(t_csqrt(CMPLXL(INFINITY, NAN)))));
189 
190 	ATF_CHECK(isnan(creall(t_csqrt(CMPLXL(-INFINITY, NAN)))));
191 	ATF_CHECK(isinf(cimagl(t_csqrt(CMPLXL(-INFINITY, NAN)))));
192 
193 	assert_equal(t_csqrt(CMPLXL(NAN, INFINITY)),
194 		     CMPLXL(INFINITY, INFINITY));
195 	assert_equal(t_csqrt(CMPLXL(NAN, -INFINITY)),
196 		     CMPLXL(INFINITY, -INFINITY));
197 
198 	assert_equal(t_csqrt(CMPLXL(0.0, NAN)), CMPLXL(NAN, NAN));
199 	assert_equal(t_csqrt(CMPLXL(-0.0, NAN)), CMPLXL(NAN, NAN));
200 	assert_equal(t_csqrt(CMPLXL(42.0, NAN)), CMPLXL(NAN, NAN));
201 	assert_equal(t_csqrt(CMPLXL(-42.0, NAN)), CMPLXL(NAN, NAN));
202 	assert_equal(t_csqrt(CMPLXL(NAN, 0.0)), CMPLXL(NAN, NAN));
203 	assert_equal(t_csqrt(CMPLXL(NAN, -0.0)), CMPLXL(NAN, NAN));
204 	assert_equal(t_csqrt(CMPLXL(NAN, 42.0)), CMPLXL(NAN, NAN));
205 	assert_equal(t_csqrt(CMPLXL(NAN, -42.0)), CMPLXL(NAN, NAN));
206 	assert_equal(t_csqrt(CMPLXL(NAN, NAN)), CMPLXL(NAN, NAN));
207 }
208 
209 /*
210  * Test whether csqrt(a + bi) works for inputs that are large enough to
211  * cause overflow in hypot(a, b) + a.  Each of the tests is scaled up to
212  * near MAX_EXP.
213  */
214 static void
215 test_overflow(int maxexp)
216 {
217 	long double a, b;
218 	long double complex result;
219 	int exp, i;
220 
221 	ATF_CHECK(maxexp > 0 && maxexp % 2 == 0);
222 
223 	for (i = 0; i < 4; i++) {
224 		exp = maxexp - 2 * i;
225 
226 		/* csqrt(115 + 252*I) == 14 + 9*I */
227 		a = ldexpl(115 * 0x1p-8, exp);
228 		b = ldexpl(252 * 0x1p-8, exp);
229 		result = t_csqrt(CMPLXL(a, b));
230 		ATF_CHECK_EQ(creall(result), ldexpl(14 * 0x1p-4, exp / 2));
231 		ATF_CHECK_EQ(cimagl(result), ldexpl(9 * 0x1p-4, exp / 2));
232 
233 		/* csqrt(-11 + 60*I) = 5 + 6*I */
234 		a = ldexpl(-11 * 0x1p-6, exp);
235 		b = ldexpl(60 * 0x1p-6, exp);
236 		result = t_csqrt(CMPLXL(a, b));
237 		ATF_CHECK_EQ(creall(result), ldexpl(5 * 0x1p-3, exp / 2));
238 		ATF_CHECK_EQ(cimagl(result), ldexpl(6 * 0x1p-3, exp / 2));
239 
240 		/* csqrt(225 + 0*I) == 15 + 0*I */
241 		a = ldexpl(225 * 0x1p-8, exp);
242 		b = 0;
243 		result = t_csqrt(CMPLXL(a, b));
244 		ATF_CHECK_EQ(creall(result), ldexpl(15 * 0x1p-4, exp / 2));
245 		ATF_CHECK_EQ(cimagl(result), 0);
246 	}
247 }
248 
249 /*
250  * Test that precision is maintained for some large squares.  Set all or
251  * some bits in the lower mantdig/2 bits, square the number, and try to
252  * recover the sqrt.  Note:
253  * 	(x + xI)**2 = 2xxI
254  */
255 static void
256 test_precision(int maxexp, int mantdig)
257 {
258 	long double b, x;
259 	long double complex result;
260 #if LDBL_MANT_DIG <= 64
261 	typedef uint64_t ldbl_mant_type;
262 #elif LDBL_MANT_DIG <= 128
263 	typedef __uint128_t ldbl_mant_type;
264 #else
265 #error "Unsupported long double format"
266 #endif
267 	ldbl_mant_type mantbits, sq_mantbits;
268 	int exp, i;
269 
270 	ATF_REQUIRE(maxexp > 0 && maxexp % 2 == 0);
271 	ATF_REQUIRE(mantdig <= LDBL_MANT_DIG);
272 	mantdig = rounddown(mantdig, 2);
273 
274 	for (exp = 0; exp <= maxexp; exp += 2) {
275 		mantbits = ((ldbl_mant_type)1 << (mantdig / 2)) - 1;
276 		for (i = 0; i < 100 &&
277 		     mantbits > ((ldbl_mant_type)1 << (mantdig / 2 - 1));
278 		     i++, mantbits--) {
279 			sq_mantbits = mantbits * mantbits;
280 			/*
281 			 * sq_mantibts is a mantdig-bit number.  Divide by
282 			 * 2**mantdig to normalize it to [0.5, 1), where,
283 			 * note, the binary power will be -1.  Raise it by
284 			 * 2**exp for the test.  exp is even.  Lower it by
285 			 * one to reach a final binary power which is also
286 			 * even.  The result should be exactly
287 			 * representable, given that mantdig is less than or
288 			 * equal to the available precision.
289 			 */
290 			b = ldexpl((long double)sq_mantbits,
291 			    exp - 1 - mantdig);
292 			x = ldexpl(mantbits, (exp - 2 - mantdig) / 2);
293 			CHECK_FPEQUAL(b, x * x * 2);
294 			result = t_csqrt(CMPLXL(0, b));
295 			CHECK_FPEQUAL(x, creall(result));
296 			CHECK_FPEQUAL(x, cimagl(result));
297 		}
298 	}
299 }
300 
301 ATF_TC_WITHOUT_HEAD(csqrt);
302 ATF_TC_BODY(csqrt, tc)
303 {
304 	/* Test csqrt() */
305 	t_csqrt = _csqrt;
306 
307 	test_finite();
308 
309 	test_zeros();
310 
311 	test_infinities();
312 
313 	test_nans();
314 
315 	test_overflow(DBL_MAX_EXP);
316 
317 	test_precision(DBL_MAX_EXP, DBL_MANT_DIG);
318 }
319 
320 ATF_TC_WITHOUT_HEAD(csqrtf);
321 ATF_TC_BODY(csqrtf, tc)
322 {
323 	/* Now test csqrtf() */
324 	t_csqrt = _csqrtf;
325 
326 	test_finite();
327 
328 	test_zeros();
329 
330 	test_infinities();
331 
332 	test_nans();
333 
334 	test_overflow(FLT_MAX_EXP);
335 
336 	test_precision(FLT_MAX_EXP, FLT_MANT_DIG);
337 }
338 
339 ATF_TC_WITHOUT_HEAD(csqrtl);
340 ATF_TC_BODY(csqrtl, tc)
341 {
342 	/* Now test csqrtl() */
343 	t_csqrt = csqrtl;
344 
345 	test_finite();
346 
347 	test_zeros();
348 
349 	test_infinities();
350 
351 	test_nans();
352 
353 	test_overflow(LDBL_MAX_EXP);
354 
355 	/* i386 is configured to use 53-bit rounding precision for long double. */
356 	test_precision(LDBL_MAX_EXP,
357 #ifndef __i386__
358 	    LDBL_MANT_DIG
359 #else
360 	    DBL_MANT_DIG
361 #endif
362 	    );
363 }
364 
365 ATF_TP_ADD_TCS(tp)
366 {
367 	ATF_TP_ADD_TC(tp, csqrt);
368 	ATF_TP_ADD_TC(tp, csqrtf);
369 	ATF_TP_ADD_TC(tp, csqrtl);
370 
371 	return (atf_no_error());
372 }
373