ldexp with directed rounding

ldexp with directed rounding

[Paul Zimmermann]

       Hi,

it seems that Newlib does not honor the rounding mode for ldexp(0x1p-1074,-1):
for rounding updwards it should return 0x1p-1074 and not 0.

Test program:

#include <stdio.h>
#include <math.h>
#include <fenv.h>

int main()
{
  volatile double x = 1.0, y;
  while (x * 0.5 != 0)
    x = x * 0.5;
  fesetround (FE_TONEAREST);
  y = ldexp (x, -1);
  printf ("FE_TONEAREST:  y=%la\n", y);
  fesetround (FE_TOWARDZERO);
  y = ldexp (x, -1);
  printf ("FE_TOWARDZERO: y=%la\n", y);
  fesetround (FE_UPWARD);
  y = ldexp (x, -1);
  printf ("FE_UPWARD:     y=%la\n", y);
  fesetround (FE_DOWNWARD);
  y = ldexp (x, -1);
  printf ("FE_DOWNWARD:   y=%la\n", y);
}

With newlib-4.3.0.20230120 on x86_64:

zimmerma@salade:~/svn/tbd/20/src/binary64$ gcc -fno-builtin check_ldexp.c /localdisk/zimmerma/${NEWLIB}/libm.a
zimmerma@salade:~/svn/tbd/20/src/binary64$ ./a.out
FE_TONEAREST:  y=0x0p+0
FE_TOWARDZERO: y=0x0p+0
FE_UPWARD:     y=0x0p+0
FE_DOWNWARD:   y=0x0p+0

With GNU libc:

zimmerma@salade:~/svn/tbd/20/src/binary64$ gcc -fno-builtin check_ldexp.c -lm
zimmerma@salade:~/svn/tbd/20/src/binary64$ ./a.out
FE_TONEAREST:  y=0x0p+0
FE_TOWARDZERO: y=0x0p+0
FE_UPWARD:     y=0x0.0000000000001p-1022
FE_DOWNWARD:   y=0x0p+0

This test was suggested by Fred Tydeman.

Paul