This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

[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."


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]