diff options
Diffstat (limited to 'math/s_clog10.c')
| -rw-r--r-- | math/s_clog10.c | 36 |
1 files changed, 27 insertions, 9 deletions
diff --git a/math/s_clog10.c b/math/s_clog10.c index cf5fb8a0b5..b733b048b9 100644 --- a/math/s_clog10.c +++ b/math/s_clog10.c @@ -44,21 +44,21 @@ __clog10 (__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; @@ -66,9 +66,27 @@ __clog10 (__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 * M_LN10) + { +#if __FLT_EVAL_METHOD__ == 0 + __real__ result = (absy2 / 2.0 - absy2 * absy2 / 4.0) * M_LOG10E; +#else + volatile double force_underflow = absy2 * absy2 / 4.0; + __real__ result = (absy2 / 2.0 - force_underflow) * M_LOG10E; +#endif + } + else + __real__ result = __log1p (absy2) * (M_LOG10E / 2.0); + } + else + { + double d = __ieee754_hypot (absx, absy); + __real__ result = __ieee754_log10 (d) - scale * M_LOG10_2; + } - __real__ result = __ieee754_log10 (d) - scale * M_LOG10_2; __imag__ result = M_LOG10E * __ieee754_atan2 (__imag__ x, __real__ x); } else |