#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, t, twoy;

  twoy = y + y;
  if (twoy != twoy + twoy) /* twoy is finite */
    t = nextafterf (y, twoy);
  else /* twoy is +Inf or -Inf */
    t = nextafter (y, (float) 0);
  t = y - t;
  if (t < (float) 0)
    t = -t;
  /* now t = ulp(y) */
  return z / t;
}

int
main ()
{
  float x = 0x3.8b7f28p-128;
  float y, z;
  mpfr_t t;
  int inex, runderflow;

  printf("GNU libc version: %s\n", gnu_get_libc_version ());
  printf("GNU libc release: %s\n", gnu_get_libc_release ());
  fesetround (FE_TOWARDZERO);
  feclearexcept (FE_ALL_EXCEPT);
  y = erff (x);
  runderflow = fetestexcept (FE_UNDERFLOW);
  mpfr_set_emin (-125 - 24 + 1);
  mpfr_set_emax (128);
  mpfr_init2 (t, 24);
  mpfr_set_flt (t, x, MPFR_RNDN);
  inex = mpfr_erf (t, t, MPFR_RNDZ);
  inex = mpfr_subnormalize (t, inex, MPFR_RNDZ);
  z = mpfr_get_flt (t, MPFR_RNDD);
  printf ("x=%a (%.9e)\n", x, x);
  printf ("libc: erff(x)=%a (%.9e) underflow=%d\n", y, y, runderflow);
  printf ("mpfr: erff(x)=%a (%.9e)\n", z, z);
  printf ("ulp difference: %f\n", ulp (y, z));
  return 0;
}