This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Remove ldbl-128ibm variants of complex math functions [committed]
- From: Joseph Myers <joseph at codesourcery dot com>
- To: <libc-alpha at sourceware dot org>
- Date: Wed, 17 Jun 2015 21:21:09 +0000
- Subject: Remove ldbl-128ibm variants of complex math functions [committed]
- Authentication-results: sourceware.org; auth=none
sysdeps/ieee754/ldbl-128ibm has its own versions of cprojl, ctanhl and
ctanl.
Having its own versions, where otherwise the math/ copies are
generally used for all floating-point formats, means they are liable
to get out of sync and not benefit from bug fixes to the generic
versions. The substantive differences (not arising from getting out
of sync and slightly different fixes for the same issues) are: long
double compat handling (also done in the ldbl-opt versions, so doesn't
require special versions for ldbl-128ibm); handling of LDBL_EPSILON
(conditionally undefined and redefined in other math/ implementations,
so doesn't justify a special version), and:
/* __gcc_qmul does not respect -0.0 so we need the following fixup. */
if ((__real__ res == 0.0L) && (__real__ x == 0.0L))
__real__ res = __real__ x;
if ((__real__ res == 0.0L) && (__imag__ x == 0.0L))
__imag__ res = __imag__ x;
But if that statement about __gcc_qmul was ever true for an old
version of that libgcc function, it's not the case for any GCC version
now supported to build glibc; there's explicit logic early in that
function (and similarly in __gcc_qdiv) to return an appropriately
signed zero if the product of the high parts is zero. So this patch
adds the special LDBL_EPSILON handling to the generic functions and
removes the ldbl-128ibm versions.
Tested for powerpc32 (compared test-ldouble.out before and after the
changes; there are slight changes to results for ctanl / ctanhl,
arising from divergence of the implementations, but nothing that
affects the overall nature of the issues shown by the testsuite, and
in particular nothing related to signs of zero results). Committed.
2015-06-17 Joseph Myers <joseph@codesourcery.com>
* math/s_ctanhl.c [LDBL_MANT_DIG == 106] (LDBL_EPSILON): Undefine
and redefine.
* math/s_ctanl.c [LDBL_MANT_DIG == 106] (LDBL_EPSILON): Undefine
and redefine.
* sysdeps/ieee754/ldbl-128ibm/s_cprojl.c: Remove file.
* sysdeps/ieee754/ldbl-128ibm/s_ctanhl.c: Likewise.
* sysdeps/ieee754/ldbl-128ibm/s_ctanl.c: Likewise.
diff --git a/math/s_ctanhl.c b/math/s_ctanhl.c
index 016ebc5..fa11db7 100644
--- a/math/s_ctanhl.c
+++ b/math/s_ctanhl.c
@@ -23,6 +23,13 @@
#include <math_private.h>
#include <float.h>
+/* To avoid spurious underflows, use this definition to treat IBM long
+ double as approximating an IEEE-style format. */
+#if LDBL_MANT_DIG == 106
+# undef LDBL_EPSILON
+# define LDBL_EPSILON 0x1p-106L
+#endif
+
__complex__ long double
__ctanhl (__complex__ long double x)
{
diff --git a/math/s_ctanl.c b/math/s_ctanl.c
index 6ad1b75..1f92755 100644
--- a/math/s_ctanl.c
+++ b/math/s_ctanl.c
@@ -23,6 +23,13 @@
#include <math_private.h>
#include <float.h>
+/* To avoid spurious underflows, use this definition to treat IBM long
+ double as approximating an IEEE-style format. */
+#if LDBL_MANT_DIG == 106
+# undef LDBL_EPSILON
+# define LDBL_EPSILON 0x1p-106L
+#endif
+
__complex__ long double
__ctanl (__complex__ long double x)
{
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_cprojl.c b/sysdeps/ieee754/ldbl-128ibm/s_cprojl.c
deleted file mode 100644
index 4243a4b..0000000
--- a/sysdeps/ieee754/ldbl-128ibm/s_cprojl.c
+++ /dev/null
@@ -1,39 +0,0 @@
-/* Compute projection of complex long double value to Riemann sphere.
- Copyright (C) 1997-2015 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_ldbl_opt.h>
-
-__complex__ long double
-__cprojl (__complex__ long double x)
-{
- if (__isinf_nsl (__real__ x) || __isinf_nsl (__imag__ x))
- {
- __complex__ long double res;
-
- __real__ res = INFINITY;
- __imag__ res = __copysignl (0.0, __imag__ x);
-
- return res;
- }
-
- return x;
-}
-long_double_symbol (libm, __cprojl, cprojl);
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_ctanhl.c b/sysdeps/ieee754/ldbl-128ibm/s_ctanhl.c
deleted file mode 100644
index 96bdccd..0000000
--- a/sysdeps/ieee754/ldbl-128ibm/s_ctanhl.c
+++ /dev/null
@@ -1,121 +0,0 @@
-/* Complex hyperbole tangent for long double. IBM extended format version.
- Copyright (C) 1997-2015 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 <float.h>
-#include <math.h>
-#include <math_ldbl_opt.h>
-
-#include <math_private.h>
-
-/* IBM long double GCC builtin sets LDBL_EPSILON == LDBL_DENORM_MIN */
-static const long double ldbl_eps = 0x1p-106L;
-
-__complex__ long double
-__ctanhl (__complex__ long double x)
-{
- __complex__ long double res;
-
- if (!isfinite (__real__ x) || !isfinite (__imag__ x))
- {
- if (isinf (__real__ x))
- {
- __real__ res = __copysignl (1.0L, __real__ x);
- __imag__ res = __copysignl (0.0L, __imag__ x);
- }
- else if (__imag__ x == 0.0)
- {
- res = x;
- }
- else
- {
- __real__ res = __nanl ("");
- __imag__ res = __nanl ("");
-
-#ifdef FE_INVALID
- if (isinf (__imag__ x))
- feraiseexcept (FE_INVALID);
-#endif
- }
- }
- else
- {
- long double sinix, cosix;
- long double den;
- const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l / 2.0L);
-
- /* tanh(x+iy) = (sinh(2x) + i*sin(2y))/(cosh(2x) + cos(2y))
- = (sinh(x)*cosh(x) + i*sin(y)*cos(y))/(sinh(x)^2 + cos(y)^2). */
-
- __sincosl (__imag__ x, &sinix, &cosix);
-
- if (fabsl (__real__ x) > t)
- {
- /* Avoid intermediate overflow when the imaginary part of
- the result may be subnormal. Ignoring negligible terms,
- the real part is +/- 1, the imaginary part is
- sin(y)*cos(y)/sinh(x)^2 = 4*sin(y)*cos(y)/exp(2x). */
- long double exp_2t = __ieee754_expl (2 * t);
- __real__ res = __copysignl (1.0L, __real__ x);
- __imag__ res = 4 * sinix * cosix;
- __real__ x = fabsl (__real__ x);
- __real__ x -= t;
- __imag__ res /= exp_2t;
- if (__real__ x > t)
- {
- /* Underflow (original real part of x has absolute value
- > 2t). */
- __imag__ res /= exp_2t;
- }
- else
- __imag__ res /= __ieee754_expl (2.0L * __real__ x);
- }
- else
- {
- long double sinhrx, coshrx;
- if (fabs (__real__ x) > LDBL_MIN)
- {
- sinhrx = __ieee754_sinhl (__real__ x);
- coshrx = __ieee754_coshl (__real__ x);
- }
- else
- {
- sinhrx = __real__ x;
- coshrx = 1.0L;
- }
-
- if (fabsl (sinhrx) > fabsl (cosix) * ldbl_eps)
- den = sinhrx * sinhrx + cosix * cosix;
- else
- den = cosix * cosix;
- __real__ res = sinhrx * (coshrx / den);
- __imag__ res = sinix * (cosix / den);
- }
- /* __gcc_qmul does not respect -0.0 so we need the following fixup. */
- if ((__real__ res == 0.0L) && (__real__ x == 0.0L))
- __real__ res = __real__ x;
-
- if ((__real__ res == 0.0L) && (__imag__ x == 0.0L))
- __imag__ res = __imag__ x;
- }
-
- return res;
-}
-long_double_symbol (libm, __ctanhl, ctanhl);
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_ctanl.c b/sysdeps/ieee754/ldbl-128ibm/s_ctanl.c
deleted file mode 100644
index fc44378..0000000
--- a/sysdeps/ieee754/ldbl-128ibm/s_ctanl.c
+++ /dev/null
@@ -1,121 +0,0 @@
-/* Complex tangent function for long double. IBM extended format version.
- Copyright (C) 1997-2015 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_ldbl_opt.h>
-#include <float.h>
-
-#include <math_private.h>
-
-/* IBM long double GCC builtin sets LDBL_EPSILON == LDBL_DENORM_MIN */
-static const long double ldbl_eps = 0x1p-106L;
-
-__complex__ long double
-__ctanl (__complex__ long double x)
-{
- __complex__ long double res;
-
- if (!isfinite (__real__ x) || !isfinite (__imag__ x))
- {
- if (isinf (__imag__ x))
- {
- __real__ res = __copysignl (0.0, __real__ x);
- __imag__ res = __copysignl (1.0, __imag__ x);
- }
- else if (__real__ x == 0.0)
- {
- res = x;
- }
- else
- {
- __real__ res = __nanl ("");
- __imag__ res = __nanl ("");
-
- if (__isinf_nsl (__real__ x))
- feraiseexcept (FE_INVALID);
- }
- }
- else
- {
- long double sinrx, cosrx;
- long double den;
- const int t = (int) ((LDBL_MAX_EXP - 1) * M_LN2l / 2.0L);
-
- /* tan(x+iy) = (sin(2x) + i*sinh(2y))/(cos(2x) + cosh(2y))
- = (sin(x)*cos(x) + i*sinh(y)*cosh(y)/(cos(x)^2 + sinh(y)^2). */
-
- __sincosl (__real__ x, &sinrx, &cosrx);
-
- if (fabsl (__imag__ x) > t)
- {
- /* Avoid intermediate overflow when the real part of the
- result may be subnormal. Ignoring negligible terms, the
- imaginary part is +/- 1, the real part is
- sin(x)*cos(x)/sinh(y)^2 = 4*sin(x)*cos(x)/exp(2y). */
- long double exp_2t = __ieee754_expl (2 * t);
-
- __imag__ res = __copysignl (1.0L, __imag__ x);
- __real__ res = 4 * sinrx * cosrx;
- __imag__ x = fabsl (__imag__ x);
- __imag__ x -= t;
- __real__ res /= exp_2t;
- if (__imag__ x > t)
- {
- /* Underflow (original imaginary part of x has absolute
- value > 2t). */
- __real__ res /= exp_2t;
- }
- else
- __real__ res /= __ieee754_expl (2.0L * __imag__ x);
- }
- else
- {
- long double sinhix, coshix;
- if (fabsl (__imag__ x) > LDBL_MIN)
- {
- sinhix = __ieee754_sinhl (__imag__ x);
- coshix = __ieee754_coshl (__imag__ x);
- }
- else
- {
- sinhix = __imag__ x;
- coshix = 1.0L;
- }
-
- if (fabsl (sinhix) > fabsl (cosrx) * ldbl_eps)
- den = cosrx * cosrx + sinhix * sinhix;
- else
- den = cosrx * cosrx;
- __real__ res = sinrx * (cosrx / den);
- __imag__ res = sinhix * (coshix / den);
- }
-
- /* __gcc_qmul does not respect -0.0 so we need the following fixup. */
- if ((__real__ res == 0.0L) && (__real__ x == 0.0L))
- __real__ res = __real__ x;
-
- if ((__real__ res == 0.0L) && (__imag__ x == 0.0L))
- __imag__ res = __imag__ x;
- }
-
- return res;
-}
-long_double_symbol (libm, __ctanl, ctanl);
--
Joseph S. Myers
joseph@codesourcery.com