Re: [PATCH] Add targetm.libm_function_max_error

[Jakub Jelinek <jakub@redhat.com>]

[-- Attachment #1: Type: text/plain, Size: 3745 bytes --]

On Wed, Apr 26, 2023 at 04:10:32PM +0000, Michael Matz wrote:
> On Wed, 26 Apr 2023, Jakub Jelinek via Gcc-patches wrote:
> 
> > For glibc I've gathered data from:
> > 4) using attached ulp-tester.c (how to invoke in file comment; tested
> >    both x86_64, ppc64, ppc64le 50M pseudo-random values in all 4 rounding
> >    modes, plus on x86_64 float/double sin/cos using libmvec - see
> >    attached libmvec-wrapper.c as well)
> 
> That ulp-tester.c file as attached here is not testing what you think it 
> does.  (1) It doesn't compile as it doesn't #define the TESTS macro in the 

Oops, (1) was from a last minute change.  Initially (what was posted already
back in November) it had no libmvec support, then I wanted to try to
add libmvec-wrapper.so as a real wrapper which would implement
sinf/sin/cosf/cos using the libmvec APIs.  Unfortunately it turns out the
libmvec APIs in certain corner cases call the scalar functions and then it
obviously crashed due to endless recursion; I guess I could have some
recursion counter and use dlsym RTLD_NEXT, but just chose to use different
names and quickly adjusted the tester and only tested it with the libmvec
stuff.

