[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 #29 from g.peterh...@t-online.de ---
(In reply to Jakub Jelinek from comment #28)
> As long as the scale is a power of two or 1.0 / power of two, I don't see
> why any version wouldn't be inaccurate.

Yes, but the constant scale_up is incorrectly selected.
scale_up = std::exp2(Type(limits::max_exponent-1)) --> ok
scale_up = std::exp2(Type(limits::max_exponent/2)) --> error
scale_up = prev_power2(sqrt_max) --> error

scale_down = std::exp2(Type(limits::min_exponent-1))
also seems to me to be more favorable.

PS:
There seems to be a problem with random numbers and std::float16_t, which is
why I use std::uniform_real_distribution<std::float128_t>. I have not yet found
out exactly where the error lies.

thx
Gero


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

        constexpr Type
                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);

        int
                exp;

        z = std::frexp(z, &exp);
        y = std::ldexp(y, -exp);
        x = std::ldexp(x, -exp);
        return std::ldexp(std::sqrt(__builtin_assoc_barrier(x*x + y*y) + z*z),
exp);
}

template <std::floating_point 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-1)),
                scale_down      = std::exp2(Type(limits::min_exponent-1)),
                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 (const bool b{z>=sqrt_min}; b && z<=sqrt_max) [[likely]]
        {
                //      no scale
                return std::sqrt(__builtin_assoc_barrier(x*x + y*y) + z*z);
        }
        else
        {
                const Type
                        scale = b ? scale_down : scale_up;

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

template <std::floating_point Type>
void    test(const size_t count, const Type min, const Type max, const Type
factor)
{
        std::random_device rd{};
        std::mt19937 gen{rd()};
        std::uniform_real_distribution<std::float128_t> dis{min, max};

        auto rnd = [&]() noexcept -> Type { return Type(dis(gen) * factor); };

        for (size_t i=0; i<count; ++i)
        {
                const Type
                        x = rnd(),
                        y = rnd(),
                        z = rnd(),
                        r0= hypot_exp(x, y, z),
                        r1= hypot_gp(x, y, z);

                if (r0 != r1) [[unlikely]]
                {
                        std::cout << "error\n";
                        return;
                }
        }
        std::cout << "ok\n";
}

int main()
{
        using value_type = std::float64_t;
        using limits = std::numeric_limits<value_type>;

        test<value_type>(1024*1024, 0.5, 1, limits::max());
        test<value_type>(1024*1024, 0, 1, limits::min());

        return EXIT_SUCCESS;
}