incorrectly rounded square root

[Paul Zimmermann <Paul.Zimmermann@inria.fr>]

       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