diff options
| -rw-r--r-- | ChangeLog | 9 | ||||
| -rw-r--r-- | math/s_cexp_template.c | 157 | ||||
| -rw-r--r-- | math/s_clog10_template.c | 124 | ||||
| -rw-r--r-- | math/s_clog_template.c | 118 | ||||
| -rw-r--r-- | math/s_cpow_template.c | 33 | ||||
| -rw-r--r-- | math/s_cproj_template.c | 44 | ||||
| -rw-r--r-- | math/s_csqrt_template.c | 165 |
7 files changed, 650 insertions, 0 deletions
diff --git a/ChangeLog b/ChangeLog index 7438924c7a..6932deb11a 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,14 @@ 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] * sysdeps/sparc/sparc32/sparcv9/fpu/multiarch/s_fdim.c: Include ldbl-opt/s_fdim.c to ensure fdiml@GLIBC_2.1 is diff --git a/math/s_cexp_template.c b/math/s_cexp_template.c new file mode 100644 index 0000000000..3a476bde3c --- /dev/null +++ b/math/s_cexp_template.c @@ -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 diff --git a/math/s_clog10_template.c b/math/s_clog10_template.c new file mode 100644 index 0000000000..8d9245bac6 --- /dev/null +++ b/math/s_clog10_template.c @@ -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 diff --git a/math/s_clog_template.c b/math/s_clog_template.c new file mode 100644 index 0000000000..b546030313 --- /dev/null +++ b/math/s_clog_template.c @@ -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 diff --git a/math/s_cpow_template.c b/math/s_cpow_template.c new file mode 100644 index 0000000000..037e575b1a --- /dev/null +++ b/math/s_cpow_template.c @@ -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 diff --git a/math/s_cproj_template.c b/math/s_cproj_template.c new file mode 100644 index 0000000000..d47f009502 --- /dev/null +++ b/math/s_cproj_template.c @@ -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 diff --git a/math/s_csqrt_template.c b/math/s_csqrt_template.c new file mode 100644 index 0000000000..1f073e7f17 --- /dev/null +++ b/math/s_csqrt_template.c @@ -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 |