From mboxrd@z Thu Jan  1 00:00:00 1970
Return-Path: <gcc-bugzilla@gcc.gnu.org>
Received: by sourceware.org (Postfix, from userid 48)
	id A3DC7385843B; Mon,  4 Mar 2024 10:45:35 +0000 (GMT)
DKIM-Filter: OpenDKIM Filter v2.11.0 sourceware.org A3DC7385843B
DKIM-Signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gcc.gnu.org;
	s=default; t=1709549135;
	bh=XHyrRSgjoucfu9IvyIVURsHMPIKVB5nIFdD/AYnNNIs=;
	h=From:To:Subject:Date:In-Reply-To:References:From;
	b=rIdABfDjdH3VlrjM2eIp0rqVnqGrgeIzjuwN73JsMhrG4JRfpyWzdNl4UfqXb6597
	 L16rEZm2jc4ilg+s7qYZ6V6DZevp8lnkhhXKYFcsLHsjPCHkHkvMC1duJTFIqWscYE
	 0iCH7rvvzTVKM3EEPc9NEPeW9g0nddoQRmz+4yn0=
From: "mkretz at gcc dot gnu.org" <gcc-bugzilla@gcc.gnu.org>
To: gcc-bugs@gcc.gnu.org
Subject: [Bug libstdc++/77776] C++17 std::hypot implementation is poor
Date: Mon, 04 Mar 2024 10:45:31 +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: mkretz 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: 
Message-ID: <bug-77776-4-NbvVrQMtS2@http.gcc.gnu.org/bugzilla/>
In-Reply-To: <bug-77776-4@http.gcc.gnu.org/bugzilla/>
References: <bug-77776-4@http.gcc.gnu.org/bugzilla/>
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: <gcc-bugs.sourceware.org>

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

--- Comment #15 from Matthias Kretz (Vir) <mkretz at gcc dot gnu.org> ---
Your implementation still needs to solve:

1. Loss of precision because of division & subsequent scaling by max. Users
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.

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; ... But=
 that
can reduce precision even further.

3. Summation of the x, y, and z squares isn't associative if you care about
precision. A high quality implementation needs to add the two lowest values
first.

Here's a precise but inefficient implementation:
(https://compiler-explorer.com/z/ocGPnsYE3)

template <typename T>
[[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<T>;
  else if (std::isnan(x) || std::isnan(y) || std::isnan(z))
    return std::__quiet_NaN_v<T>;
  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);
    }
}

AFAIK
https://gcc.gnu.org/git/?p=3Dgcc.git;a=3Dblob;f=3Dlibstdc%2B%2B-v3/include/=
experimental/bits/simd_math.h;h=3D06e7b4496f9917f886f66fbd7629700dd17e55f9;=
hb=3DHEAD#l1168
is a precise and efficient implementation. It also avoids division altogeth=
er
unless an input is subnormal.=