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

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

int
main ()
{
  float x = 0x1.0c2d26p+9;
  float 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_TOWARDZERO);
  y = ynf (42, x);
  mpfr_init2 (t, 24);
  mpfr_set_flt (t, x, MPFR_RNDN);
  mpfr_yn (t, 42, t, MPFR_RNDN);
  z = mpfr_get_flt (t, MPFR_RNDU);
  printf ("x=%a\n", x);
  printf ("libc: yn(42,x)=%a\n", y);
  printf ("mpfr: yn(42,x)=%a\n", z);
  printf ("ulp difference: %f\n", ulp (y, z));
  return 0;
}