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