From: Joseph Myers Date: Sat, 16 Mar 2013 17:50:28 +0000 (+0000) Subject: Fix spurious underflow exceptions for Bessel functions for ldbl-128 / ldbl-128ibm... X-Git-Tag: upstream/2.20~2336 X-Git-Url: http://review.tizen.org/git/?a=commitdiff_plain;h=2a185d32e830589bf9ae50f9243bb304f84b110b;p=platform%2Fupstream%2Flinaro-glibc.git Fix spurious underflow exceptions for Bessel functions for ldbl-128 / ldbl-128ibm (bug 14155). --- diff --git a/ChangeLog b/ChangeLog index e931f27..8192d5a 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,16 @@ +2013-03-16 Joseph Myers + + [BZ #14155] + * sysdeps/ieee754/ldbl-128/e_j0l.c (__ieee754_j0l): Do not compute + 1 / x and functions P and Q for arguments above 0x1p256L. + (__ieee754_y0l): Likewise. + * sysdeps/ieee754/ldbl-128/e_j1l.c (__ieee754_j1l): Likewise. + (__ieee754_y1l): Likewise. + * math/libm-test.inc (j0_test): Do not allow spurious underflows. + (j1_test): Likewise. + (y0_test): Likewise. + (y1_test): Likewise. + 2013-03-16 Thomas Schwinge * math/test-snan.c (TEST_FUNC): Add and use minus_sNaN_var diff --git a/math/libm-test.inc b/math/libm-test.inc index c85bdcb..d9df034 100644 --- a/math/libm-test.inc +++ b/math/libm-test.inc @@ -6243,8 +6243,7 @@ j0_test (void) TEST_f_f (j0, 0x1.d7ce3ap+107L, 2.775523647291230802651040996274861694514e-17L); #ifndef TEST_FLOAT - /* Bug 14155: spurious exception may occur. */ - TEST_f_f (j0, -0x1.001000001p+593L, -3.927269966354206207832593635798954916263e-90L, UNDERFLOW_EXCEPTION_OK); + TEST_f_f (j0, -0x1.001000001p+593L, -3.927269966354206207832593635798954916263e-90L); #endif END (j0); @@ -6285,8 +6284,7 @@ j1_test (void) TEST_f_f (j1, 0x1.3ffp+74L, 1.818984347516051243459364437186082741567e-12L); #ifndef TEST_FLOAT - /* Bug 14155: spurious exception may occur. */ - TEST_f_f (j1, 0x1.ff00000000002p+840L, 1.846591691699331493194965158699937660696e-127L, UNDERFLOW_EXCEPTION_OK); + TEST_f_f (j1, 0x1.ff00000000002p+840L, 1.846591691699331493194965158699937660696e-127L); #endif END (j1); @@ -10457,8 +10455,7 @@ y0_test (void) TEST_f_f (y0, 0x1.3ffp+74L, 1.818984347516051243459467456433028748678e-12L); #ifndef TEST_FLOAT - /* Bug 14155: spurious exception may occur. */ - TEST_f_f (y0, 0x1.ff00000000002p+840L, 1.846591691699331493194965158699937660696e-127L, UNDERFLOW_EXCEPTION_OK); + TEST_f_f (y0, 0x1.ff00000000002p+840L, 1.846591691699331493194965158699937660696e-127L); #endif TEST_f_f (y0, 0x1p-10L, -4.4865150767109739412411806297168793661098L); @@ -10511,8 +10508,7 @@ y1_test (void) TEST_f_f (y1, 0x1.27e204p+99L, -8.881610148467797208469612080785210013461e-16L); #ifndef TEST_FLOAT - /* Bug 14155: spurious exception may occur. */ - TEST_f_f (y1, 0x1.001000001p+593L, 3.927269966354206207832593635798954916263e-90L, UNDERFLOW_EXCEPTION_OK); + TEST_f_f (y1, 0x1.001000001p+593L, 3.927269966354206207832593635798954916263e-90L); #endif TEST_f_f (y1, 0x1p-10L, -6.5190099301063115047395187618929589514382e+02L); diff --git a/sysdeps/ieee754/ldbl-128/e_j0l.c b/sysdeps/ieee754/ldbl-128/e_j0l.c index 1b18289..9e7880c 100644 --- a/sysdeps/ieee754/ldbl-128/e_j0l.c +++ b/sysdeps/ieee754/ldbl-128/e_j0l.c @@ -700,6 +700,25 @@ __ieee754_j0l (long double x) return p; } + /* X = x - pi/4 + cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4) + = 1/sqrt(2) * (cos(x) + sin(x)) + sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4) + = 1/sqrt(2) * (sin(x) - cos(x)) + sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) + cf. Fdlibm. */ + __sincosl (xx, &s, &c); + ss = s - c; + cc = s + c; + z = -__cosl (xx + xx); + if ((s * c) < 0) + cc = z / ss; + else + ss = z / cc; + + if (xx > 0x1p256L) + return ONEOSQPI * cc / __ieee754_sqrtl (xx); + xinv = 1.0L / xx; z = xinv * xinv; if (xinv <= 0.25) @@ -761,21 +780,6 @@ __ieee754_j0l (long double x) p = 1.0L + z * p; q = z * xinv * q; q = q - 0.125L * xinv; - /* X = x - pi/4 - cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4) - = 1/sqrt(2) * (cos(x) + sin(x)) - sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4) - = 1/sqrt(2) * (sin(x) - cos(x)) - sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - cf. Fdlibm. */ - __sincosl (xx, &s, &c); - ss = s - c; - cc = s + c; - z = -__cosl (xx + xx); - if ((s * c) < 0) - cc = z / ss; - else - ss = z / cc; z = ONEOSQPI * (p * cc - q * ss) / __ieee754_sqrtl (xx); return z; } @@ -843,6 +847,25 @@ long double return p; } + /* X = x - pi/4 + cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4) + = 1/sqrt(2) * (cos(x) + sin(x)) + sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4) + = 1/sqrt(2) * (sin(x) - cos(x)) + sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) + cf. Fdlibm. */ + __sincosl (x, &s, &c); + ss = s - c; + cc = s + c; + z = -__cosl (x + x); + if ((s * c) < 0) + cc = z / ss; + else + ss = z / cc; + + if (xx > 0x1p256L) + return ONEOSQPI * ss / __ieee754_sqrtl (x); + xinv = 1.0L / xx; z = xinv * xinv; if (xinv <= 0.25) @@ -904,21 +927,6 @@ long double p = 1.0L + z * p; q = z * xinv * q; q = q - 0.125L * xinv; - /* X = x - pi/4 - cos(X) = cos(x) cos(pi/4) + sin(x) sin(pi/4) - = 1/sqrt(2) * (cos(x) + sin(x)) - sin(X) = sin(x) cos(pi/4) - cos(x) sin(pi/4) - = 1/sqrt(2) * (sin(x) - cos(x)) - sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x)) - cf. Fdlibm. */ - __sincosl (x, &s, &c); - ss = s - c; - cc = s + c; - z = -__cosl (x + x); - if ((s * c) < 0) - cc = z / ss; - else - ss = z / cc; z = ONEOSQPI * (p * ss + q * cc) / __ieee754_sqrtl (x); return z; } diff --git a/sysdeps/ieee754/ldbl-128/e_j1l.c b/sysdeps/ieee754/ldbl-128/e_j1l.c index f16343b..95e01a3 100644 --- a/sysdeps/ieee754/ldbl-128/e_j1l.c +++ b/sysdeps/ieee754/ldbl-128/e_j1l.c @@ -706,6 +706,29 @@ __ieee754_j1l (long double x) return p; } + /* X = x - 3 pi/4 + cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4) + = 1/sqrt(2) * (-cos(x) + sin(x)) + sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4) + = -1/sqrt(2) * (sin(x) + cos(x)) + cf. Fdlibm. */ + __sincosl (xx, &s, &c); + ss = -s - c; + cc = s - c; + z = __cosl (xx + xx); + if ((s * c) > 0) + cc = z / ss; + else + ss = z / cc; + + if (xx > 0x1p256L) + { + z = ONEOSQPI * cc / __ieee754_sqrtl (xx); + if (x < 0) + z = -z; + return z; + } + xinv = 1.0L / xx; z = xinv * xinv; if (xinv <= 0.25) @@ -767,20 +790,6 @@ __ieee754_j1l (long double x) p = 1.0L + z * p; q = z * q; q = q * xinv + 0.375L * xinv; - /* X = x - 3 pi/4 - cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4) - = 1/sqrt(2) * (-cos(x) + sin(x)) - sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4) - = -1/sqrt(2) * (sin(x) + cos(x)) - cf. Fdlibm. */ - __sincosl (xx, &s, &c); - ss = -s - c; - cc = s - c; - z = __cosl (xx + xx); - if ((s * c) > 0) - cc = z / ss; - else - ss = z / cc; z = ONEOSQPI * (p * cc - q * ss) / __ieee754_sqrtl (xx); if (x < 0) z = -z; @@ -850,6 +859,24 @@ __ieee754_y1l (long double x) return p; } + /* X = x - 3 pi/4 + cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4) + = 1/sqrt(2) * (-cos(x) + sin(x)) + sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4) + = -1/sqrt(2) * (sin(x) + cos(x)) + cf. Fdlibm. */ + __sincosl (xx, &s, &c); + ss = -s - c; + cc = s - c; + z = __cosl (xx + xx); + if ((s * c) > 0) + cc = z / ss; + else + ss = z / cc; + + if (xx > 0x1p256L) + return ONEOSQPI * ss / __ieee754_sqrtl (xx); + xinv = 1.0L / xx; z = xinv * xinv; if (xinv <= 0.25) @@ -911,20 +938,6 @@ __ieee754_y1l (long double x) p = 1.0L + z * p; q = z * q; q = q * xinv + 0.375L * xinv; - /* X = x - 3 pi/4 - cos(X) = cos(x) cos(3 pi/4) + sin(x) sin(3 pi/4) - = 1/sqrt(2) * (-cos(x) + sin(x)) - sin(X) = sin(x) cos(3 pi/4) - cos(x) sin(3 pi/4) - = -1/sqrt(2) * (sin(x) + cos(x)) - cf. Fdlibm. */ - __sincosl (xx, &s, &c); - ss = -s - c; - cc = s - c; - z = __cosl (xx + xx); - if ((s * c) > 0) - cc = z / ss; - else - ss = z / cc; z = ONEOSQPI * (p * ss + q * cc) / __ieee754_sqrtl (xx); return z; }