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

["H . J . Lu" <hjl at lucon dot org>]

On Tue, Nov 27, 2001 at 12:31:01AM -0800, Geoff Keating wrote:
> > Date: Tue, 27 Nov 2001 00:00:02 -0800
> > From: "H . J . Lu" <hjl@lucon.org>
> > Cc: moshier@moshier.ne.mediaone.net, libc-alpha@sourceware.cygnus.com
> 
> > On Mon, Nov 26, 2001 at 08:03:04PM -0800, Ulrich Drepper wrote:
> > > Stephen L Moshier <moshier@mediaone.net> writes:
> > > 
> > > > Yes, it is broken on i386.
> > > 
> > > The "fix" is simple: don't use the code on x86.  The code was written
> > > using accurate rounding based on the precision and not the extra-long
> > > internal format.  Uglifying and slowing down of generally used code
> > > for a platform which isn't even using the code is not acceptable.
> > 
> > Ok, the IBM code doesn't work on x86. Do we have "make check" results
> > on none x86 platforms?
> 
> Happily, I do.  On powerpc, all the tests pass except for those
> related to 'conj' in test-ifloat.  I need to look at that, but it's
> likely a compiler bug in the 2.95.2 compiler I'm using.
> 

Please try this on your Linux/PPC with the current glibc. I believe
sysdeps/ieee754/dbl-64/e_remainder.c is broken. But I don't know what
is the best way to fix it and verify the whole IBM math library is ok.
I am considering to use the old IEEE math library for MIPS.


H.J.
-----
#include <math.h>
#include <stdio.h>

int main( int argc,char * argv[ ] )
{

  double res;
  
  union {
    unsigned long long l;
    double d;
  } op1, op2;

  op1.l = 0x7fefffffffffffff;
  op2.l = 0x0000000000000001;  

  printf("%llx %llx\n", op1.l, op2.l);
  
  res = remainder(op1.d, op2.d);

  printf("%llx\n", res);

}