Index: lib/libm/src/s_csqrt.c =================================================================== RCS file: /cvs/src/lib/libm/src/s_csqrt.c,v retrieving revision 1.8 diff -u -p -r1.8 s_csqrt.c --- lib/libm/src/s_csqrt.c 12 Sep 2016 19:47:02 -0000 1.8 +++ lib/libm/src/s_csqrt.c 1 Jun 2019 09:22:46 -0000 @@ -1,134 +1,116 @@ -/* $OpenBSD: s_csqrt.c,v 1.8 2016/09/12 19:47:02 guenther Exp $ */ -/* - * Copyright (c) 2008 Stephen L. Moshier +/*- + * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * - * Permission to use, copy, modify, and distribute this software for any - * purpose with or without fee is hereby granted, provided that the above - * copyright notice and this permission notice appear in all copies. + * Copyright (c) 2007 David Schultz + * All rights reserved. * - * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES - * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF - * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR - * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES - * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN - * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF - * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. */ -/* csqrt() - * - * Complex square root - * - * - * - * SYNOPSIS: - * - * double complex csqrt(); - * double complex z, w; - * - * w = csqrt (z); - * - * - * - * DESCRIPTION: - * - * - * If z = x + iy, r = |z|, then - * - * 1/2 - * Re w = [ (r + x)/2 ] , - * - * 1/2 - * Im w = [ (r - x)/2 ] . - * - * Cancellation error in r-x or r+x is avoided by using the - * identity 2 Re w Im w = y. - * - * Note that -w is also a square root of z. The root chosen - * is always in the right half plane and Im w has the same sign as y. - * - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * DEC -10,+10 25000 3.2e-17 9.6e-18 - * IEEE -10,+10 1,000,000 2.9e-16 6.1e-17 - * - */ +#include #include #include #include +#include "math_private.h" + +/* For avoiding overflow for components >= DBL_MAX / (1 + sqrt(2)). */ +#define THRESH 0x1.a827999fcef32p+1022 + double complex csqrt(double complex z) { - double complex w; - double x, y, r, t, scale; + double complex result; + double a, b, rx, ry, scale, t; - x = creal (z); - y = cimag (z); + a = creal(z); + b = cimag(z); - if (y == 0.0) { - if (x == 0.0) { - w = 0.0 + y * I; - } - else { - r = fabs (x); - r = sqrt (r); - if (x < 0.0) { - w = 0.0 + copysign(r, y) * I; - } - else { - w = r + y * I; - } - } - return (w); - } - if (x == 0.0) { - r = fabs (y); - r = sqrt (0.5*r); - if (y > 0) - w = r + r * I; + /* Handle special cases. */ + if (z == 0) + return (CMPLX(0, b)); + if (isinf(b)) + return (CMPLX(INFINITY, b)); + if (isnan(a)) { + t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ + return (CMPLX(a + 0.0L + t, a + 0.0L + t)); /* NaN + NaN i */ + } + if (isinf(a)) { + /* + * csqrt(inf + NaN i) = inf + NaN i + * csqrt(inf + y i) = inf + 0 i + * csqrt(-inf + NaN i) = NaN +- inf i + * csqrt(-inf + y i) = 0 + inf i + */ + if (signbit(a)) + return (CMPLX(fabs(b - b), copysign(a, b))); else - w = r - r * I; - return (w); + return (CMPLX(a, copysign(b - b, b))); } - /* Rescale to avoid internal overflow or underflow. */ - if ((fabs(x) > 4.0) || (fabs(y) > 4.0)) { - x *= 0.25; - y *= 0.25; - scale = 2.0; - } - else { - x *= 1.8014398509481984e16; /* 2^54 */ - y *= 1.8014398509481984e16; - scale = 7.450580596923828125e-9; /* 2^-27 */ -#if 0 - x *= 4.0; - y *= 4.0; - scale = 0.5; -#endif + if (isnan(b)) { + t = (a - a) / (a - a); /* raise invalid */ + return (CMPLX(b + 0.0L + t, b + 0.0L + t)); /* NaN + NaN i */ + } + + /* Scale to avoid overflow. */ + if (fabs(a) >= THRESH || fabs(b) >= THRESH) { + /* + * Don't scale a or b if this might give (spurious) + * underflow. Then the unscaled value is an equivalent + * infinitesmal (or 0). + */ + if (fabs(a) >= 0x1p-1020) + a *= 0.25; + if (fabs(b) >= 0x1p-1020) + b *= 0.25; + scale = 2; + } else { + scale = 1; + } + + /* Scale to reduce inaccuracies when both components are denormal. */ + if (fabs(a) < 0x1p-1022 && fabs(b) < 0x1p-1022) { + a *= 0x1p54; + b *= 0x1p54; + scale = 0x1p-27; + } + + /* Algorithm 312, CACM vol 10, Oct 1967. */ + if (a >= 0) { + t = sqrt((a + hypot(a, b)) * 0.5); + rx = scale * t; + ry = scale * b / (2 * t); + } else { + t = sqrt((-a + hypot(a, b)) * 0.5); + rx = scale * fabs(b) / (2 * t); + ry = copysign(scale * t, b); } - w = x + y * I; - r = cabs(w); - if (x > 0) { - t = sqrt(0.5 * r + 0.5 * x); - r = scale * fabs((0.5 * y) / t); - t *= scale; - } - else { - r = sqrt( 0.5 * r - 0.5 * x ); - t = scale * fabs( (0.5 * y) / r ); - r *= scale; - } - if (y < 0) - w = t - r * I; - else - w = t + r * I; - return (w); + + return (CMPLX(rx, ry)); } + +#if LDBL_MANT_DIG == 53 +__weak_reference(csqrt, csqrtl); +#endif DEF_STD(csqrt); LDBL_MAYBE_CLONE(csqrt); Index: lib/libm/src/s_csqrtf.c =================================================================== RCS file: /cvs/src/lib/libm/src/s_csqrtf.c,v retrieving revision 1.4 diff -u -p -r1.4 s_csqrtf.c --- lib/libm/src/s_csqrtf.c 12 Sep 2016 19:47:02 -0000 1.4 +++ lib/libm/src/s_csqrtf.c 1 Jun 2019 09:22:46 -0000 @@ -1,132 +1,84 @@ -/* $OpenBSD: s_csqrtf.c,v 1.4 2016/09/12 19:47:02 guenther Exp $ */ -/* - * Copyright (c) 2008 Stephen L. Moshier +/*- + * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * - * Permission to use, copy, modify, and distribute this software for any - * purpose with or without fee is hereby granted, provided that the above - * copyright notice and this permission notice appear in all copies. + * Copyright (c) 2007 David Schultz + * All rights reserved. * - * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES - * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF - * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR - * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES - * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN - * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF - * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. */ -/* csqrtf() - * - * Complex square root - * - * - * - * SYNOPSIS: - * - * float complex csqrtf(); - * float complex z, w; - * - * w = csqrtf( z ); - * - * - * - * DESCRIPTION: - * - * - * If z = x + iy, r = |z|, then - * - * 1/2 - * Re w = [ (r + x)/2 ] , - * - * 1/2 - * Im w = [ (r - x)/2 ] . - * - * Cancellation error in r-x or r+x is avoided by using the - * identity 2 Re w Im w = y. - * - * Note that -w is also a square root of z. The root chosen - * is always in the right half plane and Im w has the same sign as y. - * - * - * - * ACCURACY: - * - * - * Relative error: - * arithmetic domain # trials peak rms - * IEEE -10,+10 1,000,000 1.8e-7 3.5e-8 - * - */ +#include #include #include +#include "math_private.h" + float complex csqrtf(float complex z) { - float complex w; - float x, y, r, t, scale; + double t; + float a, b; - x = crealf(z); - y = cimagf(z); + a = creal(z); + b = cimag(z); - if(y == 0.0f) { - if (x < 0.0f) { - w = 0.0f + copysign(sqrtf(-x), y) * I; - return (w); - } - else if (x == 0.0f) { - return (0.0f + y * I); - } - else { - w = sqrtf(x) + y * I; - return (w); - } - } - - if (x == 0.0f) { - r = fabsf(y); - r = sqrtf(0.5f*r); - if(y > 0) - w = r + r * I; + /* Handle special cases. */ + if (z == 0) + return (CMPLXF(0, b)); + if (isinf(b)) + return (CMPLXF(INFINITY, b)); + if (isnan(a)) { + t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ + return (CMPLXF(a + 0.0L + t, a + 0.0L + t)); /* NaN + NaN i */ + } + if (isinf(a)) { + /* + * csqrtf(inf + NaN i) = inf + NaN i + * csqrtf(inf + y i) = inf + 0 i + * csqrtf(-inf + NaN i) = NaN +- inf i + * csqrtf(-inf + y i) = 0 + inf i + */ + if (signbit(a)) + return (CMPLXF(fabsf(b - b), copysignf(a, b))); else - w = r - r * I; - return (w); - } - - /* Rescale to avoid internal overflow or underflow. */ - if ((fabsf(x) > 4.0f) || (fabsf(y) > 4.0f)) { - x *= 0.25f; - y *= 0.25f; - scale = 2.0f; + return (CMPLXF(a, copysignf(b - b, b))); } - else { - x *= 6.7108864e7f; /* 2^26 */ - y *= 6.7108864e7f; - scale = 1.220703125e-4f; /* 2^-13 */ -#if 0 - x *= 4.0f; - y *= 4.0f; - scale = 0.5f; -#endif + if (isnan(b)) { + t = (a - a) / (a - a); /* raise invalid */ + return (CMPLXF(b + 0.0L + t, b + 0.0L + t)); /* NaN + NaN i */ + } + + /* + * We compute t in double precision to avoid overflow and to + * provide correct rounding in nearly all cases. + * This is Algorithm 312, CACM vol 10, Oct 1967. + */ + if (a >= 0) { + t = sqrt((a + hypot(a, b)) * 0.5); + return (CMPLXF(t, b / (2 * t))); + } else { + t = sqrt((-a + hypot(a, b)) * 0.5); + return (CMPLXF(fabsf(b) / (2 * t), copysignf(t, b))); } - w = x + y * I; - r = cabsf(w); - if (x > 0) { - t = sqrtf( 0.5f * r + 0.5f * x ); - r = scale * fabsf((0.5f * y) / t); - t *= scale; - } - else { - r = sqrtf(0.5f * r - 0.5f * x); - t = scale * fabsf((0.5f * y) / r); - r *= scale; - } - - if (y < 0) - w = t - r * I; - else - w = t + r * I; - return (w); } DEF_STD(csqrtf); Index: lib/libm/src/s_csqrtl.c =================================================================== RCS file: /cvs/src/lib/libm/src/s_csqrtl.c,v retrieving revision 1.4 diff -u -p -r1.4 s_csqrtl.c --- lib/libm/src/s_csqrtl.c 12 Sep 2016 19:47:02 -0000 1.4 +++ lib/libm/src/s_csqrtl.c 1 Jun 2019 09:22:46 -0000 @@ -1,129 +1,126 @@ -/* $OpenBSD: s_csqrtl.c,v 1.4 2016/09/12 19:47:02 guenther Exp $ */ - -/* - * Copyright (c) 2008 Stephen L. Moshier +/*- + * SPDX-License-Identifier: BSD-2-Clause-FreeBSD * - * Permission to use, copy, modify, and distribute this software for any - * purpose with or without fee is hereby granted, provided that the above - * copyright notice and this permission notice appear in all copies. + * Copyright (c) 2007-2008 David Schultz + * All rights reserved. * - * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES - * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF - * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR - * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES - * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN - * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF - * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. */ -/* csqrtl() - * - * Complex square root - * - * - * - * SYNOPSIS: - * - * long double complex csqrtl(); - * long double complex z, w; - * - * w = csqrtl( z ); - * - * - * - * DESCRIPTION: - * - * - * If z = x + iy, r = |z|, then - * - * 1/2 - * Re w = [ (r + x)/2 ] , - * - * 1/2 - * Im w = [ (r - x)/2 ] . - * - * Cancellation error in r-x or r+x is avoided by using the - * identity 2 Re w Im w = y. - * - * Note that -w is also a square root of z. The root chosen - * is always in the right half plane and Im w has the same sign as y. - * - * - * - * ACCURACY: - * - * Relative error: - * arithmetic domain # trials peak rms - * IEEE -10,+10 500000 1.1e-19 3.0e-20 - * - */ +#include #include +#include #include +#include "math_private.h" + +/* + * Several thresholds require a 15-bit exponent and also the usual bias. + * s_logl.c and e_hypotl have less hard-coding but end up requiring the + * same for the exponent and more for the mantissa. + */ +#if LDBL_MAX_EXP != 0x4000 +#error "Unsupported long double format" +#endif + +/* + * Overflow must be avoided for components >= LDBL_MAX / (1 + sqrt(2)). + * The precise threshold is nontrivial to determine and spell, so use a + * lower threshold of approximaely LDBL_MAX / 4, and don't use LDBL_MAX + * to spell this since LDBL_MAX is broken on i386 (it overflows in 53-bit + * precision). + */ +#define THRESH 0x1p16382L + long double complex csqrtl(long double complex z) { - long double complex w; - long double x, y, r, t, scale; + long double complex result; + long double a, b, rx, ry, scale, t; - x = creall(z); - y = cimagl(z); + a = creall(z); + b = cimagl(z); - if (y == 0.0L) { - if (x < 0.0L) { - w = 0.0L + copysign(sqrtl(-x), y) * I; - return (w); - } - else { - w = sqrtl(x) + 0.0L * I; - return (w); - } - } - - if (x == 0.0L) { - r = fabsl(y); - r = sqrtl(0.5L * r); - if (y > 0.0L) - w = r + r * I; + /* Handle special cases. */ + if (z == 0) + return (CMPLXL(0, b)); + if (isinf(b)) + return (CMPLXL(INFINITY, b)); + if (isnan(a)) { + t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ + return (CMPLXL(a + 0.0L + t, a + 0.0L + t)); /* NaN + NaN i */ + } + if (isinf(a)) { + /* + * csqrt(inf + NaN i) = inf + NaN i + * csqrt(inf + y i) = inf + 0 i + * csqrt(-inf + NaN i) = NaN +- inf i + * csqrt(-inf + y i) = 0 + inf i + */ + if (signbit(a)) + return (CMPLXL(fabsl(b - b), copysignl(a, b))); else - w = r - r * I; - return (w); + return (CMPLXL(a, copysignl(b - b, b))); } - - /* Rescale to avoid internal overflow or underflow. */ - if ((fabsl(x) > 4.0L) || (fabsl(y) > 4.0L)) { - x *= 0.25L; - y *= 0.25L; - scale = 2.0L; - } - else { -#if 1 - x *= 7.3786976294838206464e19; /* 2^66 */ - y *= 7.3786976294838206464e19; - scale = 1.16415321826934814453125e-10; /* 2^-33 */ -#else - x *= 4.0L; - y *= 4.0L; - scale = 0.5L; -#endif + if (isnan(b)) { + t = (a - a) / (a - a); /* raise invalid */ + return (CMPLXL(b + 0.0L + t, b + 0.0L + t)); /* NaN + NaN i */ + } + + /* Scale to avoid overflow. */ + if (fabsl(a) >= THRESH || fabsl(b) >= THRESH) { + /* + * Don't scale a or b if this might give (spurious) + * underflow. Then the unscaled value is an equivalent + * infinitesmal (or 0). + */ + if (fabsl(a) >= 0x1p-16380L) + a *= 0.25; + if (fabsl(b) >= 0x1p-16380L) + b *= 0.25; + scale = 2; + } else { + scale = 1; + } + + /* Scale to reduce inaccuracies when both components are denormal. */ + if (fabsl(a) < 0x1p-16382L && fabsl(b) < 0x1p-16382L) { + a *= 0x1p64; + b *= 0x1p64; + scale = 0x1p-32; + } + + /* Algorithm 312, CACM vol 10, Oct 1967. */ + if (a >= 0) { + t = sqrtl((a + hypotl(a, b)) * 0.5); + rx = scale * t; + ry = scale * b / (2 * t); + } else { + t = sqrtl((-a + hypotl(a, b)) * 0.5); + rx = scale * fabsl(b) / (2 * t); + ry = copysignl(scale * t, b); } - w = x + y * I; - r = cabsl(w); - if (x > 0) { - t = sqrtl(0.5L * r + 0.5L * x); - r = scale * fabsl((0.5L * y) / t); - t *= scale; - } - else { - r = sqrtl(0.5L * r - 0.5L * x); - t = scale * fabsl((0.5L * y) / r); - r *= scale; - } - if (y < 0) - w = t - r * I; - else - w = t + r * I; - return (w); + + return (CMPLXL(rx, ry)); } DEF_STD(csqrtl); Index: regress/lib/libm/msun/Makefile =================================================================== RCS file: /cvs/src/regress/lib/libm/msun/Makefile,v retrieving revision 1.1.1.1 diff -u -p -r1.1.1.1 Makefile --- regress/lib/libm/msun/Makefile 21 Feb 2019 16:14:03 -0000 1.1.1.1 +++ regress/lib/libm/msun/Makefile 1 Jun 2019 09:40:10 -0000 @@ -2,6 +2,7 @@ TESTS = TESTS += conj_test +TESTS += csqrt_test TESTS += fenv_test TESTS += lrint_test Index: regress/lib/libm/msun/csqrt_test.c =================================================================== RCS file: regress/lib/libm/msun/csqrt_test.c diff -N regress/lib/libm/msun/csqrt_test.c --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ regress/lib/libm/msun/csqrt_test.c 1 Jun 2019 09:23:33 -0000 @@ -0,0 +1,377 @@ +/*- + * Copyright (c) 2007 David Schultz + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + */ + +/* + * Tests for csqrt{,f}() + */ + +#define rounddown(x,y) (((x)/(y))*(y)) + +#include + +#include + +#include +#include +#include +#include +#include + +#include "test-utils.h" + +/* + * This is a test hook that can point to csqrtl(), _csqrt(), or to _csqrtf(). + * The latter two convert to float or double, respectively, and test csqrtf() + * and csqrt() with the same arguments. + */ +static long double complex (*t_csqrt)(long double complex); + +static long double complex +_csqrtf(long double complex d) +{ + + return (csqrtf((float complex)d)); +} + +static long double complex +_csqrt(long double complex d) +{ + + return (csqrt((double complex)d)); +} + +#pragma STDC CX_LIMITED_RANGE OFF + +/* + * Compare d1 and d2 using special rules: NaN == NaN and +0 != -0. + * Fail an assertion if they differ. + */ +static void +assert_equal(long double complex d1, long double complex d2) +{ + + assert(cfpequal(d1, d2)); +} + +/* + * Test csqrt for some finite arguments where the answer is exact. + * (We do not test if it produces correctly rounded answers when the + * result is inexact, nor do we check whether it throws spurious + * exceptions.) + */ +static void +test_finite(void) +{ + static const double tests[] = { + /* csqrt(a + bI) = x + yI */ + /* a b x y */ + 0, 8, 2, 2, + 0, -8, 2, -2, + 4, 0, 2, 0, + -4, 0, 0, 2, + 3, 4, 2, 1, + 3, -4, 2, -1, + -3, 4, 1, 2, + -3, -4, 1, -2, + 5, 12, 3, 2, + 7, 24, 4, 3, + 9, 40, 5, 4, + 11, 60, 6, 5, + 13, 84, 7, 6, + 33, 56, 7, 4, + 39, 80, 8, 5, + 65, 72, 9, 4, + 987, 9916, 74, 67, + 5289, 6640, 83, 40, + 460766389075.0, 16762287900.0, 678910, 12345 + }; + /* + * We also test some multiples of the above arguments. This + * array defines which multiples we use. Note that these have + * to be small enough to not cause overflow for float precision + * with all of the constants in the above table. + */ + static const double mults[] = { + 1, + 2, + 3, + 13, + 16, + 0x1.p30, + 0x1.p-30, + }; + + double a, b; + double x, y; + unsigned i, j; + + for (i = 0; i < nitems(tests); i += 4) { + for (j = 0; j < nitems(mults); j++) { + a = tests[i] * mults[j] * mults[j]; + b = tests[i + 1] * mults[j] * mults[j]; + x = tests[i + 2] * mults[j]; + y = tests[i + 3] * mults[j]; + assert(t_csqrt(CMPLXL(a, b)) == CMPLXL(x, y)); + } + } + +} + +/* + * Test the handling of +/- 0. + */ +static void +test_zeros(void) +{ + + assert_equal(t_csqrt(CMPLXL(0.0, 0.0)), CMPLXL(0.0, 0.0)); + assert_equal(t_csqrt(CMPLXL(-0.0, 0.0)), CMPLXL(0.0, 0.0)); + assert_equal(t_csqrt(CMPLXL(0.0, -0.0)), CMPLXL(0.0, -0.0)); + assert_equal(t_csqrt(CMPLXL(-0.0, -0.0)), CMPLXL(0.0, -0.0)); +} + +/* + * Test the handling of infinities when the other argument is not NaN. + */ +static void +test_infinities(void) +{ + static const double vals[] = { + 0.0, + -0.0, + 42.0, + -42.0, + INFINITY, + -INFINITY, + }; + + unsigned i; + + for (i = 0; i < nitems(vals); i++) { + if (isfinite(vals[i])) { + assert_equal(t_csqrt(CMPLXL(-INFINITY, vals[i])), + CMPLXL(0.0, copysignl(INFINITY, vals[i]))); + assert_equal(t_csqrt(CMPLXL(INFINITY, vals[i])), + CMPLXL(INFINITY, copysignl(0.0, vals[i]))); + } + assert_equal(t_csqrt(CMPLXL(vals[i], INFINITY)), + CMPLXL(INFINITY, INFINITY)); + assert_equal(t_csqrt(CMPLXL(vals[i], -INFINITY)), + CMPLXL(INFINITY, -INFINITY)); + } +} + +/* + * Test the handling of NaNs. + */ +static void +test_nans(void) +{ + + assert(creall(t_csqrt(CMPLXL(INFINITY, NAN))) == INFINITY); + assert(isnan(cimagl(t_csqrt(CMPLXL(INFINITY, NAN))))); + + assert(isnan(creall(t_csqrt(CMPLXL(-INFINITY, NAN))))); + assert(isinf(cimagl(t_csqrt(CMPLXL(-INFINITY, NAN))))); + + assert_equal(t_csqrt(CMPLXL(NAN, INFINITY)), + CMPLXL(INFINITY, INFINITY)); + assert_equal(t_csqrt(CMPLXL(NAN, -INFINITY)), + CMPLXL(INFINITY, -INFINITY)); + + assert_equal(t_csqrt(CMPLXL(0.0, NAN)), CMPLXL(NAN, NAN)); + assert_equal(t_csqrt(CMPLXL(-0.0, NAN)), CMPLXL(NAN, NAN)); + assert_equal(t_csqrt(CMPLXL(42.0, NAN)), CMPLXL(NAN, NAN)); + assert_equal(t_csqrt(CMPLXL(-42.0, NAN)), CMPLXL(NAN, NAN)); + assert_equal(t_csqrt(CMPLXL(NAN, 0.0)), CMPLXL(NAN, NAN)); + assert_equal(t_csqrt(CMPLXL(NAN, -0.0)), CMPLXL(NAN, NAN)); + assert_equal(t_csqrt(CMPLXL(NAN, 42.0)), CMPLXL(NAN, NAN)); + assert_equal(t_csqrt(CMPLXL(NAN, -42.0)), CMPLXL(NAN, NAN)); + assert_equal(t_csqrt(CMPLXL(NAN, NAN)), CMPLXL(NAN, NAN)); +} + +/* + * Test whether csqrt(a + bi) works for inputs that are large enough to + * cause overflow in hypot(a, b) + a. Each of the tests is scaled up to + * near MAX_EXP. + */ +static void +test_overflow(int maxexp) +{ + long double a, b; + long double complex result; + int exp, i; + + assert(maxexp > 0 && maxexp % 2 == 0); + + for (i = 0; i < 4; i++) { + exp = maxexp - 2 * i; + + /* csqrt(115 + 252*I) == 14 + 9*I */ + a = ldexpl(115 * 0x1p-8, exp); + b = ldexpl(252 * 0x1p-8, exp); + result = t_csqrt(CMPLXL(a, b)); + assert(creall(result) == ldexpl(14 * 0x1p-4, exp / 2)); + assert(cimagl(result) == ldexpl(9 * 0x1p-4, exp / 2)); + + /* csqrt(-11 + 60*I) = 5 + 6*I */ + a = ldexpl(-11 * 0x1p-6, exp); + b = ldexpl(60 * 0x1p-6, exp); + result = t_csqrt(CMPLXL(a, b)); + assert(creall(result) == ldexpl(5 * 0x1p-3, exp / 2)); + assert(cimagl(result) == ldexpl(6 * 0x1p-3, exp / 2)); + + /* csqrt(225 + 0*I) == 15 + 0*I */ + a = ldexpl(225 * 0x1p-8, exp); + b = 0; + result = t_csqrt(CMPLXL(a, b)); + assert(creall(result) == ldexpl(15 * 0x1p-4, exp / 2)); + assert(cimagl(result) == 0); + } +} + +/* + * Test that precision is maintained for some large squares. Set all or + * some bits in the lower mantdig/2 bits, square the number, and try to + * recover the sqrt. Note: + * (x + xI)**2 = 2xxI + */ +static void +test_precision(int maxexp, int mantdig) +{ + long double b, x; + long double complex result; + uint64_t mantbits, sq_mantbits; + int exp, i; + + assert(maxexp > 0 && maxexp % 2 == 0); + assert(mantdig <= 64); + mantdig = rounddown(mantdig, 2); + + for (exp = 0; exp <= maxexp; exp += 2) { + mantbits = ((uint64_t)1 << (mantdig / 2 )) - 1; + for (i = 0; + i < 100 && mantbits > ((uint64_t)1 << (mantdig / 2 - 1)); + i++, mantbits--) { + sq_mantbits = mantbits * mantbits; + /* + * sq_mantibts is a mantdig-bit number. Divide by + * 2**mantdig to normalize it to [0.5, 1), where, + * note, the binary power will be -1. Raise it by + * 2**exp for the test. exp is even. Lower it by + * one to reach a final binary power which is also + * even. The result should be exactly + * representable, given that mantdig is less than or + * equal to the available precision. + */ + b = ldexpl((long double)sq_mantbits, + exp - 1 - mantdig); + x = ldexpl(mantbits, (exp - 2 - mantdig) / 2); + assert(b == x * x * 2); + result = t_csqrt(CMPLXL(0, b)); + assert(creall(result) == x); + assert(cimagl(result) == x); + } + } +} + +int +main(void) +{ + + printf("1..18\n"); + + /* Test csqrt() */ + t_csqrt = _csqrt; + + test_finite(); + printf("ok 1 - csqrt\n"); + + test_zeros(); + printf("ok 2 - csqrt\n"); + + test_infinities(); + printf("ok 3 - csqrt\n"); + + test_nans(); + printf("ok 4 - csqrt\n"); + + test_overflow(DBL_MAX_EXP); + printf("ok 5 - csqrt\n"); + + test_precision(DBL_MAX_EXP, DBL_MANT_DIG); + printf("ok 6 - csqrt\n"); + + /* Now test csqrtf() */ + t_csqrt = _csqrtf; + + test_finite(); + printf("ok 7 - csqrt\n"); + + test_zeros(); + printf("ok 8 - csqrt\n"); + + test_infinities(); + printf("ok 9 - csqrt\n"); + + test_nans(); + printf("ok 10 - csqrt\n"); + + test_overflow(FLT_MAX_EXP); + printf("ok 11 - csqrt\n"); + + test_precision(FLT_MAX_EXP, FLT_MANT_DIG); + printf("ok 12 - csqrt\n"); + + /* Now test csqrtl() */ + t_csqrt = csqrtl; + + test_finite(); + printf("ok 13 - csqrt\n"); + + test_zeros(); + printf("ok 14 - csqrt\n"); + + test_infinities(); + printf("ok 15 - csqrt\n"); + + test_nans(); + printf("ok 16 - csqrt\n"); + + test_overflow(LDBL_MAX_EXP); + printf("ok 17 - csqrt\n"); + + test_precision(LDBL_MAX_EXP, +#ifndef __i386__ + LDBL_MANT_DIG +#else + DBL_MANT_DIG +#endif + ); + printf("ok 18 - csqrt\n"); + + return (0); +}