This adds the new 3D std::hypot() functions. This implementation seems
to be faster than the naïve sqrt(x*x + y*y + z*z) implementation, or
hypot(hypot(x, y), z), and should be a bit more accurate at very large
or very small values due to reducing the arguments by the largest one.
Improvements welcome though, as this is not my forte.
The test might not be very good, but tests some small integer values
and some other values where accuracy is lost for one or other of the
alternative implementations mentioned above. If this FAILs for some
32-bit targets we might need to adjust the tolerances or the
dg-options.
* doc/xml/manual/status_cxx2017.xml: Update status.
* include/c_global/cmath (hypot): Add three-dimensional overloads.
* testsuite/26_numerics/headers/cmath/hypot.cc: New.
Tested powerpc64le-linux and x86_64-linux, committed to trunk.
commit 33fb40d1445fbc711d64a1ff8b9a0a8f092a2e7e
Author: Jonathan Wakely <jwak...@redhat.com>
Date: Tue Sep 27 15:44:45 2016 +0100
Define 3-argument overloads of std::hypot for C++17 (P0030R1)
* doc/xml/manual/status_cxx2017.xml: Update status.
* include/c_global/cmath (hypot): Add three-dimensional overloads.
* testsuite/26_numerics/headers/cmath/hypot.cc: New.
diff --git a/libstdc++-v3/doc/xml/manual/status_cxx2017.xml
b/libstdc++-v3/doc/xml/manual/status_cxx2017.xml
index 76eaaa0..4ead6b9 100644
--- a/libstdc++-v3/doc/xml/manual/status_cxx2017.xml
+++ b/libstdc++-v3/doc/xml/manual/status_cxx2017.xml
@@ -633,14 +633,13 @@ Feature-testing recommendations for C++</link>.
</row>
<row>
- <?dbhtml bgcolor="#C8B0B0" ?>
<entry> Proposal to Introduce a 3-Argument Overload to
<code>std::hypot</code> </entry>
<entry>
<link xmlns:xlink="http://www.w3.org/1999/xlink";
xlink:href="http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2015/p0030r1.pdf";>
P0030R1
</link>
</entry>
- <entry align="center"> No </entry>
+ <entry align="center"> 7 </entry>
<entry><code> __cpp_lib_hypot >= 201603 </code></entry>
</row>
diff --git a/libstdc++-v3/include/c_global/cmath
b/libstdc++-v3/include/c_global/cmath
index 6db9dee..fffa0e7 100644
--- a/libstdc++-v3/include/c_global/cmath
+++ b/libstdc++-v3/include/c_global/cmath
@@ -1455,6 +1455,46 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
return hypot(__type(__x), __type(__y));
}
+#if __cplusplus > 201402L
+#define __cpp_lib_hypot 201603
+ // [c.math.hypot3], three-dimensional hypotenuse
+
+ template<typename _Tp>
+ inline _Tp
+ __hypot3(_Tp __x, _Tp __y, _Tp __z)
+ {
+ __x = std::abs(__x);
+ __y = std::abs(__y);
+ __z = std::abs(__z);
+ if (_Tp __a = __x < __y ? __y < __z ? __z : __y : __x < __z ? __z : __x)
+ return __a * std::sqrt((__x / __a) * (__x / __a)
+ + (__y / __a) * (__y / __a)
+ + (__z / __a) * (__z / __a));
+ else
+ return {};
+ }
+
+ inline float
+ hypot(float __x, float __y, float __z)
+ { return std::__hypot3<float>(__x, __y, __z); }
+
+ inline double
+ hypot(double __x, double __y, double __z)
+ { return std::__hypot3<double>(__x, __y, __z); }
+
+ inline long double
+ hypot(long double __x, long double __y, long double __z)
+ { return std::__hypot3<long double>(__x, __y, __z); }
+
+ template<typename _Tp, typename _Up, typename _Vp>
+ typename __gnu_cxx::__promote_3<_Tp, _Up, _Vp>::__type
+ hypot(_Tp __x, _Up __y, _Vp __z)
+ {
+ using __type = typename __gnu_cxx::__promote_3<_Tp, _Up, _Vp>::__type;
+ return std::__hypot3<__type>(__x, __y, __z);
+ }
+#endif // C++17
+
#ifndef __CORRECT_ISO_CPP11_MATH_H_PROTO
constexpr int
ilogb(float __x)
diff --git a/libstdc++-v3/testsuite/26_numerics/headers/cmath/hypot.cc
b/libstdc++-v3/testsuite/26_numerics/headers/cmath/hypot.cc
new file mode 100644
index 0000000..4a6841c
--- /dev/null
+++ b/libstdc++-v3/testsuite/26_numerics/headers/cmath/hypot.cc
@@ -0,0 +1,138 @@
+// Copyright (C) 2016 Free Software Foundation, Inc.
+//
+// This file is part of the GNU ISO C++ Library. This library is free
+// software; you can redistribute it and/or modify it under the
+// terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 3, or (at your option)
+// any later version.
+
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+
+// You should have received a copy of the GNU General Public License along
+// with this library; see the file COPYING3. If not see
+// <http://www.gnu.org/licenses/>.
+
+// { dg-options "-std=gnu++17" }
+// { dg-do run { target c++1z } }
+
+#include <cmath>
+#include <type_traits>
+#if defined(__TEST_DEBUG)
+#include <iostream>
+#define VERIFY(A) \
+if (!(A)) \
+ { \
+ std::cout << "line " << __LINE__ \
+ << " max_abs_frac = " << max_abs_frac \
+ << " tolerance = " << toler \
+ << std::endl; \
+ }
+#else
+#include <testsuite_hooks.h>
+#endif
+
+using std::is_same_v;
+static_assert(is_same_v<double, decltype(std::hypot(0.0, 0.0, 0.0))>);
+static_assert(is_same_v<double, decltype(std::hypot(0.0f, 0.0, 0.0))>);
+static_assert(is_same_v<double, decltype(std::hypot(0.0, 0.0f, 0.0))>);
+static_assert(is_same_v<double, decltype(std::hypot(0.0, 0.0, 0.0f))>);
+static_assert(is_same_v<double, decltype(std::hypot(0.0f, 0.0f, 0.0))>);
+static_assert(is_same_v<double, decltype(std::hypot(0.0f, 0.0, 0))>);
+static_assert(is_same_v<long double, decltype(std::hypot(0.0f, 0.0, 0.0l))>);
+static_assert(is_same_v<long double, decltype(std::hypot(0, 0.0, 0.0l))>);
+
+template<typename T> struct testcase_hypot { T x, y, z, f0; };
+
+template<typename Tp, unsigned int Num>
+ void
+ test(const testcase_hypot<Tp> (&data)[Num], Tp toler)
+ {
+ bool test __attribute__((unused)) = true;
+ const Tp eps = std::numeric_limits<Tp>::epsilon();
+ Tp max_abs_diff = -Tp(1);
+ Tp max_abs_frac = -Tp(1);
+ unsigned int num_datum = Num;
+ for (unsigned int i = 0; i < num_datum; ++i)
+ {
+ const Tp f = std::hypot(data[i].x, data[i].y, data[i].z);
+ const Tp f0 = data[i].f0;
+ const Tp diff = f - f0;
+ if (std::abs(diff) > max_abs_diff)
+ max_abs_diff = std::abs(diff);
+ if (std::abs(f0) > Tp(10) * eps && std::abs(f) > Tp(10) * eps)
+ {
+ const Tp frac = diff / f0;
+ if (std::abs(frac) > max_abs_frac)
+ max_abs_frac = std::abs(frac);
+ }
+ }
+ VERIFY(max_abs_frac < toler);
+ }
+
+const testcase_hypot<double> data1[] = {
+ { 0.0, 0.0, 0.0, 0.0 },
+ { 0.0, 1.0, 1.0, std::sqrt(2.0) },
+ { 1.0, 1.0, 1.0, std::sqrt(3.0) },
+ { 1.0, 2.0, 2.0, 3.0 },
+ { 2.0, 3.0, 6.0, 7.0 },
+ { 1.0, 4.0, 8.0, 9.0 },
+ { 4.0, 4.0, 7.0, 9.0 },
+ { 12.0, 16.0, 21.0, 29.0 },
+ { 1e8, 1., 1e-8, 1e8 },
+ { 1., 1e8, 1e-8, 1e8 },
+ { 1e-8, 1., 1e8, 1e8 },
+ { 1e-2, 1e-4, 1e-4, 0.01000099995 },
+ { 214748364., 214748364., 214748364., 371955077.2902952 }
+};
+const double toler1 = 1e-12;
+
+const testcase_hypot<float> data2[] = {
+ { 0.0f, 0.0f, 0.0f, 0.0f },
+ { 0.0f, 1.0f, 1.0f, std::sqrt(2.0f) },
+ { 1.0f, 1.0f, 1.0f, std::sqrt(3.0f) },
+ { 1.0f, 2.0f, 2.0f, 3.0f },
+ { 2.0f, 3.0f, 6.0f, 7.0f },
+ { 1.0f, 4.0f, 8.0f, 9.0f },
+ { 4.0f, 4.0f, 7.0f, 9.0f },
+ { 12.0f, 16.0f, 21.0f, 29.0f },
+ { 1e8f, 1.f, 1e-8f, 1e8f },
+ { 1.f, 1e8f, 1e-8f, 1e8f },
+ { 1e-8f, 1.f, 1e8f, 1e8f },
+ { 1e-2f, 1e-4f, 1e-4f, 0.010001f },
+ { 214748364.f, 214748364.f, 214748364.f, 371955072.f }
+};
+const float toler2 = 1e-7f;
+
+const testcase_hypot<long double> data3[] = {
+ { 0.0l, 0.0l, 0.0l, 0.0l },
+ { 0.0l, 1.0l, 1.0l, std::sqrt(2.0l) },
+ { 1.0l, 1.0l, 1.0l, std::sqrt(3.0l) },
+ { 1.0l, 2.0l, 2.0l, 3.0l },
+ { 2.0l, 3.0l, 6.0l, 7.0l },
+ { 1.0l, 4.0l, 8.0l, 9.0l },
+ { 4.0l, 4.0l, 7.0l, 9.0l },
+ { 12.0l, 16.0l, 21.0l, 29.0l },
+ { 1e8l, 1.l, 1e-8l, 1e8l },
+ { 1.l, 1e8l, 1e-8l, 1e8l },
+ { 1e-8l, 1.l, 1e8l, 1e8l },
+ { 1e-2l, 1e-4l, 1e-4l, 0.010000999950004999375l },
+ { 2147483647.l, 2147483647.l, 2147483647.l, 3719550785.027307813987l }
+};
+const long double toler3 = 1e-16l;
+
+void
+test01()
+{
+ test(data1, toler1);
+ test(data2, toler2);
+ test(data3, toler3);
+}
+
+int
+main()
+{
+ test01();
+}
