+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]
--- /dev/null
+/* 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
--- /dev/null
+/* 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
--- /dev/null
+/* 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
--- /dev/null
+/* 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
--- /dev/null
+/* 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
--- /dev/null
+/* 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
projects / thirdparty / glibc.git / commitdiff
