This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Clean up dbl-64 rint, nearbyint
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: libc-alpha at sourceware dot org
- Cc: Richard Henderson <rth at twiddle dot net>
- Date: Tue, 13 Mar 2012 00:45:11 +0000 (UTC)
- Subject: Clean up dbl-64 rint, nearbyint
Thomas's comments about applying fixes globally to all instances
reminded me of an item on my glibc todo list to apply my
rintf/nearbyintf fixes to other instances of those functions. This
patch does so for the dbl-64 versions (ldbl-128 appears to have been
the only other affected version, and was fixed by David Miller because
there turned out to be an actual bug in the unnecessary code).
<http://sourceware.org/ml/libc-alpha/2012-02/msg00397.html> explains
the general issue here; I don't know any actual bug fixed by this
cleanup, but it still seems desirable.
Tested x86_64 --disable-multi-arch (directly for the wordsize-64
versions; the non-wordsize-64 versions were tested on x86_64 by
copying them to the wordsize-64 directory so they got used).
Richard, looking at other versions of these functions I'm suspicious
of the alpha versions of nearbyint and nearbyintf (which just add and
subtract 0x1p23 or 0x1p52 without checking the size of the argument).
At least for normal IEEE arithmetic, that's not correct (for double,
consider an argument of 0x1p52 + 1, for example, which is an integer
but whose value will be changed when 0x1p52 is added and subtracted).
If the alpha versions do indeed get such cases wrong I guess it's
evidence we need more testcases for rint and nearbyint to cover these
cases of large integer arguments.
2012-03-12 Joseph Myers <joseph@codesourcery.com>
* sysdeps/ieee754/dbl-64/s_nearbyint.c (__nearbyint): Do not
manipulate bits before adding and subtracting TWO52[sx].
* sysdeps/ieee754/dbl-64/s_rint.c (__rint): Likewise.
* sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c (__nearbyint):
Likewise.
* sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c (__rint): Likewise.
diff --git a/sysdeps/ieee754/dbl-64/s_nearbyint.c b/sysdeps/ieee754/dbl-64/s_nearbyint.c
index d780150..60afafd 100644
--- a/sysdeps/ieee754/dbl-64/s_nearbyint.c
+++ b/sysdeps/ieee754/dbl-64/s_nearbyint.c
@@ -38,18 +38,12 @@ double __nearbyint(double x)
{
fenv_t env;
int32_t i0,j0,sx;
- u_int32_t i,i1;
double w,t;
- EXTRACT_WORDS(i0,i1,x);
+ GET_HIGH_WORD(i0,x);
sx = (i0>>31)&1;
j0 = ((i0>>20)&0x7ff)-0x3ff;
- if(j0<20) {
+ if(j0<52) {
if(j0<0) {
- if(((i0&0x7fffffff)|i1)==0) return x;
- i1 |= (i0&0x0fffff);
- i0 &= 0xfffe0000;
- i0 |= ((i1|-i1)>>12)&0x80000;
- SET_HIGH_WORD(x,i0);
feholdexcept (&env);
w = TWO52[sx]+x;
t = w-TWO52[sx];
@@ -57,32 +51,11 @@ double __nearbyint(double x)
GET_HIGH_WORD(i0,t);
SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
return t;
- } else {
- i = (0x000fffff)>>j0;
- if(((i0&i)|i1)==0) return x; /* x is integral */
- i>>=1;
- if(((i0&i)|i1)!=0) {
- if (j0==19)
- i1 = 0x40000000;
- else if (j0<18)
- i0 = (i0&(~i))|((0x20000)>>j0);
- else
- {
- i0 &= ~i;
- i1 = 0x80000000;
- }
- }
}
- } else if (j0>51) {
+ } else {
if(j0==0x400) return x+x; /* inf or NaN */
else return x; /* x is integral */
- } else {
- i = ((u_int32_t)(0xffffffff))>>(j0-20);
- if((i1&i)==0) return x; /* x is integral */
- i>>=1;
- if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
}
- INSERT_WORDS(x,i0,i1);
feholdexcept (&env);
w = TWO52[sx]+x;
t = w-TWO52[sx];
diff --git a/sysdeps/ieee754/dbl-64/s_rint.c b/sysdeps/ieee754/dbl-64/s_rint.c
index c45d784..8458909 100644
--- a/sysdeps/ieee754/dbl-64/s_rint.c
+++ b/sysdeps/ieee754/dbl-64/s_rint.c
@@ -33,49 +33,22 @@ double
__rint(double x)
{
int32_t i0,j0,sx;
- u_int32_t i,i1;
double w,t;
- EXTRACT_WORDS(i0,i1,x);
+ GET_HIGH_WORD(i0,x);
sx = (i0>>31)&1;
j0 = ((i0>>20)&0x7ff)-0x3ff;
- if(j0<20) {
+ if(j0<52) {
if(j0<0) {
- if(((i0&0x7fffffff)|i1)==0) return x;
- i1 |= (i0&0x0fffff);
- i0 &= 0xfffe0000;
- i0 |= ((i1|-i1)>>12)&0x80000;
- SET_HIGH_WORD(x,i0);
w = TWO52[sx]+x;
t = w-TWO52[sx];
GET_HIGH_WORD(i0,t);
SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
return t;
- } else {
- i = (0x000fffff)>>j0;
- if(((i0&i)|i1)==0) return x; /* x is integral */
- i>>=1;
- if(((i0&i)|i1)!=0) {
- if (j0==19)
- i1 = 0x40000000;
- else if (j0<18)
- i0 = (i0&(~i))|((0x20000)>>j0);
- else
- {
- i0 &= ~i;
- i1 = 0x80000000;
- }
- }
}
- } else if (j0>51) {
+ } else {
if(j0==0x400) return x+x; /* inf or NaN */
else return x; /* x is integral */
- } else {
- i = ((u_int32_t)(0xffffffff))>>(j0-20);
- if((i1&i)==0) return x; /* x is integral */
- i>>=1;
- if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
}
- INSERT_WORDS(x,i0,i1);
w = TWO52[sx]+x;
return w-TWO52[sx];
}
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c
index a380998..a58a620 100644
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c
+++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_nearbyint.c
@@ -41,29 +41,17 @@ __nearbyint(double x)
j0 = ((i0>>52)&0x7ff)-0x3ff;
if(__builtin_expect(j0<52, 1)) {
if(j0<0) {
- if((i0&UINT64_C(0x7fffffffffffffff))==0) return x;
- uint64_t i = i0 & UINT64_C(0xfffffffffffff);
- i0 &= UINT64_C(0xfffe000000000000);
- i0 |= (((i|-i) >> 12) & UINT64_C(0x8000000000000));
- INSERT_WORDS64(x,i0);
libc_feholdexcept (&env);
double w = TWO52[sx]+x;
double t = w-TWO52[sx];
math_opt_barrier(t);
libc_fesetenv (&env);
return copysign(t, x);
- } else {
- uint64_t i = UINT64_C(0x000fffffffffffff)>>j0;
- if((i0&i)==0) return x; /* x is integral */
- i>>=1;
- if((i0&i)!=0)
- i0 = (i0&(~i))|(UINT64_C(0x4000000000000)>>j0);
}
} else {
if(j0==0x400) return x+x; /* inf or NaN */
else return x; /* x is integral */
}
- INSERT_WORDS64(x,i0);
libc_feholdexcept (&env);
double w = TWO52[sx]+x;
double t = w-TWO52[sx];
diff --git a/sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c b/sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c
index 20eef3c..87b2339 100644
--- a/sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c
+++ b/sysdeps/ieee754/dbl-64/wordsize-64/s_rint.c
@@ -38,28 +38,16 @@ __rint(double x)
j0 = ((i0>>52)&0x7ff)-0x3ff;
if(j0<52) {
if(j0<0) {
- if((i0 & UINT64_C(0x7fffffffffffffff))==0) return x;
- uint64_t i = i0 & UINT64_C(0xfffffffffffff);
- i0 &= UINT64_C(0xfffe000000000000);
- i0 |= (((i|-i) >> 12) & UINT64_C(0x8000000000000));
- INSERT_WORDS64(x,i0);
double w = TWO52[sx]+x;
double t = w-TWO52[sx];
EXTRACT_WORDS64(i0,t);
INSERT_WORDS64(t,(i0&UINT64_C(0x7fffffffffffffff))|(sx<<63));
return t;
- } else {
- uint64_t i = UINT64_C(0x000fffffffffffff)>>j0;
- if((i0&i)==0) return x; /* x is integral */
- i>>=1;
- if((i0&i)!=0)
- i0 = (i0&(~i))|(UINT64_C(0x4000000000000)>>j0);
}
} else {
if(j0==0x400) return x+x; /* inf or NaN */
else return x; /* x is integral */
}
- INSERT_WORDS64(x,i0);
double w = TWO52[sx]+x;
return w-TWO52[sx];
}
--
Joseph S. Myers
joseph@codesourcery.com