> !LIBMVEC_TEST case, and (2) it almost never progresses 'm', the status 
> variable used before the random numbers start, to beyond 1: you start with 
> nextafter(0.0, 1.0), which is the smallest subnormal number (with a ERANGE 
> error, but that's ignored), and you test for equality with THIS_MIN, the 
> smallest normal (!) number, until you start incrementing 'm'.

Oops, you're right.  Obviously trying to go through all denormals
exhaustively is a bad idea, attaching adjusted ulp-tester.c, which does that
only for 8192 smallest denormals, then just power of two denormals until
minimum normalized etc.  Also, I've adjusted it such that for the non-random
values it tests both positive and negative values.
And, for IBM double double it now attempts to construct a valid IBM double
double random value rather than sometimes invalid.

With this, results are on x86_64 are still
		sqrtf	sqrt	sqrtl	sqrtf128
tonearest	0	0	0	0
upward		0	0	0	0
downward	0	0	0	0
towardzero	0	0	0	0
		sinf	sin	sinl	sinf128
tonearest	1	1	1	1
upward		1	1	3	3
downward	1	1	3	3
towardzero	1	1	2	2
libmvec near	2	3-4
		cosf	cos	cosl	cosf128
tonearest	1	1	1	1
upward		1	1	3	3
downward	1	1	3	3
towardzero	1	1	2	2
libmvec near	2	3-4
On ppc64
		sqrtf	sqrt	sqrtl
tonearest	0	0	2.5
upward		0	0	2
downward	0	0	2
towardzero	0	0	3
		sinf	sin	sinl
tonearest	1	0	inf
upward		3	1	170141112472035618950840819900504604672
downward	2	1	5285505939519009108193147419208187904
towardzero	1	1	170141102330830817125005607926878961664
		cosf	cos	cosl
tonearest	1	0	83814836763238928169050593715445301248
upward		2	1	83563811520779333270048610559903924224
downward	3	1	34088889209772606301573232603422523392
towardzero	2	1	34088889209772606301573232603422523392
On ppc64le
		sqrtf	sqrt	sqrtl
tonearest	0	0	2.5
upward		0	0	2
downward	0	0	2
towardzero	0	0	3
		sinf	sin	sinl
tonearest	1	0	inf
upward		1	1	170141112472035618950840819900504604672
downward	1	1	5285505939519009108193147419208187904
towardzero	1	1	170141102330830817125005607926878961664
		cosf	cos	cosl
tonearest	1	0	83814836763238928169050593715445301248
upward		2	1	83563811520779333270048610559903924224
downward	3	1	34088889209772606301573232603422523392
towardzero	2	1	34088889209772606301573232603422523392

I guess I'll need to look at the IBM double double sinl/cosl results,
either it is some bug in my tester or the libm functions are useless.
But appart from the MODE_COMPOSITE_P cases, I think all the numbers are
within what the patch returns.
Even the sqrtl tonearest IBM double double case is larger than the libm ulps
(2.5 vs. 1).

	Jakub

[-- Attachment #2: ulp-tester.c --]
[-- Type: text/plain, Size: 8202 bytes --]

/* gcc -g -O2 -o ulp-tester ulp-tester.c -lmpfr -lm
   for i in tonearest upward downward towardzero; do ( ./ulp-tester 50000000 $i > /tmp/$i 2>&1 & ); done
   gcc -g -O2 -DLIBMVEC_TEST -o ulp-tester2 ulp-tester.c ./libmvec-wrapper.so -lmpfr -lm
   for j in b c d e; do for i in tonearest upward downward towardzero; do ( LIBMVEC_LEVEL=$j ./ulp-tester2 50000000 $i > /tmp/${i}-${j} 2>&1 & ); done; done
 */

#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 == 1)
	{
	  val = THIS_FUNC (nextafter) (val, THIS_LIT (1.0));
	  if (val == THIS_MIN)
	    m = 3;
	  else if (i == 16384)
	    ++m;
	}
      else if (m == 2)
	{
	  val = THIS_LIT (2.0) * val;
	  if (val >= THIS_MIN)
	    ++m;
	}
      else if (m <= 10)
	{
	  val = THIS_FUNC (nextafter) (val, THIS_LIT (1.0));
	  ++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;
      volatile THIS_TYPE val2 = val;
      do
	{
	  THIS_TYPE given = fn (val2);
	  THIS_MPFR_SET (v, val2, mrnd);
	  mpfr_fn (v, v, mrnd);
	  THIS_TYPE expected = THIS_MPFR_GET (v, mrnd);
	  if ((!__builtin_isnan (given)) != (!__builtin_isnan (expected))
	      || __builtin_isinf (given) != __builtin_isinf (expected))
	    {
	      THIS_SNPRINTF (buf1, val2);
	      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;
	    }
	  if (m <= 45 && !__builtin_signbit (val2))
	    {
	      val2 = -val;
	      ++i;
	    }
          else
	    break;
	}
      while (1);
    }
  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;
mpfr_rnd_t mrnd;

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)
    {
#if __LDBL_MANT_DIG__ == 106
      double d1, d2;
      memcpy (&d1, &i, sizeof (d1));
      uint64_t j = get_rand64 ();
      if (((j >> 52) & 0x7ff) > ((i >> 52) & 0x7ff) - 52)
	j = (j & 0x800fffffffffffffULL) | (((((i >> 52) & 0x7ff) - 54) + ((j >> 52) & 3)) << 52);
      memcpy (&d2, &j, sizeof (d2));
      ld = d1;
      ld += d2;
#else
      i = get_rand64 ();
      memcpy ((char *) &ld + 8, &i, sizeof (i));
#endif
    }
  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];
  mrnd = MPFR_RNDN;
  if (strcmp (rnd, "upward") == 0)
    {
      fesetround (FE_UPWARD);
      mrnd = MPFR_RNDU;
    }
  else if (strcmp (rnd, "downward") == 0)
    {
      fesetround (FE_DOWNWARD);
      mrnd = MPFR_RNDD;
    }
  else if (strcmp (rnd, "towardzero") == 0)
    {
      fesetround (FE_TOWARDZERO);
      mrnd = MPFR_RNDZ;
    }
#ifdef LIBMVEC_TEST
  extern float wsinf (float);
  extern float wcosf (float);
  extern double wsin (double);
  extern double wcos (double);
#define TESTS(fn) \
  testf (w##fn##f, mpfr_##fn, #fn "f", count);		\
  test (w##fn, mpfr_##fn, #fn, count);
  TESTS (sin);
  TESTS (cos);
#else
#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
}
#endif