This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix csqrt overflow/underflow (bug 13841)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: libc-alpha at sourceware dot org
- Date: Tue, 13 Mar 2012 21:30:12 +0000 (UTC)
- Subject: Fix csqrt overflow/underflow (bug 13841)
This patch fixes overflow and underflow in csqrt functions. These use
hypot internally, which inevitably means overflows for large inputs
and loss of precision for subnormal inputs, even though the final
correct results of csqrt will not overflow. The fix is to scale the
inputs appropriately before calling hypot, then scale the result
afterwards, in the potentially affected cases.
This patch depends on my hypotf fix
<http://sourceware.org/ml/libc-alpha/2012-03/msg00357.html> (that is,
some new tests will fail without that fix and it's through those
failures that I found the hypotf bug).
Tested x86_64 and x86 and ulps updated based on that testing.
2012-03-13 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.
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 39cda66..b30abe9 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
--
Joseph S. Myers
joseph@codesourcery.com