diff options
Diffstat (limited to 'math/s_clog.c')
| -rw-r--r-- | math/s_clog.c | 36 |
1 files changed, 27 insertions, 9 deletions
diff --git a/math/s_clog.c b/math/s_clog.c index 2249e86b43..e28aa5157d 100644 --- a/math/s_clog.c +++ b/math/s_clog.c @@ -41,21 +41,21 @@ __clog (__complex__ double x) { /* Neither real nor imaginary part is NaN. */ double absx = fabs (__real__ x), absy = fabs (__imag__ x); - double d; int scale = 0; + if (absx < absy) + { + double t = absx; + absx = absy; + absy = t; + } + if (absx > DBL_MAX / 2.0) { scale = -1; absx = __scalbn (absx, scale); absy = (absy >= DBL_MIN * 2.0 ? __scalbn (absy, scale) : 0.0); } - else if (absy > DBL_MAX / 2.0) - { - scale = -1; - absx = (absx >= DBL_MIN * 2.0 ? __scalbn (absx, scale) : 0.0); - absy = __scalbn (absy, scale); - } else if (absx < DBL_MIN && absy < DBL_MIN) { scale = DBL_MANT_DIG; @@ -63,9 +63,27 @@ __clog (__complex__ double x) absy = __scalbn (absy, scale); } - d = __ieee754_hypot (absx, absy); + if (absx == 1.0 && scale == 0) + { + double absy2 = absy * absy; + if (absy2 <= DBL_MIN * 2.0) + { +#if __FLT_EVAL_METHOD__ == 0 + __real__ result = absy2 / 2.0 - absy2 * absy2 / 4.0; +#else + volatile double force_underflow = absy2 * absy2 / 4.0; + __real__ result = absy2 / 2.0 - force_underflow; +#endif + } + else + __real__ result = __log1p (absy2) / 2.0; + } + else + { + double d = __ieee754_hypot (absx, absy); + __real__ result = __ieee754_log (d) - scale * M_LN2; + } - __real__ result = __ieee754_log (d) - scale * M_LN2; __imag__ result = __ieee754_atan2 (__imag__ x, __real__ x); } else |