Prepare to convert remaining _Complex functions
This patch has no function changes, except to ensure the git history correctly tracks the changes to convert the double version of these functions into a templated version.
This commit is contained in:
Paul E. Murphy
parent
d47d27d6c0
commit
1dbc54f61e
|
|
@ -1,3 +1,12 @@
|
|||
2016-08-29 Paul E. Murphy <murphyp@linux.vnet.ibm.com>
|
||||
|
||||
* s_cexp_template.c: Copy of s_cexp.c.
|
||||
* s_clog_template.c: Copy of s_clog.c.
|
||||
* s_clog10_template.c: Copy of s_clog10.c.
|
||||
* s_cpow_template.c: Copy of s_cpow.c.
|
||||
* s_cproj_template.c: Copy of s_cproj.c.
|
||||
* s_csqrt_template.c: Copy of s_csqrt.c.
|
||||
|
||||
2016-08-29 Paul E. Murphy <murphyp@linux.vnet.ibm.com>
|
||||
|
||||
[BZ #20517]
|
||||
|
|
|
|||
|
|
@ -0,0 +1,157 @@
|
|||
/* Return value of complex exponential function for double complex value.
|
||||
Copyright (C) 1997-2016 Free Software Foundation, Inc.
|
||||
This file is part of the GNU C Library.
|
||||
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
|
||||
|
||||
The GNU C Library is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU Lesser General Public
|
||||
License as published by the Free Software Foundation; either
|
||||
version 2.1 of the License, or (at your option) any later version.
|
||||
|
||||
The GNU C Library is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
Lesser General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU Lesser General Public
|
||||
License along with the GNU C Library; if not, see
|
||||
<http://www.gnu.org/licenses/>. */
|
||||
|
||||
#include <complex.h>
|
||||
#include <fenv.h>
|
||||
#include <math.h>
|
||||
#include <math_private.h>
|
||||
#include <float.h>
|
||||
|
||||
__complex__ double
|
||||
__cexp (__complex__ double x)
|
||||
{
|
||||
__complex__ double retval;
|
||||
int rcls = fpclassify (__real__ x);
|
||||
int icls = fpclassify (__imag__ x);
|
||||
|
||||
if (__glibc_likely (rcls >= FP_ZERO))
|
||||
{
|
||||
/* Real part is finite. */
|
||||
if (__glibc_likely (icls >= FP_ZERO))
|
||||
{
|
||||
/* Imaginary part is finite. */
|
||||
const int t = (int) ((DBL_MAX_EXP - 1) * M_LN2);
|
||||
double sinix, cosix;
|
||||
|
||||
if (__glibc_likely (fabs (__imag__ x) > DBL_MIN))
|
||||
{
|
||||
__sincos (__imag__ x, &sinix, &cosix);
|
||||
}
|
||||
else
|
||||
{
|
||||
sinix = __imag__ x;
|
||||
cosix = 1.0;
|
||||
}
|
||||
|
||||
if (__real__ x > t)
|
||||
{
|
||||
double exp_t = __ieee754_exp (t);
|
||||
__real__ x -= t;
|
||||
sinix *= exp_t;
|
||||
cosix *= exp_t;
|
||||
if (__real__ x > t)
|
||||
{
|
||||
__real__ x -= t;
|
||||
sinix *= exp_t;
|
||||
cosix *= exp_t;
|
||||
}
|
||||
}
|
||||
if (__real__ x > t)
|
||||
{
|
||||
/* Overflow (original real part of x > 3t). */
|
||||
__real__ retval = DBL_MAX * cosix;
|
||||
__imag__ retval = DBL_MAX * sinix;
|
||||
}
|
||||
else
|
||||
{
|
||||
double exp_val = __ieee754_exp (__real__ x);
|
||||
__real__ retval = exp_val * cosix;
|
||||
__imag__ retval = exp_val * sinix;
|
||||
}
|
||||
math_check_force_underflow_complex (retval);
|
||||
}
|
||||
else
|
||||
{
|
||||
/* If the imaginary part is +-inf or NaN and the real part
|
||||
is not +-inf the result is NaN + iNaN. */
|
||||
__real__ retval = __nan ("");
|
||||
__imag__ retval = __nan ("");
|
||||
|
||||
feraiseexcept (FE_INVALID);
|
||||
}
|
||||
}
|
||||
else if (__glibc_likely (rcls == FP_INFINITE))
|
||||
{
|
||||
/* Real part is infinite. */
|
||||
if (__glibc_likely (icls >= FP_ZERO))
|
||||
{
|
||||
/* Imaginary part is finite. */
|
||||
double value = signbit (__real__ x) ? 0.0 : HUGE_VAL;
|
||||
|
||||
if (icls == FP_ZERO)
|
||||
{
|
||||
/* Imaginary part is 0.0. */
|
||||
__real__ retval = value;
|
||||
__imag__ retval = __imag__ x;
|
||||
}
|
||||
else
|
||||
{
|
||||
double sinix, cosix;
|
||||
|
||||
if (__glibc_likely (fabs (__imag__ x) > DBL_MIN))
|
||||
{
|
||||
__sincos (__imag__ x, &sinix, &cosix);
|
||||
}
|
||||
else
|
||||
{
|
||||
sinix = __imag__ x;
|
||||
cosix = 1.0;
|
||||
}
|
||||
|
||||
__real__ retval = __copysign (value, cosix);
|
||||
__imag__ retval = __copysign (value, sinix);
|
||||
}
|
||||
}
|
||||
else if (signbit (__real__ x) == 0)
|
||||
{
|
||||
__real__ retval = HUGE_VAL;
|
||||
__imag__ retval = __nan ("");
|
||||
|
||||
if (icls == FP_INFINITE)
|
||||
feraiseexcept (FE_INVALID);
|
||||
}
|
||||
else
|
||||
{
|
||||
__real__ retval = 0.0;
|
||||
__imag__ retval = __copysign (0.0, __imag__ x);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
/* If the real part is NaN the result is NaN + iNaN unless the
|
||||
imaginary part is zero. */
|
||||
__real__ retval = __nan ("");
|
||||
if (icls == FP_ZERO)
|
||||
__imag__ retval = __imag__ x;
|
||||
else
|
||||
{
|
||||
__imag__ retval = __nan ("");
|
||||
|
||||
if (rcls != FP_NAN || icls != FP_NAN)
|
||||
feraiseexcept (FE_INVALID);
|
||||
}
|
||||
}
|
||||
|
||||
return retval;
|
||||
}
|
||||
weak_alias (__cexp, cexp)
|
||||
#ifdef NO_LONG_DOUBLE
|
||||
strong_alias (__cexp, __cexpl)
|
||||
weak_alias (__cexp, cexpl)
|
||||
#endif
|
||||
|
|
@ -0,0 +1,124 @@
|
|||
/* Compute complex base 10 logarithm.
|
||||
Copyright (C) 1997-2016 Free Software Foundation, Inc.
|
||||
This file is part of the GNU C Library.
|
||||
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
|
||||
|
||||
The GNU C Library is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU Lesser General Public
|
||||
License as published by the Free Software Foundation; either
|
||||
version 2.1 of the License, or (at your option) any later version.
|
||||
|
||||
The GNU C Library is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
Lesser General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU Lesser General Public
|
||||
License along with the GNU C Library; if not, see
|
||||
<http://www.gnu.org/licenses/>. */
|
||||
|
||||
#include <complex.h>
|
||||
#include <math.h>
|
||||
#include <math_private.h>
|
||||
#include <float.h>
|
||||
|
||||
/* log_10 (2). */
|
||||
#define M_LOG10_2 0.3010299956639811952137388947244930267682
|
||||
|
||||
/* pi * log10 (e). */
|
||||
#define M_PI_LOG10E 1.364376353841841347485783625431355770210
|
||||
|
||||
__complex__ double
|
||||
__clog10 (__complex__ double x)
|
||||
{
|
||||
__complex__ double result;
|
||||
int rcls = fpclassify (__real__ x);
|
||||
int icls = fpclassify (__imag__ x);
|
||||
|
||||
if (__glibc_unlikely (rcls == FP_ZERO && icls == FP_ZERO))
|
||||
{
|
||||
/* Real and imaginary part are 0.0. */
|
||||
__imag__ result = signbit (__real__ x) ? M_PI_LOG10E : 0.0;
|
||||
__imag__ result = __copysign (__imag__ result, __imag__ x);
|
||||
/* Yes, the following line raises an exception. */
|
||||
__real__ result = -1.0 / fabs (__real__ x);
|
||||
}
|
||||
else if (__glibc_likely (rcls != FP_NAN && icls != FP_NAN))
|
||||
{
|
||||
/* Neither real nor imaginary part is NaN. */
|
||||
double absx = fabs (__real__ x), absy = fabs (__imag__ x);
|
||||
int scale = 0;
|
||||
|
||||
if (absx < absy)
|
||||
{
|
||||
double t = absx;
|
||||
absx = absy;
|
||||
absy = t;
|
||||
}
|
||||
|
||||
if (absx > DBL_MAX / 2.0)
|
||||
{
|
||||
scale = -1;
|
||||
absx = __scalbn (absx, scale);
|
||||
absy = (absy >= DBL_MIN * 2.0 ? __scalbn (absy, scale) : 0.0);
|
||||
}
|
||||
else if (absx < DBL_MIN && absy < DBL_MIN)
|
||||
{
|
||||
scale = DBL_MANT_DIG;
|
||||
absx = __scalbn (absx, scale);
|
||||
absy = __scalbn (absy, scale);
|
||||
}
|
||||
|
||||
if (absx == 1.0 && scale == 0)
|
||||
{
|
||||
__real__ result = __log1p (absy * absy) * (M_LOG10E / 2.0);
|
||||
math_check_force_underflow_nonneg (__real__ result);
|
||||
}
|
||||
else if (absx > 1.0 && absx < 2.0 && absy < 1.0 && scale == 0)
|
||||
{
|
||||
double d2m1 = (absx - 1.0) * (absx + 1.0);
|
||||
if (absy >= DBL_EPSILON)
|
||||
d2m1 += absy * absy;
|
||||
__real__ result = __log1p (d2m1) * (M_LOG10E / 2.0);
|
||||
}
|
||||
else if (absx < 1.0
|
||||
&& absx >= 0.5
|
||||
&& absy < DBL_EPSILON / 2.0
|
||||
&& scale == 0)
|
||||
{
|
||||
double d2m1 = (absx - 1.0) * (absx + 1.0);
|
||||
__real__ result = __log1p (d2m1) * (M_LOG10E / 2.0);
|
||||
}
|
||||
else if (absx < 1.0
|
||||
&& absx >= 0.5
|
||||
&& scale == 0
|
||||
&& absx * absx + absy * absy >= 0.5)
|
||||
{
|
||||
double d2m1 = __x2y2m1 (absx, absy);
|
||||
__real__ result = __log1p (d2m1) * (M_LOG10E / 2.0);
|
||||
}
|
||||
else
|
||||
{
|
||||
double d = __ieee754_hypot (absx, absy);
|
||||
__real__ result = __ieee754_log10 (d) - scale * M_LOG10_2;
|
||||
}
|
||||
|
||||
__imag__ result = M_LOG10E * __ieee754_atan2 (__imag__ x, __real__ x);
|
||||
}
|
||||
else
|
||||
{
|
||||
__imag__ result = __nan ("");
|
||||
if (rcls == FP_INFINITE || icls == FP_INFINITE)
|
||||
/* Real or imaginary part is infinite. */
|
||||
__real__ result = HUGE_VAL;
|
||||
else
|
||||
__real__ result = __nan ("");
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
weak_alias (__clog10, clog10)
|
||||
#ifdef NO_LONG_DOUBLE
|
||||
strong_alias (__clog10, __clog10l)
|
||||
weak_alias (__clog10, clog10l)
|
||||
#endif
|
||||
|
|
@ -0,0 +1,118 @@
|
|||
/* Compute complex natural logarithm.
|
||||
Copyright (C) 1997-2016 Free Software Foundation, Inc.
|
||||
This file is part of the GNU C Library.
|
||||
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
|
||||
|
||||
The GNU C Library is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU Lesser General Public
|
||||
License as published by the Free Software Foundation; either
|
||||
version 2.1 of the License, or (at your option) any later version.
|
||||
|
||||
The GNU C Library is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
Lesser General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU Lesser General Public
|
||||
License along with the GNU C Library; if not, see
|
||||
<http://www.gnu.org/licenses/>. */
|
||||
|
||||
#include <complex.h>
|
||||
#include <math.h>
|
||||
#include <math_private.h>
|
||||
#include <float.h>
|
||||
|
||||
__complex__ double
|
||||
__clog (__complex__ double x)
|
||||
{
|
||||
__complex__ double result;
|
||||
int rcls = fpclassify (__real__ x);
|
||||
int icls = fpclassify (__imag__ x);
|
||||
|
||||
if (__glibc_unlikely (rcls == FP_ZERO && icls == FP_ZERO))
|
||||
{
|
||||
/* Real and imaginary part are 0.0. */
|
||||
__imag__ result = signbit (__real__ x) ? M_PI : 0.0;
|
||||
__imag__ result = __copysign (__imag__ result, __imag__ x);
|
||||
/* Yes, the following line raises an exception. */
|
||||
__real__ result = -1.0 / fabs (__real__ x);
|
||||
}
|
||||
else if (__glibc_likely (rcls != FP_NAN && icls != FP_NAN))
|
||||
{
|
||||
/* Neither real nor imaginary part is NaN. */
|
||||
double absx = fabs (__real__ x), absy = fabs (__imag__ x);
|
||||
int scale = 0;
|
||||
|
||||
if (absx < absy)
|
||||
{
|
||||
double t = absx;
|
||||
absx = absy;
|
||||
absy = t;
|
||||
}
|
||||
|
||||
if (absx > DBL_MAX / 2.0)
|
||||
{
|
||||
scale = -1;
|
||||
absx = __scalbn (absx, scale);
|
||||
absy = (absy >= DBL_MIN * 2.0 ? __scalbn (absy, scale) : 0.0);
|
||||
}
|
||||
else if (absx < DBL_MIN && absy < DBL_MIN)
|
||||
{
|
||||
scale = DBL_MANT_DIG;
|
||||
absx = __scalbn (absx, scale);
|
||||
absy = __scalbn (absy, scale);
|
||||
}
|
||||
|
||||
if (absx == 1.0 && scale == 0)
|
||||
{
|
||||
__real__ result = __log1p (absy * absy) / 2.0;
|
||||
math_check_force_underflow_nonneg (__real__ result);
|
||||
}
|
||||
else if (absx > 1.0 && absx < 2.0 && absy < 1.0 && scale == 0)
|
||||
{
|
||||
double d2m1 = (absx - 1.0) * (absx + 1.0);
|
||||
if (absy >= DBL_EPSILON)
|
||||
d2m1 += absy * absy;
|
||||
__real__ result = __log1p (d2m1) / 2.0;
|
||||
}
|
||||
else if (absx < 1.0
|
||||
&& absx >= 0.5
|
||||
&& absy < DBL_EPSILON / 2.0
|
||||
&& scale == 0)
|
||||
{
|
||||
double d2m1 = (absx - 1.0) * (absx + 1.0);
|
||||
__real__ result = __log1p (d2m1) / 2.0;
|
||||
}
|
||||
else if (absx < 1.0
|
||||
&& absx >= 0.5
|
||||
&& scale == 0
|
||||
&& absx * absx + absy * absy >= 0.5)
|
||||
{
|
||||
double d2m1 = __x2y2m1 (absx, absy);
|
||||
__real__ result = __log1p (d2m1) / 2.0;
|
||||
}
|
||||
else
|
||||
{
|
||||
double d = __ieee754_hypot (absx, absy);
|
||||
__real__ result = __ieee754_log (d) - scale * M_LN2;
|
||||
}
|
||||
|
||||
__imag__ result = __ieee754_atan2 (__imag__ x, __real__ x);
|
||||
}
|
||||
else
|
||||
{
|
||||
__imag__ result = __nan ("");
|
||||
if (rcls == FP_INFINITE || icls == FP_INFINITE)
|
||||
/* Real or imaginary part is infinite. */
|
||||
__real__ result = HUGE_VAL;
|
||||
else
|
||||
__real__ result = __nan ("");
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
weak_alias (__clog, clog)
|
||||
#ifdef NO_LONG_DOUBLE
|
||||
strong_alias (__clog, __clogl)
|
||||
weak_alias (__clog, clogl)
|
||||
#endif
|
||||
|
|
@ -0,0 +1,33 @@
|
|||
/* Complex power of double values.
|
||||
Copyright (C) 1997-2016 Free Software Foundation, Inc.
|
||||
This file is part of the GNU C Library.
|
||||
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
|
||||
|
||||
The GNU C Library is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU Lesser General Public
|
||||
License as published by the Free Software Foundation; either
|
||||
version 2.1 of the License, or (at your option) any later version.
|
||||
|
||||
The GNU C Library is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
Lesser General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU Lesser General Public
|
||||
License along with the GNU C Library; if not, see
|
||||
<http://www.gnu.org/licenses/>. */
|
||||
|
||||
#include <complex.h>
|
||||
#include <math.h>
|
||||
|
||||
|
||||
__complex__ double
|
||||
__cpow (__complex__ double x, __complex__ double c)
|
||||
{
|
||||
return __cexp (c * __clog (x));
|
||||
}
|
||||
weak_alias (__cpow, cpow)
|
||||
#ifdef NO_LONG_DOUBLE
|
||||
strong_alias (__cpow, __cpowl)
|
||||
weak_alias (__cpow, cpowl)
|
||||
#endif
|
||||
|
|
@ -0,0 +1,44 @@
|
|||
/* Compute projection of complex double value to Riemann sphere.
|
||||
Copyright (C) 1997-2016 Free Software Foundation, Inc.
|
||||
This file is part of the GNU C Library.
|
||||
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
|
||||
|
||||
The GNU C Library is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU Lesser General Public
|
||||
License as published by the Free Software Foundation; either
|
||||
version 2.1 of the License, or (at your option) any later version.
|
||||
|
||||
The GNU C Library is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
Lesser General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU Lesser General Public
|
||||
License along with the GNU C Library; if not, see
|
||||
<http://www.gnu.org/licenses/>. */
|
||||
|
||||
#include <complex.h>
|
||||
#include <math.h>
|
||||
#include <math_private.h>
|
||||
|
||||
|
||||
__complex__ double
|
||||
__cproj (__complex__ double x)
|
||||
{
|
||||
if (isinf (__real__ x) || isinf (__imag__ x))
|
||||
{
|
||||
__complex__ double res;
|
||||
|
||||
__real__ res = INFINITY;
|
||||
__imag__ res = __copysign (0.0, __imag__ x);
|
||||
|
||||
return res;
|
||||
}
|
||||
|
||||
return x;
|
||||
}
|
||||
weak_alias (__cproj, cproj)
|
||||
#ifdef NO_LONG_DOUBLE
|
||||
strong_alias (__cproj, __cprojl)
|
||||
weak_alias (__cproj, cprojl)
|
||||
#endif
|
||||
|
|
@ -0,0 +1,165 @@
|
|||
/* Complex square root of double value.
|
||||
Copyright (C) 1997-2016 Free Software Foundation, Inc.
|
||||
This file is part of the GNU C Library.
|
||||
Based on an algorithm by Stephen L. Moshier <moshier@world.std.com>.
|
||||
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
|
||||
|
||||
The GNU C Library is free software; you can redistribute it and/or
|
||||
modify it under the terms of the GNU Lesser General Public
|
||||
License as published by the Free Software Foundation; either
|
||||
version 2.1 of the License, or (at your option) any later version.
|
||||
|
||||
The GNU C Library is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
Lesser General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU Lesser General Public
|
||||
License along with the GNU C Library; if not, see
|
||||
<http://www.gnu.org/licenses/>. */
|
||||
|
||||
#include <complex.h>
|
||||
#include <math.h>
|
||||
#include <math_private.h>
|
||||
#include <float.h>
|
||||
|
||||
__complex__ double
|
||||
__csqrt (__complex__ double x)
|
||||
{
|
||||
__complex__ double res;
|
||||
int rcls = fpclassify (__real__ x);
|
||||
int icls = fpclassify (__imag__ x);
|
||||
|
||||
if (__glibc_unlikely (rcls <= FP_INFINITE || icls <= FP_INFINITE))
|
||||
{
|
||||
if (icls == FP_INFINITE)
|
||||
{
|
||||
__real__ res = HUGE_VAL;
|
||||
__imag__ res = __imag__ x;
|
||||
}
|
||||
else if (rcls == FP_INFINITE)
|
||||
{
|
||||
if (__real__ x < 0.0)
|
||||
{
|
||||
__real__ res = icls == FP_NAN ? __nan ("") : 0;
|
||||
__imag__ res = __copysign (HUGE_VAL, __imag__ x);
|
||||
}
|
||||
else
|
||||
{
|
||||
__real__ res = __real__ x;
|
||||
__imag__ res = (icls == FP_NAN
|
||||
? __nan ("") : __copysign (0.0, __imag__ x));
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
__real__ res = __nan ("");
|
||||
__imag__ res = __nan ("");
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if (__glibc_unlikely (icls == FP_ZERO))
|
||||
{
|
||||
if (__real__ x < 0.0)
|
||||
{
|
||||
__real__ res = 0.0;
|
||||
__imag__ res = __copysign (__ieee754_sqrt (-__real__ x),
|
||||
__imag__ x);
|
||||
}
|
||||
else
|
||||
{
|
||||
__real__ res = fabs (__ieee754_sqrt (__real__ x));
|
||||
__imag__ res = __copysign (0.0, __imag__ x);
|
||||
}
|
||||
}
|
||||
else if (__glibc_unlikely (rcls == FP_ZERO))
|
||||
{
|
||||
double r;
|
||||
if (fabs (__imag__ x) >= 2.0 * DBL_MIN)
|
||||
r = __ieee754_sqrt (0.5 * fabs (__imag__ x));
|
||||
else
|
||||
r = 0.5 * __ieee754_sqrt (2.0 * fabs (__imag__ x));
|
||||
|
||||
__real__ res = r;
|
||||
__imag__ res = __copysign (r, __imag__ x);
|
||||
}
|
||||
else
|
||||
{
|
||||
double d, r, s;
|
||||
int scale = 0;
|
||||
|
||||
if (fabs (__real__ x) > DBL_MAX / 4.0)
|
||||
{
|
||||
scale = 1;
|
||||
__real__ x = __scalbn (__real__ x, -2 * scale);
|
||||
__imag__ x = __scalbn (__imag__ x, -2 * scale);
|
||||
}
|
||||
else if (fabs (__imag__ x) > DBL_MAX / 4.0)
|
||||
{
|
||||
scale = 1;
|
||||
if (fabs (__real__ x) >= 4.0 * DBL_MIN)
|
||||
__real__ x = __scalbn (__real__ x, -2 * scale);
|
||||
else
|
||||
__real__ x = 0.0;
|
||||
__imag__ x = __scalbn (__imag__ x, -2 * scale);
|
||||
}
|
||||
else if (fabs (__real__ x) < 2.0 * DBL_MIN
|
||||
&& fabs (__imag__ x) < 2.0 * DBL_MIN)
|
||||
{
|
||||
scale = -((DBL_MANT_DIG + 1) / 2);
|
||||
__real__ x = __scalbn (__real__ x, -2 * scale);
|
||||
__imag__ x = __scalbn (__imag__ x, -2 * scale);
|
||||
}
|
||||
|
||||
d = __ieee754_hypot (__real__ x, __imag__ x);
|
||||
/* Use the identity 2 Re res Im res = Im x
|
||||
to avoid cancellation error in d +/- Re x. */
|
||||
if (__real__ x > 0)
|
||||
{
|
||||
r = __ieee754_sqrt (0.5 * (d + __real__ x));
|
||||
if (scale == 1 && fabs (__imag__ x) < 1.0)
|
||||
{
|
||||
/* Avoid possible intermediate underflow. */
|
||||
s = __imag__ x / r;
|
||||
r = __scalbn (r, scale);
|
||||
scale = 0;
|
||||
}
|
||||
else
|
||||
s = 0.5 * (__imag__ x / r);
|
||||
}
|
||||
else
|
||||
{
|
||||
s = __ieee754_sqrt (0.5 * (d - __real__ x));
|
||||
if (scale == 1 && fabs (__imag__ x) < 1.0)
|
||||
{
|
||||
/* Avoid possible intermediate underflow. */
|
||||
r = fabs (__imag__ x / s);
|
||||
s = __scalbn (s, scale);
|
||||
scale = 0;
|
||||
}
|
||||
else
|
||||
r = fabs (0.5 * (__imag__ x / s));
|
||||
}
|
||||
|
||||
if (scale)
|
||||
{
|
||||
r = __scalbn (r, scale);
|
||||
s = __scalbn (s, scale);
|
||||
}
|
||||
|
||||
math_check_force_underflow (r);
|
||||
math_check_force_underflow (s);
|
||||
|
||||
__real__ res = r;
|
||||
__imag__ res = __copysign (s, __imag__ x);
|
||||
}
|
||||
}
|
||||
|
||||
return res;
|
||||
}
|
||||
weak_alias (__csqrt, csqrt)
|
||||
#ifdef NO_LONG_DOUBLE
|
||||
strong_alias (__csqrt, __csqrtl)
|
||||
weak_alias (__csqrt, csqrtl)
|
||||
#endif
|
||||
Loading…
Reference in New Issue