Fix csqrt overflow/underflow (bug 13841).

8 files changed, 163 insertions, 7 deletions

diff --git a/ChangeLog b/ChangeLog
index 2915e96ad9..41355aeec7 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 ea56e6d832..2f38ad0307 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 191f35959b..caa3ba414d 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 76585e889c..002ea5fdc2 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 d9949c685b..6539ba2249 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 0c624c7a73..64332f67b2 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 977a3abd14..2e86ff6234 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 867e8dd41f..33efe9df51 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