[Bug libstdc++/77776] C++17 std::hypot implementation is poor

[g.peterhoff--- via Gcc-bugs]

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

--- Comment #25 from g.peterh...@t-online.de ---
Hi Matthias,
to get good results on average (all FP-types: (B)FP16..FP128,
scalar/vectorized(SIMD)/parallel/...) this algorithm seems to me (so far) to be
suitable:

template <typename Type>
inline constexpr Type   hypot_gp(Type x, Type y, Type z)        noexcept
{
        using limits = std::numeric_limits<Type>;

        constexpr Type
                sqrt_min        = std::sqrt(limits::min()),
                sqrt_max        = std::sqrt(limits::max()),
                scale_up        = std::exp2( Type(limits::max_exponent/2)),
                scale_down      = std::exp2(-Type(limits::max_exponent/2)),
                zero            = 0;

        x = std::abs(x);
        y = std::abs(y);
        z = std::abs(z);

        if (!(std::isnormal(x) && std::isnormal(y) && std::isnormal(z)))
[[unlikely]]
        {
                if (std::isinf(x) | std::isinf(y) | std::isinf(z))      return
limits::infinity();
                else if (std::isnan(x) | std::isnan(y) | std::isnan(z)) return
limits::quiet_NaN();
                else
                {
                        const bool
                                xz{x == zero},
                                yz{y == zero},
                                zz{z == zero};

                        if (xz)
                        {
                                if (yz)         return zz ? zero : z;
                                else if (zz)    return y;
                        }
                        else if (yz && zz)      return x;
                }
        }

        if (x > z) std::swap(x, z);
        if (y > z) std::swap(y, z);

        if ((z >= sqrt_min) && (z <= sqrt_max)) [[likely]]
        {
                //      no scale
                return std::sqrt(__builtin_assoc_barrier(x*x + y*y) + z*z);
        }
        else
        {
                const Type
                        scale = (z >= sqrt_min) ? scale_down : scale_up;

                x *= scale;
                y *= scale;
                z *= scale;
                return std::sqrt(__builtin_assoc_barrier(x*x + y*y) + z*z);
        }
}