Fix csqrt overflow/underflow (bug 13841).

8 files changed, 163 insertions, 7 deletions

diff --git a/ChangeLog b/ChangeLog
index 2915e96..41355ae 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,14 @@
 2012-03-14  Joseph Myers  <joseph@codesourcery.com>
 
+	[BZ #13841]
+	* math/s_csqrt.c: Include <float.h>.
+	(__csqrt): Scale large or subnormal inputs.
+	* math/s_csqrtf.c: Likewise.
+	* math/s_csqrtl.c: Likewise.
+	* math/libm-test.inc (csqrt_test): Add more tests.
+	* sysdeps/i386/fpu/libm-test-ulps: Update.
+	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
+
 	[BZ #13840]
 	* math/libm-test.inc (hypot_test): Add more tests.
 
diff --git a/NEWS b/NEWS
index ea56e6d..2f38ad0 100644
--- a/NEWS
+++ b/NEWS
@@ -14,7 +14,7 @@ Version 2.16
   10210, 10545, 10716, 11174, 11322, 11365, 11494, 12047, 13058, 13525,
   13526, 13527, 13528, 13529, 13530, 13531, 13532, 13533, 13547, 13551,
   13552, 13553, 13555, 13559, 13566, 13583, 13618, 13637, 13656, 13673,
-  13695, 13704, 13706, 13726, 13738, 13786, 13792, 13806, 13840
+  13695, 13704, 13706, 13726, 13738, 13786, 13792, 13806, 13840, 13841
 
 * ISO C11 support:
 
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 191f359..caa3ba4 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -2657,6 +2657,24 @@ csqrt_test (void)
      part).  */
   TEST_c_c (csqrt, 0, -1, M_SQRT_2_2, -M_SQRT_2_2);
 
+  TEST_c_c (csqrt, 0x1.fffffep+127L, 0x1.fffffep+127L, 2.026714405498316804978751017492482558075e+19L, 8.394925938143272988211878516208015586281e+18L);
+  TEST_c_c (csqrt, 0x1.fffffep+127L, 1.0L, 1.844674352395372953599975585936590505260e+19L, 2.710505511993121390769065968615872097053e-20L);
+  TEST_c_c (csqrt, 0x1p-149L, 0x1p-149L, 4.112805464342778798097003462770175200803e-23L, 1.703579802732953750368659735601389709551e-23L);
+  TEST_c_c (csqrt, 0x1p-147L, 0x1p-147L, 8.225610928685557596194006925540350401606e-23L, 3.407159605465907500737319471202779419102e-23L);
+
+#ifndef TEST_FLOAT
+  TEST_c_c (csqrt, 0x1.fffffffffffffp+1023L, 0x1.fffffffffffffp+1023L, 1.473094556905565378990473658199034571917e+154L, 6.101757441282702188537080005372547713595e+153L);
+  TEST_c_c (csqrt, 0x1.fffffffffffffp+1023L, 0x1p+1023L, 1.379778091031440685006200821918878702861e+154L, 3.257214233483129514781233066898042490248e+153L);
+  TEST_c_c (csqrt, 0x1p-1074L, 0x1p-1074L, 2.442109726130830256743814843868934877597e-162L, 1.011554969366634726113090867589031782487e-162L);
+  TEST_c_c (csqrt, 0x1p-1073L, 0x1p-1073L, 3.453664695497464982856905711457966660085e-162L, 1.430554756764195530630723976279903095110e-162L);
+#endif
+
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+  TEST_c_c (csqrt, 0x1.fp+16383L, 0x1.fp+16383L, 1.179514222452201722651836720466795901016e+2466L, 4.885707879516577666702435054303191575148e+2465L);
+  TEST_c_c (csqrt, 0x1.fp+16383L, 0x1p+16383L, 1.106698967236475180613254276996359485630e+2466L, 2.687568007603946993388538156299100955642e+2465L);
+  TEST_c_c (csqrt, 0x1p-16440L, 0x1p-16441L, 3.514690655930285351254618340783294558136e-2475L,  8.297059146828716918029689466551384219370e-2476L);
+#endif
+
   END (csqrt, complex);
 }
 
