aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorJoseph Myers <joseph@codesourcery.com>2013-03-16 17:50:28 +0000
committerJoseph Myers <joseph@codesourcery.com>2013-03-16 17:50:28 +0000
commit2a185d32e830589bf9ae50f9243bb304f84b110b (patch)
tree0b002b84729b6ebdb44b04b93849b4d960e3fae5
parent6cbec759de7941016b30a5e46bdef535657ed0eb (diff)
downloadglibc-2a185d32e830589bf9ae50f9243bb304f84b110b.tar
glibc-2a185d32e830589bf9ae50f9243bb304f84b110b.tar.gz
glibc-2a185d32e830589bf9ae50f9243bb304f84b110b.tar.bz2
glibc-2a185d32e830589bf9ae50f9243bb304f84b110b.zip
Fix spurious underflow exceptions for Bessel functions for ldbl-128 / ldbl-128ibm (bug 14155).
-rw-r--r--ChangeLog13
-rw-r--r--math/libm-test.inc12
-rw-r--r--sysdeps/ieee754/ldbl-128/e_j0l.c68
-rw-r--r--sysdeps/ieee754/ldbl-128/e_j1l.c69
4 files changed, 96 insertions, 66 deletions
diff --git a/ChangeLog b/ChangeLog
index e931f27c41..8192d5ae54 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,16 @@
+2013-03-16 Joseph Myers <joseph@codesourcery.com>
+
+ [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 <thomas@codesourcery.com>
* 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 c85bdcb35f..d9df034c63 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 1b18289588..9e7880c49d 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 f16343b26b..95e01a39cc 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;
}