with the test program below I get on a x86_64 processor:

GNU libc version: 2.18.90
GNU libc release: development
x=8.9357301493606989e-01
libc: y0(x)=-3.4749010007703585e-06
mpfr: y0(x)=-3.4749010007607514e-06
ulp difference: -22684.000000

I will keep track of such errors on http://www.loria.fr/~zimmerma/glibc/.

Paul Zimmermann

#include <stdio.h>
#include <math.h>
#include <mpfr.h>
#include <gnu/libc-version.h>

double
ulp (double x, double y)
{
  double z = x - y;
  int ey;
  y = frexp (y, &ey);
  y = ldexp (1.0, ey - 53);
  return z / y;
}

int
main ()
{
  double x = 8.9357301493606989e-1;
  double y, z;
  mpfr_t t;

  printf("GNU libc version: %s\n", gnu_get_libc_version ());
  printf("GNU libc release: %s\n", gnu_get_libc_release ());
  y = y0 (x);
  mpfr_init2 (t, 53);
  mpfr_set_d (t, x, MPFR_RNDN);
  mpfr_y0 (t, t, MPFR_RNDN);
  z = mpfr_get_d (t, MPFR_RNDN);
  printf ("x=%.16e\n", x);
  printf ("libc: y0(x)=%.16e\n", y); /* -3.4749010007703585e-06 */
  printf ("mpfr: y0(x)=%.16e\n", z); /* -3.4749010007607514e-06 */
  printf ("ulp difference: %f\n", ulp (y, z)); /* -22684.000000 */
  return 0;
}