This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
[PATCH] Fix loss of precision in cosh and sinh for IBM long double
- From: Andreas Schwab <schwab at linux-m68k dot org>
- To: libc-alpha at sourceware dot org
- Date: Mon, 05 Mar 2012 20:40:15 +0100
- Subject: [PATCH] Fix loss of precision in cosh and sinh for IBM long double
For IBM long double exp(-22) is still significant enough that dropping
the term causes a major loss of precision. Let's use 40 as the cutoff
point instead (like the 128 bit long double implementation).
Andreas.
2012-03-05 Andreas Schwab <schwab@linux-m68k.org>
* sysdeps/powerpc/fpu/libm-test-ulps: Update.
* sysdeps/ieee754/ldbl-128ibm/e_coshl.c: Drop exp(-x) term
only for |x| >= 40.
* sysdeps/ieee754/ldbl-128ibm/e_sinhl.c: Likewise.
diff --git c/sysdeps/ieee754/ldbl-128ibm/e_coshl.c w/sysdeps/ieee754/ldbl-128ibm/e_coshl.c
index ebc9436..569b841 100644
--- c/sysdeps/ieee754/ldbl-128ibm/e_coshl.c
+++ w/sysdeps/ieee754/ldbl-128ibm/e_coshl.c
@@ -20,9 +20,9 @@
* 2*exp(x)
*
* exp(x) + 1/exp(x)
- * ln2/2 <= x <= 22 : cosh(x) := -------------------
+ * ln2/2 <= x <= 40 : cosh(x) := -------------------
* 2
- * 22 <= x <= lnovft : cosh(x) := exp(x)/2
+ * 40 <= x <= lnovft : cosh(x) := exp(x)/2
* lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2)
* ln2ovft < x : cosh(x) := huge*huge (overflow)
*
@@ -57,13 +57,13 @@ __ieee754_coshl (long double x)
return one+(t*t)/(w+w);
}
- /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
- if (ix < 0x4036000000000000LL) {
+ /* |x| in [0.5*ln2,40], return (exp(|x|)+1/exp(|x|)/2; */
+ if (ix < 0x4044000000000000LL) {
t = __ieee754_expl(fabsl(x));
return half*t+half/t;
}
- /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
+ /* |x| in [40, log(maxdouble)] return half*exp(|x|) */
if (ix < 0x40862e42fefa39efLL) return half*__ieee754_expl(fabsl(x));
/* |x| in [log(maxdouble), overflowthresold] */
diff --git c/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c w/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c
index b8e86b1..4b53d90 100644
--- c/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c
+++ w/sysdeps/ieee754/ldbl-128ibm/e_sinhl.c
@@ -16,10 +16,10 @@
* 1. Replace x by |x| (sinh(-x) = -sinh(x)).
* 2.
* E + E/(E+1)
- * 0 <= x <= 22 : sinh(x) := --------------, E=expm1(x)
+ * 0 <= x <= 40 : sinh(x) := --------------, E=expm1(x)
* 2
*
- * 22 <= x <= lnovft : sinh(x) := exp(x)/2
+ * 40 <= x <= lnovft : sinh(x) := exp(x)/2
* lnovft <= x <= ln2ovft: sinh(x) := exp(x/2)/2 * exp(x/2)
* ln2ovft < x : sinh(x) := x*shuge (overflow)
*
@@ -48,8 +48,8 @@ __ieee754_sinhl(long double x)
h = 0.5;
if (jx<0) h = -h;
- /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
- if (ix < 0x4036000000000000LL) { /* |x|<22 */
+ /* |x| in [0,40], return sign(x)*0.5*(E+E/(E+1))) */
+ if (ix < 0x4044000000000000LL) { /* |x|<40 */
if (ix<0x3e20000000000000LL) /* |x|<2**-29 */
if(shuge+x>one) return x;/* sinhl(tiny) = tiny with inexact */
t = __expm1l(fabsl(x));
@@ -57,7 +57,7 @@ __ieee754_sinhl(long double x)
return h*(t+t/(t+one));
}
- /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
+ /* |x| in [40, log(maxdouble)] return 0.5*exp(|x|) */
if (ix < 0x40862e42fefa39efLL) return h*__ieee754_expl(fabsl(x));
/* |x| in [log(maxdouble), overflowthresold] */
diff --git c/sysdeps/powerpc/fpu/libm-test-ulps w/sysdeps/powerpc/fpu/libm-test-ulps
index b3aa7bc..f3d4dd9 100644
--- c/sysdeps/powerpc/fpu/libm-test-ulps
+++ w/sysdeps/powerpc/fpu/libm-test-ulps
@@ -431,6 +431,56 @@ Test "cos_upward (9) == -0.9111302618846769883682947111811653112463":
float: 2
ifloat: 2
+# cosh_downward
+Test "cosh_downward (22) == 1792456423.065795780980053377632656584997":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "cosh_downward (23) == 4872401723.124451300068625740569997090344":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "cosh_downward (24) == 13244561064.92173614708845674912733665919":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+# cosh_tonearest
+Test "cosh_tonearest (24) == 13244561064.92173614708845674912733665919":
+ildouble: 1
+ldouble: 1
+
+# cosh_towardzero
+Test "cosh_towardzero (22) == 1792456423.065795780980053377632656584997":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "cosh_towardzero (23) == 4872401723.124451300068625740569997090344":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "cosh_towardzero (24) == 13244561064.92173614708845674912733665919":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+# cosh_upward
+Test "cosh_upward (22) == 1792456423.065795780980053377632656584997":
+ildouble: 2
+ldouble: 2
+Test "cosh_upward (23) == 4872401723.124451300068625740569997090344":
+ildouble: 2
+ldouble: 2
+Test "cosh_upward (24) == 13244561064.92173614708845674912733665919":
+ildouble: 2
+ldouble: 2
+
# cpow
Test "Real part of: cpow (0.75 + 1.25 i, 0.0 + 1.0 i) == 0.331825439177608832276067945276730566 + 0.131338600281188544930936345230903032 i":
float: 1
@@ -934,6 +984,30 @@ Test "log2 (e) == M_LOG2El":
ildouble: 1
ldouble: 1
+# pow_downward
+Test "pow_downward (1.0625, 1.125) == 1.070582293028761362162622578677070098674":
+ildouble: 1
+ldouble: 1
+Test "pow_downward (1.5, 1.03125) == 1.519127098714743184071644334163037684948":
+float: 1
+ifloat: 1
+
+# pow_towardzero
+Test "pow_towardzero (1.0625, 1.125) == 1.070582293028761362162622578677070098674":
+ildouble: 1
+ldouble: 1
+Test "pow_towardzero (1.5, 1.03125) == 1.519127098714743184071644334163037684948":
+float: 1
+ifloat: 1
+
+# pow_upward
+Test "pow_upward (1.0625, 1.125) == 1.070582293028761362162622578677070098674":
+float: 1
+ifloat: 1
+Test "pow_upward (1.5, 1.03125) == 1.519127098714743184071644334163037684948":
+ildouble: 1
+ldouble: 1
+
# sin
Test "sin (16.0) == -0.2879033166650652947844562482186175296207":
ildouble: 2
@@ -1053,6 +1127,44 @@ Test "sinh (0.75) == 0.822316731935829980703661634446913849":
ildouble: 1
ldouble: 1
+# sinh_downward
+Test "sinh_downward (22) == 1792456423.065795780701106568345764104225":
+float: 1
+ifloat: 1
+ildouble: 2
+ldouble: 2
+Test "sinh_downward (23) == 4872401723.124451299966006944252978187305":
+float: 1
+ifloat: 1
+Test "sinh_downward (24) == 13244561064.92173614705070540368454568168":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+# sinh_towardzero
+Test "sinh_towardzero (22) == 1792456423.065795780701106568345764104225":
+float: 1
+ifloat: 1
+ildouble: 2
+ldouble: 2
+Test "sinh_towardzero (23) == 4872401723.124451299966006944252978187305":
+float: 1
+ifloat: 1
+Test "sinh_towardzero (24) == 13244561064.92173614705070540368454568168":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+# sinh_upward
+Test "sinh_upward (23) == 4872401723.124451299966006944252978187305":
+ildouble: 1
+ldouble: 1
+Test "sinh_upward (24) == 13244561064.92173614705070540368454568168":
+ildouble: 1
+ldouble: 1
+
# sqrt
Test "sqrt (0.75) == 0.866025403784438646763723170752936183":
double: 1
@@ -1569,6 +1681,26 @@ Function: "cosh":
ildouble: 1
ldouble: 1
+Function: "cosh_downward":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+Function: "cosh_tonearest":
+ildouble: 1
+ldouble: 1
+
+Function: "cosh_towardzero":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+Function: "cosh_upward":
+ildouble: 2
+ldouble: 2
+
Function: Real part of "cpow":
double: 2
float: 5
@@ -1777,6 +1909,24 @@ Function: "pow":
ildouble: 1
ldouble: 1
+Function: "pow_downward":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+Function: "pow_towardzero":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
+Function: "pow_upward":
+float: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+
Function: "sin":
ildouble: 1
ldouble: 1
@@ -1817,6 +1967,26 @@ Function: "sinh":
ildouble: 1
ldouble: 1
+Function: "sinh_downward":
+float: 1
+ifloat: 1
+ildouble: 2
+ldouble: 2
+
+Function: "sinh_tonearest":
+ildouble: 1
+ldouble: 1
+
+Function: "sinh_towardzero":
+float: 1
+ifloat: 1
+ildouble: 2
+ldouble: 2
+
+Function: "sinh_upward":
+ildouble: 1
+ldouble: 1
+
Function: "sqrt":
double: 1
idouble: 1
--
Andreas Schwab, schwab@linux-m68k.org
GPG Key fingerprint = 58CA 54C7 6D53 942B 1756 01D3 44D5 214B 8276 4ED5
"And now for something completely different."