#ifdef THIS_TYPE
static THIS_TYPE
THIS_FUNC (ulp) (THIS_TYPE val)
{
  if (__builtin_isnormal (val))
    return THIS_FUNC (ldexp) (THIS_LIT (1.0), THIS_FUNC (ilogb) (val) - THIS_MANT_DIG + 1);
  else
    return THIS_FUNC (ldexp) (THIS_LIT (1.0), THIS_MIN_EXP - THIS_MANT_DIG);
}

static void
THIS_FUNC (test) (THIS_TYPE (*fn) (THIS_TYPE),
		  int (*mpfr_fn) (mpfr_ptr, mpfr_srcptr, mpfr_rnd_t),
		  const char *name,
		  unsigned long count)
{
  char buf1[256], buf2[256], buf3[256];
  mpfr_set_default_prec (THIS_MANT_DIG);
  mpfr_t v;
  mpfr_init2 (v, THIS_MANT_DIG);
  THIS_TYPE max_ulp = THIS_LIT (0.0);
  volatile THIS_TYPE val = THIS_LIT (0.0);
  int m = 0;
  for (unsigned long i = 0; i < count; ++i)
    {
      if (m == 0)
	m = 1;
      else if (m <= 10)
	{
	  val = THIS_FUNC (nextafter) (val, THIS_LIT (1.0));
	  if ((m == 1 && val == THIS_MIN) || m > 1)
	    ++m;
	}
      else if (m == 11)
	{
	  val = THIS_LIT (1.0);
	  for (int j = 0; j < 10; j++)
	    val = THIS_FUNC (nextafter) (val, THIS_LIT (0.0));
	  ++m;
	}
      else if (m <= 32)
	{
	  val = THIS_FUNC (nextafter) (val, THIS_LIT (1.0));
	  ++m;
	}
      else if (m == 33)
	{
	  val = THIS_MAX;
	  for (int j = 0; j < 10; j++)
	    val = THIS_FUNC (nextafter) (val, THIS_LIT (0.0));
	  ++m;
	}
      else if (m <= 45)
	{
	  val = THIS_FUNC (nextafter) (val, THIS_LIT (1.0));
	  ++m;
	}
      else
	val = THIS_FUNC (get_rand) ();
      if (__builtin_isnan (val))
	continue;
      THIS_TYPE given = fn (val);
      THIS_MPFR_SET (v, val, MPFR_RNDN);
      mpfr_fn (v, v, MPFR_RNDN);
      THIS_TYPE expected = THIS_MPFR_GET (v, MPFR_RNDN);
      if ((!__builtin_isnan (given)) != (!__builtin_isnan (expected))
	  || __builtin_isinf (given) != __builtin_isinf (expected))
	{
	  THIS_SNPRINTF (buf1, val);
	  THIS_SNPRINTF (buf2, given);
	  THIS_SNPRINTF (buf3, expected);
	  printf ("%s (%s) = %s rather than %s\n", name, buf1, buf2, buf3);
	}
      else if (!__builtin_isnan (given) && !__builtin_isinf (given))
	{
	  THIS_TYPE this_ulp
	    = THIS_FUNC (fabs) (given - expected) / THIS_FUNC (ulp) (expected);
	  if (this_ulp > max_ulp)
	    max_ulp = this_ulp;
	}
    }
  printf ("%s max error %.1fulp\n", name, (float) max_ulp);
}

#undef THIS_TYPE
#undef THIS_LIT
#undef THIS_FUNC
#undef THIS_MIN_EXP
#undef THIS_MANT_DIG
#undef THIS_MIN
#undef THIS_MAX
#undef THIS_MPFR_SET
#undef THIS_MPFR_GET
#undef THIS_SNPRINTF
#else
#define _GNU_SOURCE 1
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <fenv.h>

#if defined(__FLT128_DIG__) && defined(__GLIBC_PREREQ)
#if __GLIBC_PREREQ (2, 26)
#define TEST_FLT128
#define MPFR_WANT_FLOAT128
#endif
#endif

#include <gmp.h>
#include <mpfr.h>

static long rand_n;
static int rand_c;

static uint32_t
get_rand32 (void)
{
  uint32_t ret = 0;
  if (rand_c == 0)
    {
      ret = random () & 0x7fffffff;
      rand_c = 31;
    }
  else
    ret = rand_n & (((uint32_t) 1 << rand_c) - 1);
  ret <<= (32 - rand_c);
  rand_n = random ();
  ret |= rand_n & (((uint32_t) 1 << (32 - rand_c)) - 1);
  rand_n >>= (32 - rand_c);
  rand_c = 31 - (32 - rand_c);
  return ret;
}

static uint64_t
get_rand64 (void)
{
  return (((uint64_t) get_rand32 ()) << 32) | get_rand32 ();
}

