This is the mail archive of the
libc-alpha@sources.redhat.com
mailing list for the glibc project.
Re: Is sysdeps/ieee754/dbl-64/e_remainder.c broken?
- From: Carsten Langgaard <carstenl at mips dot com>
- To: moshier at moshier dot ne dot mediaone dot net
- Cc: "H . J . Lu" <hjl at lucon dot org>, libc-alpha at sourceware dot cygnus dot com
- Date: Mon, 26 Nov 2001 15:54:07 +0100
- Subject: Re: Is sysdeps/ieee754/dbl-64/e_remainder.c broken?
- References: <Pine.LNX.4.20.0111242113280.13575-100000@moshier.ne.mediaone.net>
Stephen L Moshier wrote:
> Yes, it is broken on i386.
I also think it is broken for the rest of the world (MIPS, Sparc, etc....),
but it's a different problem though.
>
> 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) */
--
_ _ ____ ___ Carsten Langgaard Mailto:carstenl@mips.com
|\ /|||___)(___ MIPS Denmark Direct: +45 4486 5527
| \/ ||| ____) Lautrupvang 4B Switch: +45 4486 5555
TECHNOLOGIES 2750 Ballerup Fax...: +45 4486 5556
Denmark http://www.mips.com