On Thu, Apr 27, 2023 at 10:59:47AM +0200, Jakub Jelinek via Gcc-patches wrote:
> 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).
The first really large error I see is for sinl with
x/2gx &val
0x748160ed90d9425b 0xefd8b811d6293294
i.e.
1.5926552660973502228303666578452949e+253
with most significant double being
1.5926552660973502e+253
and low double
-5.9963639272208416e+230
Now, 0x748 - 0x6fd is 75, which is much larger than
53, so the number has precision larger than 106 bits.
given is
-0.4025472157704263326278375983156912
and expected (mpfr computed)
-0.46994008859023245970759964236618727
But if I try on x86_64:
#define _GNU_SOURCE
#include <math.h>
int
main ()
{
_Float128 f, f2, f3, f4;
double d, d2;
f = 1.5926552660973502228303666578452949e+253f128;
d = 1.5926552660973502e+253;
f2 = d;
f2 += -5.9963639272208416e+230;
f3 = sinf128 (f);
f4 = sinf128 (f2);
d2 = sin (d);
return 0;
}
where I think f2 is what matches most closely the 106 bit precision value,
(gdb) p f
$7 = 1.5926552660973502228303666578452949e+253
(gdb) p f2
$8 = 1.59265526609735022283036665784527174e+253
(gdb) p f3
$9 = -0.277062522218693980443596385112227247
(gdb) p f4
$10 = -0.402547215770426332627837598315693221
and f4 is much closer to the given than to expected.
On the other side, GCC will really work only with the assumption
the numbers have 106-bit precision, so shouldn't care much about
exact precision in between the range boundaries.
I think I'll for now just trust for IBM double double
the ulps files rather than mpfr.
Jakub
