Fix csqrt overflow/underflow (bug 13841).

This commit is contained in:
Joseph Myers 2012-03-14 11:53:32 +00:00
parent 6278569b6c
commit e456826d7a
8 changed files with 163 additions and 7 deletions

View File

@ -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.

2
NEWS
View File

@ -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:

View File

@ -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);
}

View File

@ -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);
}

View File

@ -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);
}

View File

@ -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);
}

View File

@ -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

View File

@ -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,9 +2062,19 @@ 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
idouble: 1