incorrectly rounded square root

[Jeff Johnston]

Hi Paul,

Not all platforms supply a machine implementation of fesetround.  Which
platform are you using?

-- Jeff J.

On Tue, May 4, 2021 at 4:09 AM Paul Zimmermann <Paul.Zimmermann@inria.fr>
wrote:

>        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
>
>