Project

General

Profile

Download (29.3 KB)

Task #6746 » csqrt.diff

Moritz Buhl, 06/02/2019 11:51 AM

View differences:

lib/libm/src/s_csqrt.c 1 Jun 2019 09:22:46 -0000
/* $OpenBSD: s_csqrt.c,v 1.8 2016/09/12 19:47:02 guenther Exp $ */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
/*-
* 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 <das@FreeBSD.ORG>
* 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 <sys/cdefs.h>
#include <complex.h>
#include <float.h>
#include <math.h>
#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);
lib/libm/src/s_csqrtf.c 1 Jun 2019 09:22:46 -0000
/* $OpenBSD: s_csqrtf.c,v 1.4 2016/09/12 19:47:02 guenther Exp $ */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
/*-
* 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 <das@FreeBSD.ORG>
* 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 <sys/cdefs.h>
#include <complex.h>
#include <math.h>
#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);
lib/libm/src/s_csqrtl.c 1 Jun 2019 09:22:46 -0000
/* $OpenBSD: s_csqrtl.c,v 1.4 2016/09/12 19:47:02 guenther Exp $ */
/*
* Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
/*-
* 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 <das@FreeBSD.ORG>
* 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 <sys/cdefs.h>
#include <complex.h>
#include <float.h>
#include <math.h>
#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);
regress/lib/libm/msun/Makefile 1 Jun 2019 09:40:10 -0000
TESTS =
TESTS += conj_test
TESTS += csqrt_test
TESTS += fenv_test
TESTS += lrint_test
regress/lib/libm/msun/csqrt_test.c 1 Jun 2019 09:23:33 -0000
/*-
* Copyright (c) 2007 David Schultz <das@FreeBSD.org>
* 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 <sys/cdefs.h>
#include <sys/param.h>
#include <assert.h>
#include <complex.h>
#include <float.h>
#include <math.h>
#include <stdio.h>
#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);
}
(3-3/4)