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
/* 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
