This is the mail archive of the libc-alpha@sources.redhat.com 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]

Re: Is sysdeps/ieee754/dbl-64/e_remainder.c broken?


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




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