From: Joseph Myers Date: Wed, 14 Mar 2012 11:53:32 +0000 (+0000) Subject: Fix csqrt overflow/underflow (bug 13841). X-Git-Tag: upstream/2.30~11445 X-Git-Url: http://review.tizen.org/git/?a=commitdiff_plain;h=e456826d7a539fb322bb9719297bd386eded8e32;p=external%2Fglibc.git Fix csqrt overflow/underflow (bug 13841). --- diff --git a/ChangeLog b/ChangeLog index 2915e96..41355ae 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,14 @@ 2012-03-14 Joseph Myers + [BZ #13841] + * math/s_csqrt.c: Include . + (__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 . Contributed by Ulrich Drepper , 1997. @@ -21,7 +21,7 @@ #include #include #include - +#include __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 . Contributed by Ulrich Drepper , 1997. @@ -21,7 +21,7 @@ #include #include #include - +#include __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 . Contributed by Ulrich Drepper , 1997. @@ -21,7 +21,7 @@ #include #include #include - +#include __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