From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: by sourceware.org (Postfix, from userid 48) id 4E3843858C66; Mon, 4 Mar 2024 11:12:29 +0000 (GMT) DKIM-Filter: OpenDKIM Filter v2.11.0 sourceware.org 4E3843858C66 DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gcc.gnu.org; s=default; t=1709550749; bh=2Ud1B994FfRn6VxdySZh+ZAzRC/yxESWUKAp09FKcbc=; h=From:To:Subject:Date:In-Reply-To:References:From; b=yIsDrQyo1eJx/4klZhwdVst2eGRbFutCm+GnsxYvYvld4+NIUvTgjqd5zy9PaU/BF IJgomhkeBneg+atiNY6yVDIEQW3Bc56xtPCUjmknx41qgbE8xnyAbVF6ll4nMWZMJx lTUL0K5yrO21yDIsf5t0AyKKGHVNBlAMaNcGzIbw= From: "jakub at gcc dot gnu.org" To: gcc-bugs@gcc.gnu.org Subject: [Bug libstdc++/77776] C++17 std::hypot implementation is poor Date: Mon, 04 Mar 2024 11:12:28 +0000 X-Bugzilla-Reason: CC X-Bugzilla-Type: changed X-Bugzilla-Watch-Reason: None X-Bugzilla-Product: gcc X-Bugzilla-Component: libstdc++ X-Bugzilla-Version: 7.0 X-Bugzilla-Keywords: X-Bugzilla-Severity: normal X-Bugzilla-Who: jakub at gcc dot gnu.org X-Bugzilla-Status: ASSIGNED X-Bugzilla-Resolution: X-Bugzilla-Priority: P3 X-Bugzilla-Assigned-To: emsr at gcc dot gnu.org X-Bugzilla-Target-Milestone: --- X-Bugzilla-Flags: X-Bugzilla-Changed-Fields: cc Message-ID: In-Reply-To: References: Content-Type: text/plain; charset="UTF-8" Content-Transfer-Encoding: quoted-printable X-Bugzilla-URL: http://gcc.gnu.org/bugzilla/ Auto-Submitted: auto-generated MIME-Version: 1.0 List-Id: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=3D77776 Jakub Jelinek changed: What |Removed |Added ---------------------------------------------------------------------------- CC| |jakub at gcc dot gnu.org --- Comment #16 from Jakub Jelinek --- (In reply to Matthias Kretz (Vir) from comment #15) > Your implementation still needs to solve: >=20 > 1. Loss of precision because of division & subsequent scaling by max. Use= rs > comparing std::hypot(x, y, z) against a simple std::sqrt(x * x + y * y + = z * > z) might wonder why they want to use std::hypot if it's less precise. >=20 > 2. Relatively high cost (in latency and throughput) because of the three > divisions. You could replace it with scale =3D 1/max; x *=3D scale; ... B= ut that > can reduce precision even further. >=20 > 3. Summation of the x, y, and z squares isn't associative if you care abo= ut > precision. A high quality implementation needs to add the two lowest valu= es > first. >=20 > Here's a precise but inefficient implementation: > (https://compiler-explorer.com/z/ocGPnsYE3) >=20 > template > [[gnu::optimize("-fno-unsafe-math-optimizations")]] > T > hypot3(T x, T y, T z) > { > x =3D std::abs(x); > y =3D std::abs(y); > z =3D std::abs(z); > if (std::isinf(x) || std::isinf(y) || std::isinf(z)) > return std::__infinity_v; > else if (std::isnan(x) || std::isnan(y) || std::isnan(z)) > return std::__quiet_NaN_v; > else if (x =3D=3D y && y =3D=3D z) > return x * std::sqrt(T(3)); > else if (z =3D=3D 0 && y =3D=3D 0) > return x; > else if (x =3D=3D 0 && z =3D=3D 0) > return y; > else if (x =3D=3D 0 && y =3D=3D 0) > return z; > else > { > T hi =3D std::max(std::max(x, y), z); > T lo0 =3D std::min(std::max(x, y), z); > T lo1 =3D std::min(x, y); > int e =3D 0; > hi =3D std::frexp(hi, &e); > lo0 =3D std::ldexp(lo0, -e); > lo1 =3D std::ldexp(lo1, -e); > T lo =3D lo0 * lo0 + lo1 * lo1; > return std::ldexp(std::sqrt(hi * hi + lo), e); > } > } >=20 > AFAIK > https://gcc.gnu.org/git/?p=3Dgcc.git;a=3Dblob;f=3Dlibstdc%2B%2B-v3/includ= e/ > experimental/bits/simd_math.h;h=3D06e7b4496f9917f886f66fbd7629700dd17e55f= 9; > hb=3DHEAD#l1168 is a precise and efficient implementation. It also avoids > division altogether unless an input is subnormal. What glibc does there for the 2 argument hypot is after handling the non-fi= nite cases finds the minimum and maximum and uses just normal multiplication, addition + sqrt for the common case where maximum isn't too large and minim= um isn't too small. 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.=