static float
get_randf (void)
{
  uint32_t i = get_rand32 ();
  float f;
  memcpy (&f, &i, sizeof (f));
  return f;
}

static double
get_rand (void)
{
  uint64_t i = get_rand64 ();
  double d;
  memcpy (&d, &i, sizeof (d));
  return d;
}

static long double
get_randl (void)
{
  long double ld;
  uint64_t i = get_rand64 ();
  memcpy (&ld, &i, sizeof (i));
  if (sizeof (long double) == 12)
    {
      uint32_t j = get_rand32 ();
      memcpy ((char *) &ld + 8, &j, sizeof (j));
    }
  else if (sizeof (long double) == 16)
    {
      i = get_rand64 ();
      memcpy ((char *) &ld + 8, &i, sizeof (i));
    }
  return ld;
}

#ifdef TEST_FLT128
static long double
get_randf128 (void)
{
  _Float128 f128;
  uint64_t i = get_rand64 ();
  memcpy (&f128, &i, sizeof (i));
  i = get_rand64 ();
  memcpy ((char *) &f128 + 8, &i, sizeof (i));
  return f128;
}
#endif

#define THIS_TYPE float
#define THIS_LIT(v) v##f
#define THIS_FUNC(v) v##f
#define THIS_MIN_EXP __FLT_MIN_EXP__
#define THIS_MANT_DIG __FLT_MANT_DIG__
#define THIS_MIN __FLT_MIN__
#define THIS_MAX __FLT_MAX__
#define THIS_MPFR_SET mpfr_set_flt
#define THIS_MPFR_GET mpfr_get_flt
#define THIS_SNPRINTF(buf, x) snprintf ((buf), sizeof (buf), "%a", (x));
#include "ulp-tester.c"

#define THIS_TYPE double
#define THIS_LIT(v) v
#define THIS_FUNC(v) v
#define THIS_MIN_EXP __DBL_MIN_EXP__
#define THIS_MANT_DIG __DBL_MANT_DIG__
#define THIS_MIN __DBL_MIN__
#define THIS_MAX __DBL_MAX__
#define THIS_MPFR_SET mpfr_set_d
#define THIS_MPFR_GET mpfr_get_d
#define THIS_SNPRINTF(buf, x) snprintf ((buf), sizeof (buf), "%a", (x));
#include "ulp-tester.c"

#define THIS_TYPE long double
#define THIS_LIT(v) v##L
#define THIS_FUNC(v) v##l
#define THIS_MIN_EXP __LDBL_MIN_EXP__
#define THIS_MANT_DIG __LDBL_MANT_DIG__
#define THIS_MIN __LDBL_MIN__
#define THIS_MAX __LDBL_MAX__
#define THIS_MPFR_SET mpfr_set_ld
#define THIS_MPFR_GET mpfr_get_ld
#define THIS_SNPRINTF(buf, x) snprintf ((buf), sizeof (buf), "%La", (x));
#include "ulp-tester.c"

#ifdef TEST_FLT128
#define THIS_TYPE _Float128
#define THIS_LIT(v) v##F128
#define THIS_FUNC(v) v##f128
#define THIS_MIN_EXP __FLT128_MIN_EXP__
#define THIS_MANT_DIG __FLT128_MANT_DIG__
#define THIS_MIN __FLT128_MIN__
#define THIS_MAX __FLT128_MAX__
#define THIS_MPFR_SET mpfr_set_float128
#define THIS_MPFR_GET mpfr_get_float128
#define THIS_SNPRINTF(buf, x) strfromf128 ((buf), sizeof (buf), "%a", (x));
#include "ulp-tester.c"
#else
#define testf128(fn, mpfr_fn, name, count) do { } while (0)
#endif

int
main (int argc, const char **argv)
{
  const char *arg;
  char *endptr;
  (void) argc;
  if (argc <= 1)
    arg = "";
  else
    arg = argv[1];
  unsigned long count = strtoul (arg, &endptr, 10);
  if (endptr == arg)
    {
      fprintf (stderr, "ulp-tester number_of_iterations rnd\n");
      return 1;
    }
  const char *rnd = "tonearest";
  if (argc >= 3)
    rnd = argv[2];
  if (strcmp (rnd, "upward") == 0)
    fesetround (FE_UPWARD);
  else if (strcmp (rnd, "downward") == 0)
    fesetround (FE_DOWNWARD);
  else if (strcmp (rnd, "towardzero") == 0)
    fesetround (FE_TOWARDZERO);
#define TESTS(fn) \
  testf (fn##f, mpfr_##fn, #fn "f", count);		\
  test (fn, mpfr_##fn, #fn, count);			\
  testl (fn##l, mpfr_##fn, #fn "l", count);		\
  testf128 (fn##f128, mpfr_##fn, #fn "f128", count)
  TESTS (sqrt);
  TESTS (sin);
  TESTS (cos);
  TESTS (exp10);
}
#endif