[Bug libquadmath/105101] incorrect rounding for sqrtq

[joseph at codesourcery dot com]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=105101

--- Comment #20 from joseph at codesourcery dot com <joseph at codesourcery dot com> ---
On Sat, 11 Jun 2022, already5chosen at yahoo dot com via Gcc-bugs wrote:

> On MSYS2 _Float128 and __float128 appears to be mostly the same thing, mapped
> to the same library routines with significant difference that _Float128 is not
> accessible from C++. Since all my test benches are written in C++ I can't even
> validate that what I wrote above is 100% true.
> 
> Also according to my understanding of glibc docs (not the clearest piece of
> text that I ever read) a relevant square root routine should be named
> sqrtf128().
> Unfortunately, nothing like that appears to be present in either math.h or in
> library. Am I doing something wrong?

The function should be sqrtf128 (present in glibc 2.26 and later on 
x86_64, x86, powerpc64le, ia64).  I don't know about support in non-glibc 
C libraries.

> Right now, there are only two [gcc] platforms with hw binary128 - IBM POWER and
> IBM z. I am not sure about the later, but the former has xssqrtqp instruction
> which is likely the right way to do sqrtq()/sqrtf128() on this platform. If z
> is the same, which sound likely, then implementation based on binary128
> mul/add/fma by now has no use cases at all.

That may well be the case for sqrt.

> > fma is a particularly tricky case because it *is* required to be correctly 
> > rounding, in all rounding modes, and correct rounding implies correct 
> > exceptions, *and* correct exceptions for fma includes getting right the 
> > architecture-specific choice of whether tininess is detected before or 
> > after rounding.
> 
> I suspect that by strict IEEE-754 rules sqrt() is the same as fma(), i.e. you
> have to calculate a root to infinite precision and then to round with
> accordance to current mode.

Yes, but fma has the extra complication of the architecture-specific 
tininess detection rules (to quote IEEE 754, "The implementer shall choose 
how tininess is detected [i.e. from the two options listed immediately 
above this quote in IEEE 754], but shall detect tininess in the same way 
for all operations in radix two"), which doesn't apply to sqrt because 
sqrt results can never underflow.

(I expect the existing soft-fp version of fma in glibc to be rather better 
optimized than the soft-fp version of sqrt.)

> I don't quite or understand why you say that. First, I don't remember using any
> double math in the variant of sqrtq() posted here. So, may be, you meant
> division?
> Here I use doable math, but IMO, I use it in a way that never causes
> exceptions.

Yes, the double division.  If that division can ever be inexact when the 
result of the square root itself is exact, you have a problem (which in 
turn could lead to actual incorrectly rounded results from other functions 
such as the square root operations rounding to a narrower type, when the 
round-to-odd versions of those functions are used, because the 
round-to-odd technique relies on correct "inexact" exceptions).