diff options
| author | Joseph Myers <joseph@codesourcery.com> | 2012-07-05 11:02:13 +0000 |
|---|---|---|
| committer | Joseph Myers <joseph@codesourcery.com> | 2012-07-05 11:02:13 +0000 |
| commit | cdfe2c5eb3703ed964cbfdb6906b21ace2956385 (patch) | |
| tree | c9e8f3452c8450524666f1287f5d8e0fc72a2875 /math/s_csqrt.c | |
| parent | 704bc4594dc1fad46831823627749fa10924b41d (diff) | |
| download | glibc-cdfe2c5eb3703ed964cbfdb6906b21ace2956385.tar.gz glibc-cdfe2c5eb3703ed964cbfdb6906b21ace2956385.tar.xz glibc-cdfe2c5eb3703ed964cbfdb6906b21ace2956385.zip | |
Fix csqrt underflow (bugs 14157, 14331).
Diffstat (limited to 'math/s_csqrt.c')
| -rw-r--r-- | math/s_csqrt.c | 26 |
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) |