This is the mail archive of the
libc-alpha@sources.redhat.com
mailing list for the glibc project.
Is sysdeps/ieee754/dbl-64/e_remainder.c broken?
- From: Stephen L Moshier <moshier at mediaone dot net>
- To: "H . J . Lu" <hjl at lucon dot org>, Carsten Langgaard <carstenl at mips dot com>
- Cc: libc-alpha at sourceware dot cygnus dot com
- Date: Sat, 24 Nov 2001 21:18:03 -0500 (EST)
- Subject: Is sysdeps/ieee754/dbl-64/e_remainder.c broken?
- Reply-to: moshier at moshier dot ne dot mediaone dot net
Yes, it is broken on i386.
According to the libc manual,
`drem (6.5, 2.3)' returns `-0.4',
which is `6.5' minus `6.9'.
but the IBM version of e_remainder.c gives a bogus answer on i386
due to an extra-precise fpu register problem.
Here is an exercise program:
-----
double __ieee754_remainder(double, double);
main()
{
double a;
a = __ieee754_remainder(6.5, 2.3);
printf ("%.16e\n", a);
}
-----
Here is an ugly patch for that particular problem.
* sysdeps/ieee754/dbl-64/e_remainder.c (__ieee754_remainder):
Cure an extra-precise fpu register symptom.
*** e_remainder.c 2001/05/12 14:31:44 1.1
--- e_remainder.c 2001/11/25 02:03:11
***************
*** 42,47 ****
--- 42,48 ----
double __ieee754_remainder(double x, double y)
{
double z,d,xx;
+ volatile double r;
#if 0
double yy;
#endif
*************** double __ieee754_remainder(double x, dou
*** 61,67 ****
if ((kx-0x01500000)<ky) {
z=x/t.x;
v.i[HIGH_HALF]=t.i[HIGH_HALF];
! d=(z+big.x)-big.x;
xx=(x-d*v.x)-d*(t.x-v.x);
if (d-z!=0.5&&d-z!=-0.5) return (xx!=0)?xx:((x>0)?ZERO.x:nZERO.x);
else {
--- 62,70 ----
if ((kx-0x01500000)<ky) {
z=x/t.x;
v.i[HIGH_HALF]=t.i[HIGH_HALF];
! r = z + big.x;
! r = r - big.x;
! d = r;
xx=(x-d*v.x)-d*(t.x-v.x);
if (d-z!=0.5&&d-z!=-0.5) return (xx!=0)?xx:((x>0)?ZERO.x:nZERO.x);
else {
*************** double __ieee754_remainder(double x, dou
*** 83,89 ****
z=u.x*r.x;
w.i[HIGH_HALF]=n+l;
ww.i[HIGH_HALF]=(n1)?n1+l:n1;
! d=(z+big.x)-big.x;
u.x=(u.x-d*w.x)-d*ww.x;
l=(u.i[HIGH_HALF]&0x7ff00000)-nn;
}
--- 86,94 ----
z=u.x*r.x;
w.i[HIGH_HALF]=n+l;
ww.i[HIGH_HALF]=(n1)?n1+l:n1;
! r = z + big.x;
! r = r - big.x;
! d = r;
u.x=(u.x-d*w.x)-d*ww.x;
l=(u.i[HIGH_HALF]&0x7ff00000)-nn;
}
*************** double __ieee754_remainder(double x, dou
*** 91,103 ****
w.i[HIGH_HALF]=n;
ww.i[HIGH_HALF]=n1;
z=u.x*r.x;
! d=(z+big.x)-big.x;
u.x=(u.x-d*w.x)-d*ww.x;
if (ABS(u.x)<0.5*t.x) return (u.x!=0)?u.x:((x>0)?ZERO.x:nZERO.x);
else
if (ABS(u.x)>0.5*t.x) return (d>z)?u.x+t.x:u.x-t.x;
else
! {z=u.x/t.x; d=(z+big.x)-big.x; return ((u.x-d*w.x)-d*ww.x);}
}
} /* (kx<0x7ff00000&&ky<0x7ff00000&&ky>=0x03500000) */
--- 96,116 ----
w.i[HIGH_HALF]=n;
ww.i[HIGH_HALF]=n1;
z=u.x*r.x;
! r = z + big.x;
! r = r - big.x;
! d = r;
u.x=(u.x-d*w.x)-d*ww.x;
if (ABS(u.x)<0.5*t.x) return (u.x!=0)?u.x:((x>0)?ZERO.x:nZERO.x);
else
if (ABS(u.x)>0.5*t.x) return (d>z)?u.x+t.x:u.x-t.x;
else
! {
! z=u.x/t.x;
! r = z + big.x;
! r = r - big.x;
! d = r;
! return ((u.x-d*w.x)-d*ww.x);
! }
}
} /* (kx<0x7ff00000&&ky<0x7ff00000&&ky>=0x03500000) */