Improve clog, clog10 handling of values with real or imaginary part 1 (bug 13629).

7 files changed, 191 insertions, 54 deletions

diff --git a/math/libm-test.inc b/math/libm-test.inc
index 60abf0e47f..c2fb50cec2 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -2464,6 +2464,26 @@ clog_test (void)
 # endif
 #endif
 
+  TEST_c_c (clog, 1.0L, 0x1.234566p-10L, 6.172834701221959432440126967147726538097e-7L, 1.111110564353742042376451655136933182201e-3L);
+  TEST_c_c (clog, -1.0L, 0x1.234566p-20L, 5.886877547844618300918562490463748605537e-13L, 3.141591568520436206990380699322226378452L);
+  TEST_c_c (clog, 0x1.234566p-30L, 1.0L, 5.614163921211322622623353961365728040115e-19L, 1.570796325735258575254858696548386439740L);
+  TEST_c_c (clog, -0x1.234566p-40L, -1.0L, 5.354083939753840089583620652120903838944e-25L, -1.570796326795931422008642456283782656359L);
+  TEST_c_c (clog, 0x1.234566p-50L, 1.0L, 5.106052341226425256332038420428899201070e-31L, 1.570796326794895608681734464330528755366L);
+  TEST_c_c (clog, 0x1.234566p-60L, 1.0L, 4.869510976053643471080816669875627875933e-37L, 1.570796326794896618244456860363082279319L);
+  TEST_c_c (clog, 0x1p-62L, 1.0L, 2.350988701644575015937473074444491355582e-38L, 1.570796326794896619014481257142650555297L);
+  TEST_c_c (clog, 0x1p-63L, 1.0L, 5.877471754111437539843682686111228389059e-39L, 1.570796326794896619122901474391200998698L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (clog, 0x1p-64L, 1.0L, 1.469367938527859384960920671527807097271e-39L, 1.570796326794896619177111583015476220398L, UNDERFLOW_EXCEPTION_FLOAT);
+#ifndef TEST_FLOAT
+  TEST_c_c (clog, 0x1p-510L, 1.0L, 4.450147717014402766180465434664808128438e-308L, 1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION_LDOUBLE_IBM);
+  TEST_c_c (clog, 0x1p-511L, 1.0L, 1.112536929253600691545116358666202032110e-308L, 1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (clog, 0x1p-512L, 1.0L, 2.781342323134001728862790896665505080274e-309L, 1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION_DOUBLE);
+#endif
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+  TEST_c_c (clog, 0x1p-8190L, 1.0L, 6.724206286224187012525355634643505205196e-4932L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (clog, 0x1p-8191L, 1.0L, 1.681051571556046753131338908660876301299e-4932L, 1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (clog, 0x1p-8192L, 1.0L, 4.202628928890116882828347271652190753248e-4933L, 1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION);
+#endif
+
   END (clog, complex);
 }
 
@@ -2593,6 +2613,26 @@ clog10_test (void)
 # endif
 #endif
 
+  TEST_c_c (clog10, 1.0L, 0x1.234566p-10L, 2.680828048441605163181684680300513080769e-7L, 4.825491868832381486767558728169977751564e-4L);
+  TEST_c_c (clog10, -1.0L, 0x1.234566p-20L, 2.556638434669064077889576526006849923281e-13L, 1.364375882602207106407956770293808181427L);
+  TEST_c_c (clog10, 0x1.234566p-30L, 1.0L, 2.438200411482400072282924063740535840474e-19L, 6.821881764607257184291586401763604544928e-1L);
+  TEST_c_c (clog10, -0x1.234566p-40L, -1.0L, 2.325249110681915353442924915876654139373e-25L, -6.821881769213700828789403802671540158935e-1L);
+  TEST_c_c (clog10, 0x1.234566p-50L, 1.0L, 2.217530356103816369479108963807448194409e-31L, 6.821881769209202348667823902864283966959e-1L);
+  TEST_c_c (clog10, 0x1.234566p-60L, 1.0L, 2.114801746467415208319767917450504756866e-37L, 6.821881769209206733143018621078368211515e-1L);
+  TEST_c_c (clog10, 0x1p-61L, 1.0L, 4.084085680564517578238994467153626207224e-38L, 6.821881769209206735545466044044889962925e-1L);
+  TEST_c_c (clog10, 0x1p-62L, 1.0L, 1.021021420141129394559748616788406551878e-38L, 6.821881769209206736487192085600834406988e-1L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (clog10, 0x1p-63L, 1.0L, 2.552553550352823486399371541971016379740e-39L, 6.821881769209206736958055106378806629019e-1L, UNDERFLOW_EXCEPTION_FLOAT);
+#ifndef TEST_FLOAT
+  TEST_c_c (clog10, 0x1p-509L, 1.0L, 7.730698388614835910296270976605350994446e-308L, 6.821881769209206737428918127156778851051e-1L, UNDERFLOW_EXCEPTION_LDOUBLE_IBM);
+  TEST_c_c (clog10, 0x1p-510L, 1.0L, 1.932674597153708977574067744151337748612e-308L, 6.821881769209206737428918127156778851051e-1L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (clog10, 0x1p-511L, 1.0L, 4.831686492884272443935169360378344371529e-309L, 6.821881769209206737428918127156778851051e-1L, UNDERFLOW_EXCEPTION_DOUBLE);
+#endif
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+  TEST_c_c (clog10, 0x1p-8189L, 1.0L, 1.168114274114528946314738738025008370069e-4931L, 6.821881769209206737428918127156778851051e-1L);
+  TEST_c_c (clog10, 0x1p-8190L, 1.0L, 2.920285685286322365786846845062520925172e-4932L, 6.821881769209206737428918127156778851051e-1L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (clog10, 0x1p-8191L, 1.0L, 7.300714213215805914467117112656302312931e-4933L, 6.821881769209206737428918127156778851051e-1L, UNDERFLOW_EXCEPTION);
+#endif
+
   END (clog10, complex);
 }
 
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
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
diff --git a/math/s_clog10f.c b/math/s_clog10f.c
index 5a0c370de5..eb1b895aec 100644
--- a/math/s_clog10f.c
+++ b/math/s_clog10f.c
@@ -44,21 +44,21 @@ __clog10f (__complex__ float x)
     {
       /* Neither real nor imaginary part is NaN.  */
       float absx = fabsf (__real__ x), absy = fabsf (__imag__ x);
-      float d;
       int scale = 0;
 
+      if (absx < absy)
+	{
+	  float t = absx;
+	  absx = absy;
+	  absy = t;
+	}
+
       if (absx > FLT_MAX / 2.0f)
 	{
 	  scale = -1;
 	  absx = __scalbnf (absx, scale);
 	  absy = (absy >= FLT_MIN * 2.0f ? __scalbnf (absy, scale) : 0.0f);
 	}
-      else if (absy > FLT_MAX / 2.0f)
-	{
-	  scale = -1;
-	  absx = (absx >= FLT_MIN * 2.0f ? __scalbnf (absx, scale) : 0.0f);
-	  absy = __scalbnf (absy, scale);
-	}
       else if (absx < FLT_MIN && absy < FLT_MIN)
 	{
 	  scale = FLT_MANT_DIG;
@@ -66,9 +66,29 @@ __clog10f (__complex__ float x)
 	  absy = __scalbnf (absy, scale);
 	}
 
-      d = __ieee754_hypotf (absx, absy);
+      if (absx == 1.0f && scale == 0)
+	{
+	  float absy2 = absy * absy;
+	  if (absy2 <= FLT_MIN * 2.0f * (float) M_LN10)
+	    {
+#if __FLT_EVAL_METHOD__ == 0
+	      __real__ result
+		= (absy2 / 2.0f - absy2 * absy2 / 4.0f) * (float) M_LOG10E;
+#else
+	      volatile float force_underflow = absy2 * absy2 / 4.0f;
+	      __real__ result
+		= (absy2 / 2.0f - force_underflow) * (float) M_LOG10E;
+#endif
+	    }
+	  else
+	    __real__ result = __log1pf (absy2) * ((float) M_LOG10E / 2.0f);
+	}
+      else
+	{
+	  float d = __ieee754_hypotf (absx, absy);
+	  __real__ result = __ieee754_log10f (d) - scale * M_LOG10_2f;
+	}
 
-      __real__ result = __ieee754_log10f (d) - scale * M_LOG10_2f;
       __imag__ result = M_LOG10E * __ieee754_atan2f (__imag__ x, __real__ x);
     }
   else
diff --git a/math/s_clog10l.c b/math/s_clog10l.c
index 3e0ceaa2a8..2a380f6f3b 100644
--- a/math/s_clog10l.c
+++ b/math/s_clog10l.c
@@ -44,21 +44,21 @@ __clog10l (__complex__ long double x)
     {
       /* Neither real nor imaginary part is NaN.  */
       long double absx = fabsl (__real__ x), absy = fabsl (__imag__ x);
-      long double d;
       int scale = 0;
 
+      if (absx < absy)
+	{
+	  long double t = absx;
+	  absx = absy;
+	  absy = t;
+	}
+
       if (absx > LDBL_MAX / 2.0L)
 	{
 	  scale = -1;
 	  absx = __scalbnl (absx, scale);
 	  absy = (absy >= LDBL_MIN * 2.0L ? __scalbnl (absy, scale) : 0.0L);
 	}
-      else if (absy > LDBL_MAX / 2.0L)
-	{
-	  scale = -1;
-	  absx = (absx >= LDBL_MIN * 2.0L ? __scalbnl (absx, scale) : 0.0L);
-	  absy = __scalbnl (absy, scale);
-	}
       else if (absx < LDBL_MIN && absy < LDBL_MIN)
 	{
 	  scale = LDBL_MANT_DIG;
@@ -66,9 +66,21 @@ __clog10l (__complex__ long double x)
 	  absy = __scalbnl (absy, scale);
 	}
 
-      d = __ieee754_hypotl (absx, absy);
+      if (absx == 1.0L && scale == 0)
+	{
+	  long double absy2 = absy * absy;
+	  if (absy2 <= LDBL_MIN * 2.0L * M_LN10l)
+	    __real__ result
+	      = (absy2 / 2.0L - absy2 * absy2 / 4.0L) * M_LOG10El;
+	  else
+	    __real__ result = __log1pl (absy2) * (M_LOG10El / 2.0L);
+	}
+      else
+	{
+	  long double d = __ieee754_hypotl (absx, absy);
+	  __real__ result = __ieee754_log10l (d) - scale * M_LOG10_2l;
+	}
 
-      __real__ result = __ieee754_log10l (d) - scale * M_LOG10_2l;
       __imag__ result = M_LOG10El * __ieee754_atan2l (__imag__ x, __real__ x);
     }
   else
diff --git a/math/s_clogf.c b/math/s_clogf.c
index e20e594964..088730cfd2 100644
--- a/math/s_clogf.c
+++ b/math/s_clogf.c
@@ -41,21 +41,21 @@ __clogf (__complex__ float x)
     {
       /* Neither real nor imaginary part is NaN.  */
       float absx = fabsf (__real__ x), absy = fabsf (__imag__ x);
-      float d;
       int scale = 0;
 
+      if (absx < absy)
+	{
+	  float t = absx;
+	  absx = absy;
+	  absy = t;
+	}
+
       if (absx > FLT_MAX / 2.0f)
 	{
 	  scale = -1;
 	  absx = __scalbnf (absx, scale);
 	  absy = (absy >= FLT_MIN * 2.0f ? __scalbnf (absy, scale) : 0.0f);
 	}
-      else if (absy > FLT_MAX / 2.0f)
-	{
-	  scale = -1;
-	  absx = (absx >= FLT_MIN * 2.0f ? __scalbnf (absx, scale) : 0.0f);
-	  absy = __scalbnf (absy, scale);
-	}
       else if (absx < FLT_MIN && absy < FLT_MIN)
 	{
 	  scale = FLT_MANT_DIG;
@@ -63,9 +63,27 @@ __clogf (__complex__ float x)
 	  absy = __scalbnf (absy, scale);
 	}
 
-      d = __ieee754_hypotf (absx, absy);
+      if (absx == 1.0f && scale == 0)
+	{
+	  float absy2 = absy * absy;
+	  if (absy2 <= FLT_MIN * 2.0f)
+	    {
+#if __FLT_EVAL_METHOD__ == 0
+	      __real__ result = absy2 / 2.0f - absy2 * absy2 / 4.0f;
+#else
+	      volatile float force_underflow = absy2 * absy2 / 4.0f;
+	      __real__ result = absy2 / 2.0f - force_underflow;
+#endif
+	    }
+	  else
+	    __real__ result = __log1pf (absy2) / 2.0f;
+	}
+      else
+	{
+	  float d = __ieee754_hypotf (absx, absy);
+	  __real__ result = __ieee754_logf (d) - scale * (float) M_LN2;
+	}
 
-      __real__ result = __ieee754_logf (d) - scale * (float) M_LN2;
       __imag__ result = __ieee754_atan2f (__imag__ x, __real__ x);
     }
   else
diff --git a/math/s_clogl.c b/math/s_clogl.c
index a393b88a5f..d98a3c02c1 100644
--- a/math/s_clogl.c
+++ b/math/s_clogl.c
@@ -41,21 +41,21 @@ __clogl (__complex__ long double x)
     {
       /* Neither real nor imaginary part is NaN.  */
       long double absx = fabsl (__real__ x), absy = fabsl (__imag__ x);
-      long double d;
       int scale = 0;
 
+      if (absx < absy)
+	{
+	  long double t = absx;
+	  absx = absy;
+	  absy = t;
+	}
+
       if (absx > LDBL_MAX / 2.0L)
 	{
 	  scale = -1;
 	  absx = __scalbnl (absx, scale);
 	  absy = (absy >= LDBL_MIN * 2.0L ? __scalbnl (absy, scale) : 0.0L);
 	}
-      else if (absy > LDBL_MAX / 2.0L)
-	{
-	  scale = -1;
-	  absx = (absx >= LDBL_MIN * 2.0L ? __scalbnl (absx, scale) : 0.0L);
-	  absy = __scalbnl (absy, scale);
-	}
       else if (absx < LDBL_MIN && absy < LDBL_MIN)
 	{
 	  scale = LDBL_MANT_DIG;
@@ -63,9 +63,20 @@ __clogl (__complex__ long double x)
 	  absy = __scalbnl (absy, scale);
 	}
 
-      d = __ieee754_hypotl (absx, absy);
+      if (absx == 1.0L && scale == 0)
+	{
+	  long double absy2 = absy * absy;
+	  if (absy2 <= LDBL_MIN * 2.0L)
+	    __real__ result = absy2 / 2.0L - absy2 * absy2 / 4.0L;
+	  else
+	    __real__ result = __log1pl (absy2) / 2.0L;
+	}
+      else
+	{
+	  long double d = __ieee754_hypotl (absx, absy);
+	  __real__ result = __ieee754_logl (d) - scale * M_LN2l;
+	}
 
-      __real__ result = __ieee754_logl (d) - scale * M_LN2l;
       __imag__ result = __ieee754_atan2l (__imag__ x, __real__ x);
     }
   else