This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix csqrt spurious underflows (bug 18371) [committed]
- From: Joseph Myers <joseph at codesourcery dot com>
- To: <libc-alpha at sourceware dot org>
- Date: Tue, 23 Jun 2015 16:02:25 +0000
- Subject: Fix csqrt spurious underflows (bug 18371) [committed]
- Authentication-results: sourceware.org; auth=none
The csqrt implementations in glibc can cause spurious underflows in
some cases as a side-effect of the scaling for large arguments (when
underflow is correct for the square root of the argument that was
scaled down to avoid overflow, but not for the original argument).
This patch arranges to avoid the underflowing intermediate computation
(eliminating a multiplication in 0.5 in the problem cases where a
subsequent scaling by 2 would follow).
Tested for x86_64 and x86 and ulps updated accordingly (only needed
for x86). Committed.
(auto-libm-test-out diffs omitted below.)
2015-06-23 Joseph Myers <joseph@codesourcery.com>
[BZ #18371]
* math/s_csqrt.c (__csqrt): Avoid multiplication by 0.5 where
intermediate but not final result might underflow.
* math/s_csqrtf.c (__csqrtf): Likewise.
* math/s_csqrtl.c (__csqrtl): Likewise.
* math/auto-libm-test-in: Add more tests of csqrt.
* math/auto-libm-test-out: Regenerated.
* sysdeps/i386/fpu/libm-test-ulps: Update.
diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index cf7fc35..d1502ad 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -1033,6 +1033,13 @@ csqrt 0x1p-16494 0x1p+16383
csqrt 0x1.0000000000000000000000000001p-16382 0x1.0000000000000000000000000001p-16382
csqrt -0x1.0000000000000000000000000001p-16382 -0x1.0000000000000000000000000001p-16382
+csqrt 0x0.ffp128 0x1.1p-61
+csqrt -0x0.ffp128 0x1.1p-61
+csqrt 0x0.ffp1024 0x1.1p-509
+csqrt -0x0.ffp1024 0x1.1p-509
+csqrt 0x0.ffp16384 0x1.1p-8189
+csqrt -0x0.ffp16384 0x1.1p-8189
+
ctan 0 0
ctan 0 -0
ctan -0 0
diff --git a/math/s_csqrt.c b/math/s_csqrt.c
index 66505af..75da2e7 100644
--- a/math/s_csqrt.c
+++ b/math/s_csqrt.c
@@ -118,12 +118,28 @@ __csqrt (__complex__ double x)
if (__real__ x > 0)
{
r = __ieee754_sqrt (0.5 * (d + __real__ x));
- s = 0.5 * (__imag__ x / r);
+ 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));
- r = fabs (0.5 * (__imag__ x / s));
+ 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)
diff --git a/math/s_csqrtf.c b/math/s_csqrtf.c
index c9a800e..7f45cc1 100644
--- a/math/s_csqrtf.c
+++ b/math/s_csqrtf.c
@@ -118,12 +118,28 @@ __csqrtf (__complex__ float x)
if (__real__ x > 0)
{
r = __ieee754_sqrtf (0.5f * (d + __real__ x));
- s = 0.5f * (__imag__ x / r);
+ if (scale == 1 && fabsf (__imag__ x) < 1.0f)
+ {
+ /* Avoid possible intermediate underflow. */
+ s = __imag__ x / r;
+ r = __scalbnf (r, scale);
+ scale = 0;
+ }
+ else
+ s = 0.5f * (__imag__ x / r);
}
else
{
s = __ieee754_sqrtf (0.5f * (d - __real__ x));
- r = fabsf (0.5f * (__imag__ x / s));
+ if (scale == 1 && fabsf (__imag__ x) < 1.0f)
+ {
+ /* Avoid possible intermediate underflow. */
+ r = fabsf (__imag__ x / s);
+ s = __scalbnf (s, scale);
+ scale = 0;
+ }
+ else
+ r = fabsf (0.5f * (__imag__ x / s));
}
if (scale)
diff --git a/math/s_csqrtl.c b/math/s_csqrtl.c
index de23f8d..10988ba 100644
--- a/math/s_csqrtl.c
+++ b/math/s_csqrtl.c
@@ -118,12 +118,28 @@ __csqrtl (__complex__ long double x)
if (__real__ x > 0)
{
r = __ieee754_sqrtl (0.5L * (d + __real__ x));
- s = 0.5L * (__imag__ x / r);
+ if (scale == 1 && fabsl (__imag__ x) < 1.0L)
+ {
+ /* Avoid possible intermediate underflow. */
+ s = __imag__ x / r;
+ r = __scalbnl (r, scale);
+ scale = 0;
+ }
+ else
+ s = 0.5L * (__imag__ x / r);
}
else
{
s = __ieee754_sqrtl (0.5L * (d - __real__ x));
- r = fabsl (0.5L * (__imag__ x / s));
+ if (scale == 1 && fabsl (__imag__ x) < 1.0L)
+ {
+ /* Avoid possible intermediate underflow. */
+ r = fabsl (__imag__ x / s);
+ s = __scalbnl (s, scale);
+ scale = 0;
+ }
+ else
+ r = fabsl (0.5L * (__imag__ x / s));
}
if (scale)
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 3716a55..f319749 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -1108,6 +1108,8 @@ ildouble: 2
ldouble: 2
Function: Imaginary part of "csqrt":
+double: 1
+idouble: 1
ildouble: 2
ldouble: 2
--
Joseph S. Myers
joseph@codesourcery.com