https://gcc.gnu.org/bugzilla/show_bug.cgi?id=105101
--- Comment #8 from Steve Kargl <sgk at troutmask dot apl.washington.edu> --- On Sat, Apr 09, 2022 at 10:23:39AM +0000, jakub at gcc dot gnu.org wrote: > https://gcc.gnu.org/bugzilla/show_bug.cgi?id=105101 > > --- Comment #6 from Jakub Jelinek <jakub at gcc dot gnu.org> --- > (In reply to Thomas Koenig from comment #5) > > There is another, much worse, problem, reported and analyzed by "Michael S" > > on comp.arch. The code has > > > > #ifdef HAVE_SQRTL > > { > > long double xl = (long double) x; > > if (xl <= LDBL_MAX && xl >= LDBL_MIN) > > { > > /* Use long double result as starting point. */ > > y = (__float128) sqrtl (xl); > > > > /* One Newton iteration. */ > > y -= 0.5q * (y - x / y); > > return y; > > } > > } > > #endif > > > > which assumes that long double has a higher precision than > > normal double. On x86_64, this depends o the settings of the > > FPU flags, so a number like 0x1.06bc82f7b9d71dfcbddf2358a0eap-1024 > > is corrected with 32 ULP of error because there is only a single > > round of Newton iterations if the FPU flags are set to normal precision. > > That is only a problem on OSes that do that, I think mainly BSDs, no? > On Linux it should be fine (well, still not 0.5ulp precise, but not as bad as > when sqrtl is just double precision precise). > i686-*-freebsd sets the FPU to have 53 bits of precision for long double. It has the usual exponent range of an Intel 80-bit extended double.
