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]

k_sinl bug fix for ldbl-128 and ldbl-128ibm


This *untested* patch fixes the problem I noticed while preparing the
ldbl-96 sinl/cosl patch.  Some code in __kernel_sinl, which is trying
to compute the sine of x+y, replaces x by fabsl (x), then changes the
x+y pair to an h+l pair (where h has a precomputed sine and cosine,
and |l| <= 1/256.0).  If x was negated, this should be allowed for
when y is included in the computation of the new low part l - but it
isn't allowed for, meaning l is off by |2y|.

This won't make a difference of more than an ulp or two; the evidence
that it is a useful fix is that the corresponding ldbl-96 fix reduced
the number of ulps updates needed in that patch on x86 from 30 to 20.

Comments?  I can do some testing if needed, but some other people may
be more conveniently set up to test ldbl-128 and ldbl-128ibm than I
am.  (It is of course possible that in particular cases a more
accurate internal calculation leads to a less accurate final result
and so some ulps baseline changes may be needed.)

2012-03-16  Joseph Myers  <joseph@codesourcery.com>

	* sysdeps/ieee754/ldbl-128/k_sinl.c (__kernel_sinl): Negate y when
	computing low part if x was negated.
	* sysdeps/ieee754/ldbl-128ibm/k_sinl.c (__kernel_sinl): Likewise.

diff --git a/sysdeps/ieee754/ldbl-128/k_sinl.c b/sysdeps/ieee754/ldbl-128/k_sinl.c
index 6eaf878..1f0ca8c 100644
--- a/sysdeps/ieee754/ldbl-128/k_sinl.c
+++ b/sysdeps/ieee754/ldbl-128/k_sinl.c
@@ -116,7 +116,7 @@ __kernel_sinl(long double x, long double y, int iy)
 
       SET_LDOUBLE_WORDS64(h, ((u_int64_t)hix) << 32, 0);
       if (iy)
-	l = y - (h - x);
+	l = (ix < 0 ? -y : y) - (h - x);
       else
 	l = x - h;
       z = l * l;
diff --git a/sysdeps/ieee754/ldbl-128ibm/k_sinl.c b/sysdeps/ieee754/ldbl-128ibm/k_sinl.c
index 484b65f..94aba0b 100644
--- a/sysdeps/ieee754/ldbl-128ibm/k_sinl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/k_sinl.c
@@ -134,7 +134,7 @@ __kernel_sinl(long double x, long double y, int iy)
 */
       SET_LDOUBLE_WORDS64(h, ((u_int64_t)hix) << 32, 0);
       if (iy)
-	l = y - (h - x);
+	l = (ix < 0 ? -y : y) - (h - x);
       else
 	l = x - h;
       z = l * l;

-- 
Joseph S. Myers
joseph@codesourcery.com


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