incorrectly rounded square root

Paul Zimmermann <[email protected]>
Newsgroups gmane.comp.lib.newlib
Message-ID <[email protected]>
       Hi,

according to IEEE 754, the square root function should be correctly rounded
for all rounding modes. I noticed this is not the case in Newlib:

$ cat test_sqrt.c
#include <stdio.h>
#include <math.h>
#include <fenv.h>

#ifdef NEWLIB
int errno;
int* __errno () { return &errno; }
#endif

int main()
{
  int rnd[4] = { FE_TONEAREST, FE_TOWARDZERO, FE_UPWARD, FE_DOWNWARD };
  char Rnd[4] = "NZUD";
  float x = 0x1.ff07fep+127f;
  float y;
  for (int i = 0; i < 4; i++)
  {
    fesetround (rnd[i]);
    y = sqrtf (x);
    printf ("RND%c: %a\n", Rnd[i], y);
  }
}

$ gcc -DNEWLIB -fno-builtin test_sqrt.c /localdisk/zimmerma/newlib-4.1.0/libm.a -lm
$ ./a.out
RNDN: 0x1.ff83fp+63
RNDZ: 0x1.ff83fp+63
RNDU: 0x1.ff83fp+63
RNDD: 0x1.ff83fp+63

The RNDZ and RNDD results are wrong. With glibc I get:

$ gcc -fno-builtin test_sqrt.c -lm
$ ./a.out
RNDN: 0x1.ff83fp+63
RNDZ: 0x1.ff83eep+63
RNDU: 0x1.ff83fp+63
RNDD: 0x1.ff83eep+63

Best regards,
Paul Zimmermann
lmpx.com only provides a reader for public news (NNTP) servers. It is not affiliated with the servers or forums shown here and is not responsible for the content of articles, which is written by their respective authors.