Fix csqrt overflow/underflow (bug 13841).

diff --git a/ChangeLog b/ChangeLog
index 2915e96ad9182478502b2ffd76cc9d9b9911d8da..41355aeec746121781be1dc046a37f67233e0345 100644 (file)
--- 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 ea56e6d832e2a5d254d10ab3ba47732e59905784..2f38ad0307a47cccc53f19584b11d68c0bfe6fb0 100644 (file)
--- 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 191f35959b43625ad89556ca546e3512cdd231d2..caa3ba414d5e87600e0d974ec8d02f909b172fbc 100644 (file)
--- 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 76585e889c932e1ef400e6bb347b0154cc3965fa..002ea5fdc2e20f6da3da68747c216489a7e697f7 100644 (file)
--- 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 d9949c685b0c2b8c3338ae622d7296e00f39b485..6539ba2249aa124d35b3610d530723028e04a31d 100644 (file)
--- 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 0c624c7a7379b1edcc562f07bb8b6bea48e076ae..64332f67b2a84d966e2056cd1383954056f11bdb 100644 (file)
--- 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 977a3abd141b3e7ee33547a408c64036e12448ed..2e86ff6234a3a30668b8ea9ea7d4f31632b0a750 100644 (file)
--- 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 867e8dd41fd82b1e2e38d1d3bcfb8dd99fc3d424..33efe9df51c539b3c7105da3ec8e120a39999340 100644 (file)
--- 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