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