[Bug libquadmath/114623] New: sqrt

[Bug libquadmath/114623] New: sqrt

[g.peterhoff@t-online.de]

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

            Bug ID: 114623
           Summary: sqrt
           Product: gcc
           Version: unknown
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: libquadmath
          Assignee: unassigned at gcc dot gnu.org
          Reporter: g.peterhoff@t-online.de
  Target Milestone: ---

Hello,
sqrt does not work for std::numeric_limits<__float128>::max().
I have not checked other (special) values, perhaps the problem also occurs
there.

Please see https://godbolt.org/z/bx8or94v7

regards
Gero

[Bug libquadmath/114623] sqrt

[pinskia at gcc dot gnu.org]

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

--- Comment #1 from Andrew Pinski <pinskia at gcc dot gnu.org> ---
Please don't just link against godbolt and attach the preprocessed source.

[Bug libquadmath/114623] sqrtq and std::numeric_limits<__float128>::max()

[g.peterhoff@t-online.de]

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

--- Comment #2 from g.peterhoff@t-online.de ---
#include <boost/cstdfloat.hpp>
#include <string_view>
#include <charconv>
#include <iostream>
#include <limits>
#include <array>
#include <cmath>

void print_hex(const std::float128_t value)
{
        std::array<char, 1024>
                buffer{};
        const std::to_chars_result
                result{std::to_chars(buffer.data(),
buffer.data()+buffer.size(), value, std::chars_format::hex)};

        std::cout << std::string_view{buffer.data(), result.ptr} << std::endl;
}

template <typename Type>
void    print_sqrt_max_hex()
{
        using limits = std::numeric_limits<Type>;

        print_hex(std::sqrt(limits::max()));
}

int main()
{
        print_sqrt_max_hex<std::float128_t>();
        print_sqrt_max_hex<__float128>();

        return EXIT_SUCCESS;
}

gets
1.ffffffffffffffffffffffffffffp+8191
1p+8192

[Bug libquadmath/114623] sqrtq and std::numeric_limits<__float128>::max()

[pinskia at gcc dot gnu.org]

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

--- Comment #3 from Andrew Pinski <pinskia at gcc dot gnu.org> ---
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=83800#c4

Basically libquadmath should be using sqrt128f if it exist for sqrtq instead of
doing a version itself ... But libquadmath is normally only used for Fortran
really ...

[Bug libquadmath/114623] sqrtq and std::numeric_limits<__float128>::max()

[g.peterhoff@t-online.de]

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

--- Comment #4 from g.peterhoff@t-online.de ---
That is precisely the design error of C/C++/etc. There should be no
float/double/long double/__float128/etc, but *only* floatN_t. Then there
wouldn't be these discrepancies (if necessary you have to emulate by SW).
But that's just my humble opinion ... and now we have to face reality and make
the best of it. 
One step might be to put std::float128_t and __float128 on a common/uniform
code base :-)

cu
Gero

[Bug libquadmath/114623] sqrtq and std::numeric_limits<__float128>::max()

[jakub at gcc dot gnu.org]

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

Jakub Jelinek <jakub at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |jakub at gcc dot gnu.org

--- Comment #5 from Jakub Jelinek <jakub at gcc dot gnu.org> ---
(In reply to Andrew Pinski from comment #3)
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=83800#c4
> 
> Basically libquadmath should be using sqrt128f if it exist for sqrtq instead
> of doing a version itself ...

That doesn't make any sense.
sqrtf128 is only available in recent glibcs, not in any other C library.
libquadmath by design contains copies of the (old) glibc IEEE754 quad code, so
that
it provides the support independently from the C library.  Has been doing that
before
any *f128 support has been added to glibc, and needs to do it even now to
support other C libraries.
The sqrt* case was an exception, the glibc sysdeps/ieee754/ldbl-128/ didn't
contain any sqrt implementation, as arches which supported ldbl-128 as long
double at that point had hw sqrt support, so libquadmath has just an
approximation.
What could be done is use the soft-fp implementation like is done even in glibc
for x86.  We have the soft-fp code in gcc, but would need to arrange the right
sfp-machine.h etc.

[Bug libquadmath/114623] sqrtq and std::numeric_limits<__float128>::max()

[jakub at gcc dot gnu.org]

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

--- Comment #6 from Jakub Jelinek <jakub at gcc dot gnu.org> ---
(In reply to g.peterhoff from comment #4)
> That is precisely the design error of C/C++/etc. There should be no
> float/double/long double/__float128/etc, but *only* floatN_t.

If you don't want to use float/double/long double and just floatN_t, then just
use it.
And, it makes no sense to try to use __float128 at all when float128_t is
supported.

[Bug libquadmath/114623] sqrtq and std::numeric_limits<__float128>::max()

[jakub at gcc dot gnu.org]

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

--- Comment #7 from Jakub Jelinek <jakub at gcc dot gnu.org> ---
Created attachment 57900
  --> https://gcc.gnu.org/bugzilla/attachment.cgi?id=57900&action=edit
gcc14-pr114623.patch

Anyway, here is an untested patch to use soft-fp implementation for sqrtq for
the positive finite arguments to make them correctly rounded.

[Bug libquadmath/114623] sqrtq and std::numeric_limits<__float128>::max()

[cvs-commit at gcc dot gnu.org]

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

--- Comment #8 from GCC Commits <cvs-commit at gcc dot gnu.org> ---
The master branch has been updated by Jakub Jelinek <jakub@gcc.gnu.org>:

https://gcc.gnu.org/g:481ba4fb5fce8257f5dbeb994dac2748c0237fa2

commit r14-9853-g481ba4fb5fce8257f5dbeb994dac2748c0237fa2
Author: Jakub Jelinek <jakub@redhat.com>
Date:   Tue Apr 9 08:17:25 2024 +0200

    libquadmath: Use soft-fp for sqrtq finite positive arguments [PR114623]

    sqrt should be 0.5ulp precise, but the current implementation is less
    precise than that.
    The following patch uses the soft-fp code (like e.g. glibc for x86) for it
    if possible.  I didn't want to replicate the libgcc infrastructure for
    choosing the right sfp-machine.h, so the patch just uses a single generic
    implementation.  As the code is used solely for the finite positive
arguments,
    it shouldn't generate NaNs (so the exact form of canonical QNaN/SNaN is
    irrelevant), and sqrt for these shouldn't produce underflows/overflows
either,
    for < 1.0 arguments it always returns larger values than the argument and
for
    > 1.0 smaller values than the argument.

    2024-04-09  Jakub Jelinek  <jakub@redhat.com>

            PR libquadmath/114623
            * sfp-machine.h: New file.
            * math/sqrtq.c: Include from libgcc/soft-fp also soft-fp.h and
quad.h
            if possible.
            (USE_SOFT_FP): Define in that case.
            (sqrtq): Use soft-fp based implementation for the finite positive
            arguments if possible.