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

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

int
main ()
{
  long double x = 0xe.4c120b3ec76c99dp-4L;
  long 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 ());
  fesetround (FE_UPWARD);
  y = y0l (x);
  mpfr_init2 (t, 64);
  mpfr_set_ld (t, x, MPFR_RNDN);
  mpfr_y0 (t, t, MPFR_RNDN);
  z = mpfr_get_ld (t, MPFR_RNDN);
  printf ("x=%La\n", x);
  printf ("libc: y0(x)=%La\n", y);
  printf ("mpfr: y0(x)=%La\n", z);
  printf ("ulp difference: %f\n", ulp (y, z));
  return 0;
}