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