https://gcc.gnu.org/bugzilla/show_bug.cgi?id=105101
Jakub Jelinek <jakub at gcc dot gnu.org> changed:
What |Removed |Added
----------------------------------------------------------------------------
CC| |jakub at gcc dot gnu.org,
| |jsm28 at gcc dot gnu.org
--- Comment #1 from Jakub Jelinek <jakub at gcc dot gnu.org> ---
Trying
#include <math.h>
#include <stdio.h>
int
main ()
{
volatile long double ld = 0x2.4e5658b5f6d2a1ec21e3829dd7f8p+0L;
ld = sqrtl (ld);
printf ("%.50La\n", ld);
return 0;
}
on s390x-linux shows:
0x1.84bfedcba255aeaf7e99184e6f3b0000000000000000000000p+0
but that is implemented using a hw instruction.
#include <quadmath.h>
#include <stdio.h>
int
main ()
{
volatile __float128 ld = 0x2.4e5658b5f6d2a1ec21e3829dd7f8p+0Q;
ld = sqrtq (ld);
char buf[70];
quadmath_snprintf (buf, 69, "%.50Qa", ld);
printf ("%s\n", buf);
return 0;
}
on x86_64 indeed prints ...6f3a
The first testcase also prints ...6f3b on powerpc64le-linux, but again
that doesn't say much, because powerpc64le-linux __ieee754_sqrtf128
comes from sysdeps/powerpc/powerpc64/le/fpu/e_sqrtf128.c where in my case it is
implemented using soft-fp (and when glibc is configured for power9 or later
using hw instruction).
Seems libquadmath sqrtq.c isn't based on any glibc code, instead it uses sqrt
or sqrtl for initial approximation and then
/* Two Newton iterations. */
y -= 0.5q * (y - x / y);
y -= 0.5q * (y - x / y);
(or for sqrtl one iteration).
I bet that doesn't and can't guarantee correct rounding.
So, do we need to switch to soft-fp based implementation for it (we have a copy
already in libgcc/soft-fp), or do we have other options?
