This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix ldbl-128 cos range reduction near pi/2 (bug 15429)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: <libc-alpha at sourceware dot org>
- Date: Thu, 9 May 2013 19:23:55 +0000
- Subject: Fix ldbl-128 cos range reduction near pi/2 (bug 15429)
Bug 15429 is inaccuracy of ldbl-128 cos near pi/2. The cause is
inaccurate range reduction in that area; I don't know why the code
uses only 93 bits for the high part of pi/2, but you need at least 226
bits to get reasonable accuracy there, and subtraction of a 113-bit
high part is always going to be exact in the relevant range. This
patch changes the code to use 113 bits for both high and low parts.
Tested mips64 to confirm the large ulps no longer appear.
Notes:
(a) The low part y[1] of the reduced value may not be very accurate -
for that, you'd need 339 bits - but the low part isn't particularly
significant to the final result so this shouldn't increase errors by
more than 1ulp.
(b) I'm testing a similar patch for the ldbl-128ibm case (bug 15359).
2013-05-09 Joseph Myers <joseph@codesourcery.com>
[BZ #15429]
* sysdeps/ieee754/ldbl-128/e_rem_pio2l.c (c): Use 113 bits for
high part of pi/2.
(__ieee754_rem_pio2l): Update comments.
diff --git a/sysdeps/ieee754/ldbl-128/e_rem_pio2l.c b/sysdeps/ieee754/ldbl-128/e_rem_pio2l.c
index 84846fd..ee856ac 100644
--- a/sysdeps/ieee754/ldbl-128/e_rem_pio2l.c
+++ b/sysdeps/ieee754/ldbl-128/e_rem_pio2l.c
@@ -184,13 +184,13 @@ static const int32_t two_over_pi[] = {
};
static const long double c[] = {
-/* 93 bits of pi/2 */
+/* 113 bits of pi/2 */
#define PI_2_1 c[0]
- 1.57079632679489661923132169155131424e+00L, /* 3fff921fb54442d18469898cc5100000 */
+ 0x1.921fb54442d18469898cc51701b8p+0L,
/* pi/2 - PI_2_1 */
#define PI_2_1t c[1]
- 8.84372056613570112025531863263659260e-29L, /* 3fa1c06e0e68948127044533e63a0106 */
+ 0x3.9a252049c1114cf98e804177d4c8p-116L,
};
int32_t __ieee754_rem_pio2l(long double x, long double *y)
@@ -213,7 +213,7 @@ int32_t __ieee754_rem_pio2l(long double x, long double *y)
{
if (hx > 0)
{
- /* 113 + 93 bit PI is ok */
+ /* 113 + 113 bit PI is ok */
z = x - PI_2_1;
y[0] = z - PI_2_1t;
y[1] = (z - y[0]) - PI_2_1t;
@@ -221,7 +221,7 @@ int32_t __ieee754_rem_pio2l(long double x, long double *y)
}
else
{
- /* 113 + 93 bit PI is ok */
+ /* 113 + 113 bit PI is ok */
z = x + PI_2_1;
y[0] = z + PI_2_1t;
y[1] = (z - y[0]) + PI_2_1t;
--
Joseph S. Myers
joseph@codesourcery.com