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

["g.peterhoff@t-online.de" <gcc-bugzilla@gcc.gnu.org>]

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

--- Comment #19 from g.peterhoff@t-online.de ---
> So, no need to use frexp/ldexp, just comparisons of hi above against sqrt of
> (max finite / 3), in that case scale by multiplying all 3 args by some
> appropriate scale constant, and similarly otherwise if lo1 is too small by
> some large scale.

I don't really know. With frexp/ldexp you probably get the highest accuracy
(even if it is probably slower) instead of doing it manually. The problem is to
determine suitable scaling factors and to adjust the (return)values
accordingly. I have implemented both cases.

Error
* In the case (x==y && y==z), x*std::sqrt(T(3)) must not simply be returned, as
this can lead to an overflow (inf).

Generally
* Instead of using fmin/fmax to determine the values hi,lo1,lo0, it is better
to sort x,y,z. This is faster and clearer and no additional variables need to
be introduced.
* It also makes sense to consider the case (x==0 && y==0 && z==0).

Optimizations
* You were probably wondering why I wrote "if (std::isinf(x) | std::isinf(y) |
std::isinf(z))", for example. This is intentional. The problem is that gcc
almost always produces branch code for logical operations, so *a lot* of
conditional jumps. By using arithmetic operations, so instead of || && just |
&, I can get it to generate only actually necessary conditional jumps or
cmoves. branchfree code is always better.


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

        constexpr T
                zero = 0;

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

        if (std::isinf(x) | std::isinf(y) | std::isinf(z))  [[unlikely]]
                return limits::infinity();
        if (std::isnan(x) | std::isnan(y) | std::isnan(z))      [[unlikely]]
                return limits::quiet_NaN();
        if ((x==zero) & (y==zero) & (z==zero))  [[unlikely]]
                return zero;
        if ((y==zero) & (z==zero))      [[unlikely]]
                return x;
        if ((x==zero) & (z==zero))      [[unlikely]]
                return y;
        if ((x==zero) & (y==zero))      [[unlikely]]
                return z;

        auto sort = [](T& a, T& b, T& c)        constexpr noexcept -> void
        {
                if (a > b) std::swap(a, b);
                if (b > c) std::swap(b, c);
                if (a > b) std::swap(a, b);
        };

        sort(x, y, z);  //      x <= y <= z

        int
                exp = 0;

        z = std::frexp(z, &exp);
        y = std::ldexp(y, -exp);
        x = std::ldexp(x, -exp);

        T
                sum = x*x + y*y;

        sum += z*z;
        return std::ldexp(std::sqrt(sum), exp);
}

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

        auto prev_power2 = [](const T value)    constexpr noexcept -> T
        {
                return std::exp2(std::floor(std::log2(value)));
        };

        constexpr T
                sqrtmax         = std::sqrt(limits::max()),
                scale_up        = prev_power2(sqrtmax),
                scale_down      = T(1) / scale_up,
                zero            = 0;

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

        if (std::isinf(x) | std::isinf(y) | std::isinf(z))  [[unlikely]]
                return limits::infinity();
        if (std::isnan(x) | std::isnan(y) | std::isnan(z))      [[unlikely]]
                return limits::quiet_NaN();
        if ((x==zero) & (y==zero) & (z==zero))  [[unlikely]]
                return zero;
        if ((y==zero) & (z==zero))      [[unlikely]]
                return x;
        if ((x==zero) & (z==zero))      [[unlikely]]
                return y;
        if ((x==zero) & (y==zero))      [[unlikely]]
                return z;

        auto sort = [](T& a, T& b, T& c)        constexpr noexcept -> void
        {
                if (a > b) std::swap(a, b);
                if (b > c) std::swap(b, c);
                if (a > b) std::swap(a, b);
        };

        sort(x, y, z);  //      x <= y <= z

        const T
                scale = (z > sqrtmax) ? scale_down : (z < 1) ? scale_up : 1;

        x *= scale;
        y *= scale;
        z *= scale;

        T
                sum = x*x + y*y;

        sum += z*z;
        return std::sqrt(sum) / scale;
}


regards
Gero