Update copyright notices with scripts/update-copyrights

[thirdparty/glibc.git] / math / s_clog10.c

diff --git a/math/s_clog10.c b/math/s_clog10.c
index cf5fb8a0b5ca77b4b4a89ee5510317c079bc8c0e..0274db36172bbee90b6d2bb0ce294c898595e0e9 100644 (file)
--- a/math/s_clog10.c
+++ b/math/s_clog10.c
@@ -1,5 +1,5 @@
 /* Compute complex base 10 logarithm.
-   Copyright (C) 1997-2012 Free Software Foundation, Inc.
+   Copyright (C) 1997-2014 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
@@ -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,47 @@ __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 if (absx > 1.0 && absx < 2.0 && absy < 1.0 && scale == 0)
+       {
+         double d2m1 = (absx - 1.0) * (absx + 1.0);
+         if (absy >= DBL_EPSILON)
+           d2m1 += absy * absy;
+         __real__ result = __log1p (d2m1) * (M_LOG10E / 2.0);
+       }
+      else if (absx < 1.0
+              && absx >= 0.75
+              && absy < DBL_EPSILON / 2.0
+              && scale == 0)
+       {
+         double d2m1 = (absx - 1.0) * (absx + 1.0);
+         __real__ result = __log1p (d2m1) * (M_LOG10E / 2.0);
+       }
+      else if (absx < 1.0 && (absx >= 0.75 || absy >= 0.5) && scale == 0)
+       {
+         double d2m1 = __x2y2m1 (absx, absy);
+         __real__ result = __log1p (d2m1) * (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