#define _GNU_SOURCE /* to define expf128 */
#include <stdio.h>
#include <stdint.h>
#include <math.h>
#include <assert.h>
#include <gmp.h>
#define MPFR_WANT_FLOAT128
#include <mpfr.h>

/* if TIMINGS is not defined, check accuracy */
// #define TIMINGS

/* number of tests */
#define N 1000000UL

/* Tests are performed with number in [XMIN,XMAX] */
#define XMIN -9.99514646814214379832262663921344q
#define XMAX -XMIN

/* function to check */

#ifdef USE_GLIBC
#define FUNCTION expf128 /* for the glibc version */
#include <gnu/libc-version.h>
#else
#define FUNCTION my_expf128
#endif

extern _Float128 my_expf128 (const _Float128 x);

static void
print_float128 (_Float128 x)
{
  mpfr_t t;
  mpfr_init2 (t, 113);
  mpfr_set_float128 (t, x, MPFR_RNDN);
  mpfr_printf ("%Re", t);
  mpfr_clear (t);
}

#ifdef TIMINGS
static void
check_speed (void)
{
  _Float128 x = XMIN;
  _Float128 h = (XMAX - XMIN) / N;
  _Float128 s = 0;
  for (unsigned long i = 0; i < N; i++)
    {
      s += FUNCTION (x);
      x += h;
    }
  printf ("s="); print_float128 (s); printf ("\n");
}
#endif

/* return exp(x) using MPFR, assuming no denormal */
_Float128
correct_exp (_Float128 x, mpfr_t y)
{
  mpfr_set_float128 (y, x, MPFR_RNDN);
  mpfr_exp (y, y, MPFR_RNDN);
  return mpfr_get_float128 (y, MPFR_RNDN);
}

#ifndef TIMINGS
static unsigned long
ulp_error (_Float128 y, _Float128 z)
{
  _Float128 err, ulp;
  assert (isinfl ((long double) z) == 0);
  err = y - z;
  ulp = nextafterf128 (z, y) - z;
  err = fabsf128 (err / ulp);
  return (err >= (_Float128) ULONG_MAX) ? ULONG_MAX : (unsigned long) err;
}

static void
check_accuracy (void)
{
  _Float128 x = XMIN, y, z, maxx, maxy, maxz;
  _Float128 h = (XMAX - XMIN) / N;
  mpfr_t t;
  unsigned long correct = 0;
  unsigned long maxerr = 0, err;

  mpfr_init2 (t, 113);
  for (unsigned long i = 0; i < N; i++)
    {
      y = FUNCTION (x);
      z = correct_exp (x, t);
      if (y == z)
        correct ++;
      else
        {
          err = ulp_error (y, z);
          if (err > maxerr)
            {
              maxerr = err;
              maxx = x;
              maxy = y;
              maxz = z;
            }
        }
      x += h;
    }
  mpfr_clear (t);
  printf ("correct roundings: %lu/%lu max err=%lu ulp(s)\n",
          correct, N, maxerr);
  if (maxerr > 0)
    {
      printf ("maximal error for\nx=");
      print_float128 (maxx);
      printf ("\ny=");
      print_float128 (maxy);
      printf ("\nz=");
      print_float128 (maxz);
      printf ("\n");
    }
}
#endif

int
main()
{
  /* check that long long unsigned int has 64 bits (for _mulx_u64) */
  assert (sizeof (long long unsigned int) == sizeof (uint64_t));
#ifdef USE_GLIBC
  printf("GNU libc version: %s\n", gnu_get_libc_version ());
  printf("GNU libc release: %s\n", gnu_get_libc_release ());
#endif
#ifndef TIMINGS
  check_accuracy ();
#else
  check_speed ();
#endif
}