diff --git a/math/s_csqrt.c b/math/s_csqrt.c
index 76585e8..002ea5f 100644
--- a/math/s_csqrt.c
+++ b/math/s_csqrt.c
@@ -1,5 +1,5 @@
 /* Complex square root of double value.
-   Copyright (C) 1997, 1998, 2005, 2011 Free Software Foundation, Inc.
+   Copyright (C) 1997-2012 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.
@@ -21,7 +21,7 @@
 #include <complex.h>
 #include <math.h>
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ double
 __csqrt (__complex__ double x)
@@ -83,6 +83,22 @@ __csqrt (__complex__ double x)
       else
 	{
 	  double d, r, s;
+	  int scale = 0;
+
+	  if (fabs (__real__ x) > DBL_MAX / 2.0
+	      || fabs (__imag__ x) > DBL_MAX / 2.0)
+	    {
+	      scale = 1;
+	      __real__ x = __scalbn (__real__ x, -2 * scale);
+	      __imag__ x = __scalbn (__imag__ x, -2 * scale);
+	    }
+	  else if (fabs (__real__ x) < DBL_MIN
+		   && fabs (__imag__ x) < DBL_MIN)
+	    {
+	      scale = -(DBL_MANT_DIG / 2);
+	      __real__ x = __scalbn (__real__ x, -2 * scale);
+	      __imag__ x = __scalbn (__imag__ x, -2 * scale);
+	    }
 
 	  d = __ieee754_hypot (__real__ x, __imag__ x);
 	  /* Use the identity   2  Re res  Im res = Im x
@@ -98,6 +114,12 @@ __csqrt (__complex__ double x)
 	      r = fabs ((0.5 * __imag__ x) / s);
 	    }
 
+	  if (scale)
+	    {
+	      r = __scalbn (r, scale);
+	      s = __scalbn (s, scale);
+	    }
+
 	  __real__ res = r;
 	  __imag__ res = __copysign (s, __imag__ x);
 	}
diff --git a/math/s_csqrtf.c b/math/s_csqrtf.c
index d9949c6..6539ba2 100644
--- a/math/s_csqrtf.c
+++ b/math/s_csqrtf.c
@@ -1,5 +1,5 @@
 /* Complex square root of float value.
-   Copyright (C) 1997, 1998, 2005, 2011 Free Software Foundation, Inc.
+   Copyright (C) 1997-2012 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.
@@ -21,7 +21,7 @@
 #include <complex.h>
 #include <math.h>
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ float
 __csqrtf (__complex__ float x)
@@ -83,6 +83,22 @@ __csqrtf (__complex__ float x)
       else
 	{
 	  float d, r, s;
+	  int scale = 0;
+
+	  if (fabsf (__real__ x) > FLT_MAX / 2.0f
+	      || fabsf (__imag__ x) > FLT_MAX / 2.0f)
+	    {
+	      scale = 1;
+	      __real__ x = __scalbnf (__real__ x, -2 * scale);
+	      __imag__ x = __scalbnf (__imag__ x, -2 * scale);
+	    }
+	  else if (fabsf (__real__ x) < FLT_MIN
+		   && fabsf (__imag__ x) < FLT_MIN)
+	    {
+	      scale = -(FLT_MANT_DIG / 2);
+	      __real__ x = __scalbnf (__real__ x, -2 * scale);
+	      __imag__ x = __scalbnf (__imag__ x, -2 * scale);
+	    }
 
 	  d = __ieee754_hypotf (__real__ x, __imag__ x);
 	  /* Use the identity   2  Re res  Im res = Im x
@@ -98,6 +114,12 @@ __csqrtf (__complex__ float x)
 	      r = fabsf ((0.5f * __imag__ x) / s);
 	    }
 
+	  if (scale)
+	    {
+	      r = __scalbnf (r, scale);
+	      s = __scalbnf (s, scale);
+	    }
+
 	  __real__ res = r;
 	  __imag__ res = __copysignf (s, __imag__ x);
 	}
diff --git a/math/s_csqrtl.c b/math/s_csqrtl.c
index 0c624c7..64332f6 100644
--- a/math/s_csqrtl.c
+++ b/math/s_csqrtl.c
@@ -1,5 +1,5 @@
 /* Complex square root of long double value.
-   Copyright (C) 1997, 1998, 2005, 2011 Free Software Foundation, Inc.
+   Copyright (C) 1997-2012 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.
@@ -21,7 +21,7 @@
 #include <complex.h>
 #include <math.h>
 #include <math_private.h>
-
+#include <float.h>
 
 __complex__ long double
 __csqrtl (__complex__ long double x)
@@ -83,6 +83,22 @@ __csqrtl (__complex__ long double x)
       else
 	{
 	  long double d, r, s;
+	  int scale = 0;
+
+	  if (fabsl (__real__ x) > LDBL_MAX / 2.0L
+	      || fabsl (__imag__ x) > LDBL_MAX / 2.0L)
+	    {
+	      scale = 1;
+	      __real__ x = __scalbnl (__real__ x, -2 * scale);
+	      __imag__ x = __scalbnl (__imag__ x, -2 * scale);
+	    }
+	  else if (fabsl (__real__ x) < LDBL_MIN
+		   && fabsl (__imag__ x) < LDBL_MIN)
+	    {
+	      scale = -(LDBL_MANT_DIG / 2);
+	      __real__ x = __scalbnl (__real__ x, -2 * scale);
+	      __imag__ x = __scalbnl (__imag__ x, -2 * scale);
+	    }
 
 	  d = __ieee754_hypotl (__real__ x, __imag__ x);
 	  /* Use the identity   2  Re res  Im res = Im x
@@ -98,6 +114,12 @@ __csqrtl (__complex__ long double x)
 	      r = fabsl ((0.5L * __imag__ x) / s);
 	    }
 
+	  if (scale)
+	    {
+	      r = __scalbnl (r, scale);
+	      s = __scalbnl (s, scale);
+	    }
+
 	  __real__ res = r;
 	  __imag__ res = __copysignl (s, __imag__ x);
 	}
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 977a3ab..2e86ff6 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -805,6 +805,26 @@ Test "Imaginary part of: csinh (0.75 + 1.25 i) == 0.2592948545511627791533498306
 float: 1
 ifloat: 1
 
+# csqrt
+Test "Imaginary part of: csqrt (0x1.fffffffffffffp+1023 + 0x1p+1023 i) == 1.379778091031440685006200821918878702861e+154 + 3.257214233483129514781233066898042490248e+153 i":
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: csqrt (0x1.fp+16383 + 0x1.fp+16383 i) == 1.179514222452201722651836720466795901016e+2466 + 4.885707879516577666702435054303191575148e+2465 i":
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: csqrt (0x1p-1073 + 0x1p-1073 i) == 3.453664695497464982856905711457966660085e-162 + 1.430554756764195530630723976279903095110e-162 i":
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: csqrt (0x1p-1074 + 0x1p-1074 i) == 2.442109726130830256743814843868934877597e-162 + 1.011554969366634726113090867589031782487e-162 i":
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: csqrt (0x1p-147 + 0x1p-147 i) == 8.225610928685557596194006925540350401606e-23 + 3.407159605465907500737319471202779419102e-23 i":
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: csqrt (0x1p-149 + 0x1p-149 i) == 4.112805464342778798097003462770175200803e-23 + 1.703579802732953750368659735601389709551e-23 i":
+ildouble: 1
+ldouble: 1
+
 # ctan
 Test "Real part of: ctan (-2 - 3 i) == 0.376402564150424829275122113032269084e-2 - 1.00323862735360980144635859782192726 i":
 double: 1
@@ -2054,6 +2074,10 @@ ifloat: 1
 ildouble: 2
 ldouble: 2
 
+Function: Imaginary part of "csqrt":
+ildouble: 1
+ldouble: 1
+
 Function: Real part of "ctan":
 double: 1
 idouble: 1
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index 867e8dd..33efe9d 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -848,6 +848,35 @@ ifloat: 1
 Test "Real part of: csqrt (-2 - 3 i) == 0.89597747612983812471573375529004348 - 1.6741492280355400404480393008490519 i":
 float: 1
 ifloat: 1
+Test "Imaginary part of: csqrt (0x1.fffffep+127 + 1.0 i) == 1.844674352395372953599975585936590505260e+19 + 2.710505511993121390769065968615872097053e-20 i":
+float: 1
+ifloat: 1
+Test "Real part of: csqrt (0x1.fffffffffffffp+1023 + 0x1.fffffffffffffp+1023 i) == 1.473094556905565378990473658199034571917e+154 + 6.101757441282702188537080005372547713595e+153 i":
+double: 1
+idouble: 1
+Test "Imaginary part of: csqrt (0x1.fffffffffffffp+1023 + 0x1.fffffffffffffp+1023 i) == 1.473094556905565378990473658199034571917e+154 + 6.101757441282702188537080005372547713595e+153 i":
+double: 1
+idouble: 1
+Test "Imaginary part of: csqrt (0x1.fffffffffffffp+1023 + 0x1p+1023 i) == 1.379778091031440685006200821918878702861e+154 + 3.257214233483129514781233066898042490248e+153 i":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: csqrt (0x1.fp+16383 + 0x1.fp+16383 i) == 1.179514222452201722651836720466795901016e+2466 + 4.885707879516577666702435054303191575148e+2465 i":
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: csqrt (0x1p-1073 + 0x1p-1073 i) == 3.453664695497464982856905711457966660085e-162 + 1.430554756764195530630723976279903095110e-162 i":
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: csqrt (0x1p-1074 + 0x1p-1074 i) == 2.442109726130830256743814843868934877597e-162 + 1.011554969366634726113090867589031782487e-162 i":
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: csqrt (0x1p-147 + 0x1p-147 i) == 8.225610928685557596194006925540350401606e-23 + 3.407159605465907500737319471202779419102e-23 i":
+ildouble: 1
+ldouble: 1
+Test "Imaginary part of: csqrt (0x1p-149 + 0x1p-149 i) == 4.112805464342778798097003462770175200803e-23 + 1.703579802732953750368659735601389709551e-23 i":
+ildouble: 1
+ldouble: 1
 
 # ctan
 Test "Real part of: ctan (-2 - 3 i) == 0.376402564150424829275122113032269084e-2 - 1.00323862735360980144635859782192726 i":
@@ -2033,8 +2062,18 @@ ildouble: 2
 ldouble: 2
 
 Function: Real part of "csqrt":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+
+Function: Imaginary part of "csqrt":
+double: 1
 float: 1
+idouble: 1
 ifloat: 1
+ildouble: 1
+ldouble: 1
 
 Function: Real part of "ctan":
 double: 1