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 __FBSDID("$FreeBSD$");
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 */
49 static long double complex (*t_csqrt)(long double complex);
50
51 static long double complex
_csqrtf(long double complex d)52 _csqrtf(long double complex d)
53 {
54
55 return (csqrtf((float complex)d));
56 }
57
58 static long double complex
_csqrt(long double complex d)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 */
71 static void
assert_equal(long double complex d1,long double complex d2)72 assert_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 */
84 static void
test_finite(void)85 test_finite(void)
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 unsigned 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 */
145 static void
test_zeros(void)146 test_zeros(void)
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 */
158 static void
test_infinities(void)159 test_infinities(void)
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 unsigned 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 */
189 static void
test_nans(void)190 test_nans(void)
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. Each of the tests is scaled up to
218 * near MAX_EXP.
219 */
220 static void
test_overflow(int maxexp)221 test_overflow(int maxexp)
222 {
223 long double a, b;
224 long double complex result;
225 int exp, i;
226
227 assert(maxexp > 0 && maxexp % 2 == 0);
228
229 for (i = 0; i < 4; i++) {
230 exp = maxexp - 2 * i;
231
232 /* csqrt(115 + 252*I) == 14 + 9*I */
233 a = ldexpl(115 * 0x1p-8, exp);
234 b = ldexpl(252 * 0x1p-8, exp);
235 result = t_csqrt(CMPLXL(a, b));
236 assert(creall(result) == ldexpl(14 * 0x1p-4, exp / 2));
237 assert(cimagl(result) == ldexpl(9 * 0x1p-4, exp / 2));
238
239 /* csqrt(-11 + 60*I) = 5 + 6*I */
240 a = ldexpl(-11 * 0x1p-6, exp);
241 b = ldexpl(60 * 0x1p-6, exp);
242 result = t_csqrt(CMPLXL(a, b));
243 assert(creall(result) == ldexpl(5 * 0x1p-3, exp / 2));
244 assert(cimagl(result) == ldexpl(6 * 0x1p-3, exp / 2));
245
246 /* csqrt(225 + 0*I) == 15 + 0*I */
247 a = ldexpl(225 * 0x1p-8, exp);
248 b = 0;
249 result = t_csqrt(CMPLXL(a, b));
250 assert(creall(result) == ldexpl(15 * 0x1p-4, exp / 2));
251 assert(cimagl(result) == 0);
252 }
253 }
254
255 /*
256 * Test that precision is maintained for some large squares. Set all or
257 * some bits in the lower mantdig/2 bits, square the number, and try to
258 * recover the sqrt. Note:
259 * (x + xI)**2 = 2xxI
260 */
261 static void
test_precision(int maxexp,int mantdig)262 test_precision(int maxexp, int mantdig)
263 {
264 long double b, x;
265 long double complex result;
266 uint64_t mantbits, sq_mantbits;
267 int exp, i;
268
269 assert(maxexp > 0 && maxexp % 2 == 0);
270 assert(mantdig <= 64);
271 mantdig = rounddown(mantdig, 2);
272
273 for (exp = 0; exp <= maxexp; exp += 2) {
274 mantbits = ((uint64_t)1 << (mantdig / 2 )) - 1;
275 for (i = 0;
276 i < 100 && mantbits > ((uint64_t)1 << (mantdig / 2 - 1));
277 i++, mantbits--) {
278 sq_mantbits = mantbits * mantbits;
279 /*
280 * sq_mantibts is a mantdig-bit number. Divide by
281 * 2**mantdig to normalize it to [0.5, 1), where,
282 * note, the binary power will be -1. Raise it by
283 * 2**exp for the test. exp is even. Lower it by
284 * one to reach a final binary power which is also
285 * even. The result should be exactly
286 * representable, given that mantdig is less than or
287 * equal to the available precision.
288 */
289 b = ldexpl((long double)sq_mantbits,
290 exp - 1 - mantdig);
291 x = ldexpl(mantbits, (exp - 2 - mantdig) / 2);
292 assert(b == x * x * 2);
293 result = t_csqrt(CMPLXL(0, b));
294 assert(creall(result) == x);
295 assert(cimagl(result) == x);
296 }
297 }
298 }
299
300 int
main(void)301 main(void)
302 {
303
304 printf("1..18\n");
305
306 /* Test csqrt() */
307 t_csqrt = _csqrt;
308
309 test_finite();
310 printf("ok 1 - csqrt\n");
311
312 test_zeros();
313 printf("ok 2 - csqrt\n");
314
315 test_infinities();
316 printf("ok 3 - csqrt\n");
317
318 test_nans();
319 printf("ok 4 - csqrt\n");
320
321 test_overflow(DBL_MAX_EXP);
322 printf("ok 5 - csqrt\n");
323
324 test_precision(DBL_MAX_EXP, DBL_MANT_DIG);
325 printf("ok 6 - csqrt\n");
326
327 /* Now test csqrtf() */
328 t_csqrt = _csqrtf;
329
330 test_finite();
331 printf("ok 7 - csqrt\n");
332
333 test_zeros();
334 printf("ok 8 - csqrt\n");
335
336 test_infinities();
337 printf("ok 9 - csqrt\n");
338
339 test_nans();
340 printf("ok 10 - csqrt\n");
341
342 test_overflow(FLT_MAX_EXP);
343 printf("ok 11 - csqrt\n");
344
345 test_precision(FLT_MAX_EXP, FLT_MANT_DIG);
346 printf("ok 12 - csqrt\n");
347
348 /* Now test csqrtl() */
349 t_csqrt = csqrtl;
350
351 test_finite();
352 printf("ok 13 - csqrt\n");
353
354 test_zeros();
355 printf("ok 14 - csqrt\n");
356
357 test_infinities();
358 printf("ok 15 - csqrt\n");
359
360 test_nans();
361 printf("ok 16 - csqrt\n");
362
363 test_overflow(LDBL_MAX_EXP);
364 printf("ok 17 - csqrt\n");
365
366 test_precision(LDBL_MAX_EXP,
367 #ifndef __i386__
368 LDBL_MANT_DIG
369 #else
370 DBL_MANT_DIG
371 #endif
372 );
373 printf("ok 18 - csqrt\n");
374
375 return (0);
376 }
377