aboutsummaryrefslogtreecommitdiff
path: root/math/s_csqrt.c
diff options
context:
space:
mode:
authorJoseph Myers <joseph@codesourcery.com>2012-07-05 11:02:13 +0000
committerJoseph Myers <joseph@codesourcery.com>2012-07-05 11:02:13 +0000
commitcdfe2c5eb3703ed964cbfdb6906b21ace2956385 (patch)
treec9e8f3452c8450524666f1287f5d8e0fc72a2875 /math/s_csqrt.c
parent704bc4594dc1fad46831823627749fa10924b41d (diff)
downloadglibc-cdfe2c5eb3703ed964cbfdb6906b21ace2956385.tar
glibc-cdfe2c5eb3703ed964cbfdb6906b21ace2956385.tar.gz
glibc-cdfe2c5eb3703ed964cbfdb6906b21ace2956385.tar.bz2
glibc-cdfe2c5eb3703ed964cbfdb6906b21ace2956385.zip
Fix csqrt underflow (bugs 14157, 14331).
Diffstat (limited to 'math/s_csqrt.c')
-rw-r--r--math/s_csqrt.c26
1 files changed, 19 insertions, 7 deletions
diff --git a/math/s_csqrt.c b/math/s_csqrt.c
index 002ea5fdc2..f4d0f998e5 100644
--- a/math/s_csqrt.c
+++ b/math/s_csqrt.c
@@ -75,7 +75,11 @@ __csqrt (__complex__ double x)
}
else if (__builtin_expect (rcls == FP_ZERO, 0))
{
- double r = __ieee754_sqrt (0.5 * fabs (__imag__ x));
+ 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);
@@ -85,13 +89,21 @@ __csqrt (__complex__ double x)
double d, r, s;
int scale = 0;
- if (fabs (__real__ x) > DBL_MAX / 2.0
- || fabs (__imag__ x) > DBL_MAX / 2.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) < DBL_MIN
&& fabs (__imag__ x) < DBL_MIN)
{
@@ -105,13 +117,13 @@ __csqrt (__complex__ double x)
to avoid cancellation error in d +/- Re x. */
if (__real__ x > 0)
{
- r = __ieee754_sqrt (0.5 * d + 0.5 * __real__ x);
- s = (0.5 * __imag__ x) / r;
+ r = __ieee754_sqrt (0.5 * (d + __real__ x));
+ s = 0.5 * (__imag__ x / r);
}
else
{
- s = __ieee754_sqrt (0.5 * d - 0.5 * __real__ x);
- r = fabs ((0.5 * __imag__ x) / s);
+ s = __ieee754_sqrt (0.5 * (d - __real__ x));
+ r = fabs (0.5 * (__imag__ x / s));
}
if (scale)