public inbox for gcc-bugs@sourceware.org
help / color / mirror / Atom feed
From: "g.peterhoff@t-online.de" <gcc-bugzilla@gcc.gnu.org>
To: gcc-bugs@gcc.gnu.org
Subject: [Bug libstdc++/77776] C++17 std::hypot implementation is poor
Date: Wed, 10 Apr 2024 18:08:18 +0000 [thread overview]
Message-ID: <bug-77776-4-Gt6hPS09CM@http.gcc.gnu.org/bugzilla/> (raw)
In-Reply-To: <bug-77776-4@http.gcc.gnu.org/bugzilla/>
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=77776
--- Comment #29 from g.peterhoff@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;
}
prev parent reply other threads:[~2024-04-10 18:08 UTC|newest]
Thread overview: 20+ messages / expand[flat|nested] mbox.gz Atom feed top
[not found] <bug-77776-4@http.gcc.gnu.org/bugzilla/>
2021-05-04 12:31 ` rguenth at gcc dot gnu.org
2024-02-29 18:23 ` g.peterhoff@t-online.de
2024-03-01 12:27 ` redi at gcc dot gnu.org
2024-03-02 15:52 ` g.peterhoff@t-online.de
2024-03-04 3:49 ` de34 at live dot cn
2024-03-04 10:45 ` mkretz at gcc dot gnu.org
2024-03-04 11:12 ` jakub at gcc dot gnu.org
2024-03-04 17:14 ` mkretz at gcc dot gnu.org
2024-03-04 20:21 ` jakub at gcc dot gnu.org
2024-03-06 2:47 ` g.peterhoff@t-online.de
2024-03-06 9:43 ` mkretz at gcc dot gnu.org
2024-03-06 12:35 ` redi at gcc dot gnu.org
2024-03-12 11:31 ` mkretz at gcc dot gnu.org
2024-03-19 0:09 ` g.peterhoff@t-online.de
2024-03-25 10:06 ` mkretz at gcc dot gnu.org
2024-04-05 1:08 ` g.peterhoff@t-online.de
2024-04-05 1:23 ` g.peterhoff@t-online.de
2024-04-10 16:41 ` g.peterhoff@t-online.de
2024-04-10 16:48 ` jakub at gcc dot gnu.org
2024-04-10 18:08 ` g.peterhoff@t-online.de [this message]
Reply instructions:
You may reply publicly to this message via plain-text email
using any one of the following methods:
* Save the following mbox file, import it into your mail client,
and reply-to-all from there: mbox
Avoid top-posting and favor interleaved quoting:
https://en.wikipedia.org/wiki/Posting_style#Interleaved_style
* Reply using the --to, --cc, and --in-reply-to
switches of git-send-email(1):
git send-email \
--in-reply-to=bug-77776-4-Gt6hPS09CM@http.gcc.gnu.org/bugzilla/ \
--to=gcc-bugzilla@gcc.gnu.org \
--cc=gcc-bugs@gcc.gnu.org \
/path/to/YOUR_REPLY
https://kernel.org/pub/software/scm/git/docs/git-send-email.html
* If your mail client supports setting the In-Reply-To header
via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line
before the message body.
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).