diff options
| author | Samuel Thibault <samuel.thibault@ens-lyon.org> | 2016-08-20 19:50:45 +0200 |
|---|---|---|
| committer | Samuel Thibault <samuel.thibault@ens-lyon.org> | 2016-08-20 19:50:45 +0200 |
| commit | 4dd9e35bfd35d3138bc44169baba098005bad51e (patch) | |
| tree | a4939c43a9c3fe00eb27f023e14acc5e1fe8808c /math/s_csqrt.c | |
| parent | bd42a4599d1b6f77bcfe1e4f67b7cbd9e1cb2dfd (diff) | |
| parent | f76453c31593957fec1a99b986bfa5506618b79c (diff) | |
Merge commit 'refs/top-bases/t/bigmem' into t/bigmem
Diffstat (limited to 'math/s_csqrt.c')
| -rw-r--r-- | math/s_csqrt.c | 28 |
1 files changed, 22 insertions, 6 deletions
diff --git a/math/s_csqrt.c b/math/s_csqrt.c index 2bbd2b8adf..75da2e7449 100644 --- a/math/s_csqrt.c +++ b/math/s_csqrt.c @@ -1,5 +1,5 @@ /* Complex square root of double value. - Copyright (C) 1997-2014 Free Software Foundation, Inc. + Copyright (C) 1997-2015 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. @@ -30,7 +30,7 @@ __csqrt (__complex__ double x) int rcls = fpclassify (__real__ x); int icls = fpclassify (__imag__ x); - if (__builtin_expect (rcls <= FP_INFINITE || icls <= FP_INFINITE, 0)) + if (__glibc_unlikely (rcls <= FP_INFINITE || icls <= FP_INFINITE)) { if (icls == FP_INFINITE) { @@ -59,7 +59,7 @@ __csqrt (__complex__ double x) } else { - if (__builtin_expect (icls == FP_ZERO, 0)) + if (__glibc_unlikely (icls == FP_ZERO)) { if (__real__ x < 0.0) { @@ -73,7 +73,7 @@ __csqrt (__complex__ double x) __imag__ res = __copysign (0.0, __imag__ x); } } - else if (__builtin_expect (rcls == FP_ZERO, 0)) + else if (__glibc_unlikely (rcls == FP_ZERO)) { double r; if (fabs (__imag__ x) >= 2.0 * DBL_MIN) @@ -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) |