[libcxx-commits] [libcxx] [libc++][math] Mathematical special functions: Implementing `std::legendre` (PR #92175)
via libcxx-commits
libcxx-commits at lists.llvm.org
Tue May 14 13:55:04 PDT 2024
https://github.com/PaulXiCao created https://github.com/llvm/llvm-project/pull/92175
Implementing the Legendre polynomials `std::legendre, std::legendref, std::legendrel` which are part of C++17's mathematical special functions.
I started out from this abandoned merge request: https://reviews.llvm.org/D58876.
Documentation in `libcxx/docs/Status/` is not yet included. This will be added after the integration of #89982.
>From e6ca851c3832b9a33c8386ad76676cebe7567143 Mon Sep 17 00:00:00 2001
From: Paul <{ID}+{username}@users.noreply.github.com>
Date: Sun, 12 May 2024 21:48:25 +0200
Subject: [PATCH 1/2] implementation + unit
---
libcxx/include/CMakeLists.txt | 1 +
libcxx/include/__math/special_functions.h | 86 ++++++++++++++
libcxx/include/cmath | 9 ++
libcxx/modules/std/cmath.inc | 2 +
.../std/numerics/c.math/legendre.pass.cpp | 111 ++++++++++++++++++
5 files changed, 209 insertions(+)
create mode 100644 libcxx/include/__math/special_functions.h
create mode 100644 libcxx/test/std/numerics/c.math/legendre.pass.cpp
diff --git a/libcxx/include/CMakeLists.txt b/libcxx/include/CMakeLists.txt
index 01e9c247560ca..02c46d87759d4 100644
--- a/libcxx/include/CMakeLists.txt
+++ b/libcxx/include/CMakeLists.txt
@@ -516,6 +516,7 @@ set(files
__math/remainder.h
__math/roots.h
__math/rounding_functions.h
+ __math/special_functions.h
__math/traits.h
__math/trigonometric_functions.h
__mbstate_t.h
diff --git a/libcxx/include/__math/special_functions.h b/libcxx/include/__math/special_functions.h
new file mode 100644
index 0000000000000..69e4e0e34bbff
--- /dev/null
+++ b/libcxx/include/__math/special_functions.h
@@ -0,0 +1,86 @@
+// -*- C++ -*-
+//===----------------------------------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef _LIBCPP___MATH_SPECIAL_FUNCTIONS_H
+#define _LIBCPP___MATH_SPECIAL_FUNCTIONS_H
+
+#include <__config>
+#include <__math/abs.h>
+#include <__math/copysign.h>
+#include <__math/traits.h>
+#include <__type_traits/enable_if.h>
+#include <__type_traits/is_integral.h>
+#include <limits>
+#include <stdexcept>
+
+#if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
+# pragma GCC system_header
+#endif
+
+_LIBCPP_BEGIN_NAMESPACE_STD
+
+#if _LIBCPP_STD_VER >= 17
+
+template <class _Real>
+_LIBCPP_HIDE_FROM_ABI _Real __legendre(unsigned __n, _Real __x) {
+ // The Legendre polynomial H_n(x).
+ // The implementation is based on the recurrence formula: (n+1) P_{n+1}(x) = (2n+1) x P_n(x) - n P_{n-1}(x).
+ // Press, William H., et al. Numerical recipes 3rd edition: The art of scientific computing.
+ // Cambridge university press, 2007, p. 183.
+
+ // NOLINTBEGIN(readability-identifier-naming)
+ if (__math::isnan(__x))
+ return __x;
+
+ if (__math::fabs(__x) > 1)
+ throw std::domain_error("Argument `x` of Legendre function is out of range: `|x| <= 1`.");
+
+ _Real __P_0{1};
+ if (__n == 0)
+ return __P_0;
+
+ _Real __P_n_prev = __P_0;
+ _Real __P_n = __x;
+ for (unsigned __i = 1; __i < __n; ++__i) {
+ _Real __P_n_next = ((2 * __i + 1) * __x * __P_n - __i * __P_n_prev) / static_cast<_Real>(__i + 1);
+ __P_n_prev = __P_n;
+ __P_n = __P_n_next;
+ }
+
+ if (!__math::isfinite(__P_n)) {
+ // Overflow occured. Two possible cases:
+ // n is odd: return infinity of the same sign as x.
+ // n is even: return +Inf
+ _Real __inf = std::numeric_limits<_Real>::infinity();
+ return (__n & 1) ? __math::copysign(__inf, __x) : __inf;
+ }
+ return __P_n;
+ // NOLINTEND(readability-identifier-naming)
+}
+
+inline _LIBCPP_HIDE_FROM_ABI double legendre(unsigned __n, double __x) { return std::__legendre(__n, __x); }
+
+inline _LIBCPP_HIDE_FROM_ABI float legendre(unsigned __n, float __x) { return std::__legendre(__n, __x); }
+
+inline _LIBCPP_HIDE_FROM_ABI long double legendre(unsigned __n, long double __x) { return std::__legendre(__n, __x); }
+
+inline _LIBCPP_HIDE_FROM_ABI float legendref(unsigned __n, float __x) { return std::legendre(__n, __x); }
+
+inline _LIBCPP_HIDE_FROM_ABI long double legendrel(unsigned __n, long double __x) { return std::legendre(__n, __x); }
+
+template <class _Integer, std::enable_if_t<std::is_integral_v<_Integer>, int> = 0>
+_LIBCPP_HIDE_FROM_ABI double legendre(unsigned __n, _Integer __x) {
+ return std::legendre(__n, static_cast<double>(__x));
+}
+
+#endif // _LIBCPP_STD_VER >= 17
+
+_LIBCPP_END_NAMESPACE_STD
+
+#endif // _LIBCPP___MATH_SPECIAL_FUNCTIONS_H
diff --git a/libcxx/include/cmath b/libcxx/include/cmath
index dd194bbb55896..cffd5efdaf7a2 100644
--- a/libcxx/include/cmath
+++ b/libcxx/include/cmath
@@ -216,6 +216,14 @@ int ilogb (arithmetic x);
int ilogbf(float x);
int ilogbl(long double x);
+double legendre(unsigned n, double x); // C++17
+float legendre(unsigned n, float x); // C++17
+long double legendre(unsigned n, long double x); // C++17
+float legendref(unsigned n, float x); // C++17
+long double legendrel(unsigned n, long double x); // C++17
+template <class Integer>
+double legendre(unsigned n, Integer x); // C++17
+
floating_point lgamma (arithmetic x);
float lgammaf(float x);
long double lgammal(long double x);
@@ -315,6 +323,7 @@ constexpr long double lerp(long double a, long double b, long double t) noexcept
#include <limits>
#include <version>
+#include <__math/special_functions.h>
#include <math.h>
#ifndef _LIBCPP_MATH_H
diff --git a/libcxx/modules/std/cmath.inc b/libcxx/modules/std/cmath.inc
index a463c1e3ccf86..74b3477e7e880 100644
--- a/libcxx/modules/std/cmath.inc
+++ b/libcxx/modules/std/cmath.inc
@@ -344,12 +344,14 @@ export namespace std {
using std::laguerre;
using std::laguerref;
using std::laguerrel;
+#endif
// [sf.cmath.legendre], Legendre polynomials
using std::legendre;
using std::legendref;
using std::legendrel;
+#if 0
// [sf.cmath.riemann.zeta], Riemann zeta function
using std::riemann_zeta;
using std::riemann_zetaf;
diff --git a/libcxx/test/std/numerics/c.math/legendre.pass.cpp b/libcxx/test/std/numerics/c.math/legendre.pass.cpp
new file mode 100644
index 0000000000000..06a3e731c844e
--- /dev/null
+++ b/libcxx/test/std/numerics/c.math/legendre.pass.cpp
@@ -0,0 +1,111 @@
+//===----------------------------------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+// UNSUPPORTED: c++03, c++11, c++14
+
+// <cmath>
+
+// double legendre(unsigned n, double x);
+// float legendre(unsigned n, float x);
+// long double legendre(unsigned n, long double x);
+// float legendref(unsigned n, float x);
+// long double legendrel(unsigned n, long double x);
+// template <class Integer>
+// double legendre(unsigned n, Integer x);
+
+#include <cassert>
+#include <cmath>
+#include <limits>
+
+template <class T>
+void testLegendreNaNPropagation() {
+ const unsigned MaxN = 127;
+ const T x = std::numeric_limits<T>::quiet_NaN();
+ for (unsigned n = 0; n <= MaxN; ++n) {
+ assert(std::isnan(std::legendre(n, x)));
+ }
+}
+
+template <class T>
+void testLegendreNotNaN(const T x) {
+ assert(!std::isnan(x));
+ const unsigned MaxN = 127;
+ for (unsigned n = 0; n <= MaxN; ++n) {
+ assert(!std::isnan(std::legendre(n, x)));
+ }
+}
+
+template <class T>
+void testLegendreThrows(const T x) {
+#ifndef _LIBCPP_NO_EXCEPTIONS
+ const unsigned MaxN = 127;
+ for (unsigned n = 0; n <= MaxN; ++n) {
+ bool Throws = false;
+ try {
+ std::legendre(n, x);
+ } catch (const std::domain_error&) {
+ Throws = true;
+ }
+ assert(Throws);
+ }
+#endif // _LIBCPP_NO_EXCEPTIONS
+}
+
+template <class T>
+void testLegendreAnalytic(const T x, const T AbsTolerance, const T RelTolerance) {
+ assert(!std::isnan(x));
+ const auto compareFloatingPoint = [AbsTolerance, RelTolerance](const T Result, const T ExpectedResult) {
+ if (std::isinf(ExpectedResult) && std::isinf(Result))
+ return true;
+
+ if (std::isnan(ExpectedResult) || std::isnan(Result))
+ return false;
+
+ const T Tolerance = AbsTolerance + std::abs(ExpectedResult) * RelTolerance;
+ return std::abs(Result - ExpectedResult) < Tolerance;
+ };
+
+ const auto l0 = [](T) { return T(1); };
+ const auto l1 = [](T y) { return y; };
+ const auto l2 = [](T y) { return (T(3) * y * y - T(1)) / T(2); };
+ const auto l3 = [](T y) { return (T(5) * y * y - T(3)) * y / T(2); };
+ const auto l4 = [](T y) { return (T(35) * y * y * y * y - T(30) * y * y + T(3)) / T(8); };
+ const auto l5 = [](T y) { return (T(63) * y * y * y * y - T(70) * y * y + T(15)) * y / T(8); };
+ const auto l6 = [](T y) {
+ const T y2 = y * y;
+ return (T(231) * y2 * y2 * y2 - T(315) * y2 * y2 + T(105) * y2 - T(5)) / T(16);
+ };
+
+ assert(compareFloatingPoint(std::legendre(0, x), l0(x)));
+ assert(compareFloatingPoint(std::legendre(1, x), l1(x)));
+ assert(compareFloatingPoint(std::legendre(2, x), l2(x)));
+ assert(compareFloatingPoint(std::legendre(3, x), l3(x)));
+ assert(compareFloatingPoint(std::legendre(4, x), l4(x)));
+ assert(compareFloatingPoint(std::legendre(5, x), l5(x)));
+ assert(compareFloatingPoint(std::legendre(6, x), l6(x)));
+}
+
+template <class T>
+void testLegendre(const T AbsTolerance, const T RelTolerance) {
+ testLegendreNaNPropagation<T>();
+ testLegendreThrows<T>(T(-5));
+ testLegendreThrows<T>(T(5));
+
+ const T Samples[] = {T(-1.0), T(-0.5), T(-0.1), T(0.0), T(0.1), T(0.5), T(1.0)};
+
+ for (T x : Samples) {
+ testLegendreNotNaN(x);
+ testLegendreAnalytic(x, AbsTolerance, RelTolerance);
+ }
+}
+
+int main() {
+ testLegendre<float>(1e-6f, 1e-6f);
+ testLegendre<double>(1e-9, 1e-9);
+ testLegendre<long double>(1e-12, 1e-12);
+}
>From 80a845aa08e9b39c08b885da6212e9cf9a5e10d5 Mon Sep 17 00:00:00 2001
From: Paul <{ID}+{username}@users.noreply.github.com>
Date: Sun, 12 May 2024 23:29:51 +0200
Subject: [PATCH 2/2] structured tests; fixed implementation
---
libcxx/include/__math/special_functions.h | 9 +-
.../std/numerics/c.math/legendre.pass.cpp | 206 ++++++++++++------
2 files changed, 138 insertions(+), 77 deletions(-)
diff --git a/libcxx/include/__math/special_functions.h b/libcxx/include/__math/special_functions.h
index 69e4e0e34bbff..e507b243dbf03 100644
--- a/libcxx/include/__math/special_functions.h
+++ b/libcxx/include/__math/special_functions.h
@@ -39,7 +39,7 @@ _LIBCPP_HIDE_FROM_ABI _Real __legendre(unsigned __n, _Real __x) {
return __x;
if (__math::fabs(__x) > 1)
- throw std::domain_error("Argument `x` of Legendre function is out of range: `|x| <= 1`.");
+ std::__throw_domain_error("Argument `x` of Legendre function is out of range: `|x| <= 1`.");
_Real __P_0{1};
if (__n == 0)
@@ -53,13 +53,6 @@ _LIBCPP_HIDE_FROM_ABI _Real __legendre(unsigned __n, _Real __x) {
__P_n = __P_n_next;
}
- if (!__math::isfinite(__P_n)) {
- // Overflow occured. Two possible cases:
- // n is odd: return infinity of the same sign as x.
- // n is even: return +Inf
- _Real __inf = std::numeric_limits<_Real>::infinity();
- return (__n & 1) ? __math::copysign(__inf, __x) : __inf;
- }
return __P_n;
// NOLINTEND(readability-identifier-naming)
}
diff --git a/libcxx/test/std/numerics/c.math/legendre.pass.cpp b/libcxx/test/std/numerics/c.math/legendre.pass.cpp
index 06a3e731c844e..aee7850f0a5e8 100644
--- a/libcxx/test/std/numerics/c.math/legendre.pass.cpp
+++ b/libcxx/test/std/numerics/c.math/legendre.pass.cpp
@@ -18,94 +18,162 @@
// template <class Integer>
// double legendre(unsigned n, Integer x);
+#include <array>
#include <cassert>
#include <cmath>
#include <limits>
+#include "type_algorithms.h"
+
+inline constexpr unsigned g_max_n = 128;
+
template <class T>
-void testLegendreNaNPropagation() {
- const unsigned MaxN = 127;
- const T x = std::numeric_limits<T>::quiet_NaN();
- for (unsigned n = 0; n <= MaxN; ++n) {
- assert(std::isnan(std::legendre(n, x)));
- }
+std::array<T, 7> sample_points() {
+ return {-1.0, -0.5, -0.1, 0.0, 0.1, 0.5, 1.0};
}
-template <class T>
-void testLegendreNotNaN(const T x) {
- assert(!std::isnan(x));
- const unsigned MaxN = 127;
- for (unsigned n = 0; n <= MaxN; ++n) {
- assert(!std::isnan(std::legendre(n, x)));
+template <class Real>
+class CompareFloatingValues {
+private:
+ Real tol;
+
+public:
+ CompareFloatingValues() {
+ if (std::is_same_v<Real, float>)
+ tol = 1e-6f;
+ else if (std::is_same_v<Real, double>)
+ tol = 1e-9;
+ else
+ tol = 1e-12l;
}
-}
-template <class T>
-void testLegendreThrows(const T x) {
-#ifndef _LIBCPP_NO_EXCEPTIONS
- const unsigned MaxN = 127;
- for (unsigned n = 0; n <= MaxN; ++n) {
- bool Throws = false;
- try {
- std::legendre(n, x);
- } catch (const std::domain_error&) {
- Throws = true;
- }
- assert(Throws);
+ bool operator()(Real result, Real expected) const {
+ if (std::isinf(expected) && std::isinf(result))
+ return result == expected;
+
+ if (std::isnan(expected) || std::isnan(result))
+ return false;
+
+ return std::abs(result - expected) < (tol + std::abs(expected) * tol);
+ }
+};
+
+template <class Real>
+void test() {
+ { // checks if NaNs are reported correctly (i.e. output == input for input == NaN)
+ using nl = std::numeric_limits<Real>;
+ for (Real NaN : {nl::quiet_NaN(), nl::signaling_NaN()})
+ for (unsigned n = 0; n < g_max_n; ++n)
+ assert(std::isnan(std::legendre(n, NaN)));
}
+
+ { // simple sample points for n=0..127 should not produce NaNs.
+ for (Real x : sample_points<Real>())
+ for (unsigned n = 0; n < g_max_n; ++n)
+ assert(!std::isnan(std::legendre(n, x)));
+ }
+
+ { // For x with abs(x) > 1 an domain_error exception should be thrown.
+#ifndef _LIBCPP_NO_EXCEPTIONS
+ for (double absX : {2.0, 7.77, 42.42, std::numeric_limits<double>::infinity()})
+ for (Real x : {-absX, absX})
+ for (unsigned n = 0; n < g_max_n; ++n) {
+ bool throws = false;
+ try {
+ std::legendre(n, x);
+ } catch (const std::domain_error&) {
+ throws = true;
+ }
+ assert(throws);
+ }
#endif // _LIBCPP_NO_EXCEPTIONS
-}
+ }
-template <class T>
-void testLegendreAnalytic(const T x, const T AbsTolerance, const T RelTolerance) {
- assert(!std::isnan(x));
- const auto compareFloatingPoint = [AbsTolerance, RelTolerance](const T Result, const T ExpectedResult) {
- if (std::isinf(ExpectedResult) && std::isinf(Result))
- return true;
+ { // check against analytic polynoms for order n=0..6
+ for (Real x : sample_points<Real>()) {
+ const auto l0 = [](Real) -> Real { return 1; };
+ const auto l1 = [](Real y) -> Real { return y; };
+ const auto l2 = [](Real y) -> Real { return (3 * y * y - 1) / 2; };
+ const auto l3 = [](Real y) -> Real { return (5 * y * y - 3) * y / 2; };
+ const auto l4 = [](Real y) -> Real { return (35 * std::pow(y, 4) - 30 * y * y + 3) / 8; };
+ const auto l5 = [](Real y) -> Real { return (63 * std::pow(y, 4) - 70 * y * y + 15) * y / 8; };
+ const auto l6 = [](Real y) -> Real {
+ return (231 * std::pow(y, 6) - 315 * std::pow(y, 4) + 105 * y * y - 5) / 16;
+ };
+
+ const CompareFloatingValues<Real> compare;
+ assert(compare(std::legendre(0, x), l0(x)));
+ assert(compare(std::legendre(1, x), l1(x)));
+ assert(compare(std::legendre(2, x), l2(x)));
+ assert(compare(std::legendre(3, x), l3(x)));
+ assert(compare(std::legendre(4, x), l4(x)));
+ assert(compare(std::legendre(5, x), l5(x)));
+ assert(compare(std::legendre(6, x), l6(x)));
+ }
+ }
- if (std::isnan(ExpectedResult) || std::isnan(Result))
- return false;
+ { // checks std::legendref for bitwise equality with std::legendre(unsigned, float)
+ if constexpr (std::is_same_v<Real, float>)
+ for (unsigned n = 0; n < g_max_n; ++n)
+ for (float x : sample_points<float>())
+ assert(std::legendre(n, x) == std::legendref(n, x));
+ }
- const T Tolerance = AbsTolerance + std::abs(ExpectedResult) * RelTolerance;
- return std::abs(Result - ExpectedResult) < Tolerance;
- };
-
- const auto l0 = [](T) { return T(1); };
- const auto l1 = [](T y) { return y; };
- const auto l2 = [](T y) { return (T(3) * y * y - T(1)) / T(2); };
- const auto l3 = [](T y) { return (T(5) * y * y - T(3)) * y / T(2); };
- const auto l4 = [](T y) { return (T(35) * y * y * y * y - T(30) * y * y + T(3)) / T(8); };
- const auto l5 = [](T y) { return (T(63) * y * y * y * y - T(70) * y * y + T(15)) * y / T(8); };
- const auto l6 = [](T y) {
- const T y2 = y * y;
- return (T(231) * y2 * y2 * y2 - T(315) * y2 * y2 + T(105) * y2 - T(5)) / T(16);
- };
-
- assert(compareFloatingPoint(std::legendre(0, x), l0(x)));
- assert(compareFloatingPoint(std::legendre(1, x), l1(x)));
- assert(compareFloatingPoint(std::legendre(2, x), l2(x)));
- assert(compareFloatingPoint(std::legendre(3, x), l3(x)));
- assert(compareFloatingPoint(std::legendre(4, x), l4(x)));
- assert(compareFloatingPoint(std::legendre(5, x), l5(x)));
- assert(compareFloatingPoint(std::legendre(6, x), l6(x)));
-}
+ { // checks std::legendrel for bitwise equality with std::legendre(unsigned, long double)
+ if constexpr (std::is_same_v<Real, long double>)
+ for (unsigned n = 0; n < g_max_n; ++n)
+ for (long double x : sample_points<long double>()) {
+ assert(std::legendre(n, x) == std::legendrel(n, x));
+ assert(std::legendre(n, x) == std::legendrel(n, x));
+ }
+ }
-template <class T>
-void testLegendre(const T AbsTolerance, const T RelTolerance) {
- testLegendreNaNPropagation<T>();
- testLegendreThrows<T>(T(-5));
- testLegendreThrows<T>(T(5));
+ { // evaluation at x=1: P_n(1) = 1
+ const CompareFloatingValues<Real> compare;
+ for (unsigned n = 0; n < g_max_n; ++n)
+ assert(compare(std::legendre(n, Real{1}), 1));
+ }
- const T Samples[] = {T(-1.0), T(-0.5), T(-0.1), T(0.0), T(0.1), T(0.5), T(1.0)};
+ { // evaluation at x=-1: P_n(-1) = (-1)^n
+ const CompareFloatingValues<Real> compare;
+ for (unsigned n = 0; n < g_max_n; ++n)
+ assert(compare(std::legendre(n, Real{-1}), std::pow(-1, n)));
+ }
- for (T x : Samples) {
- testLegendreNotNaN(x);
- testLegendreAnalytic(x, AbsTolerance, RelTolerance);
+ { // evaluation at x=0:
+ // P_{2n }(0) = (-1)^n (2n-1)!! / (2n)!!
+ // P_{2n+1}(0) = 0
+ const CompareFloatingValues<Real> compare;
+ Real doubleFactorials{1};
+ for (unsigned n = 0; n < g_max_n; ++n) {
+ if (n & 1) // odd
+ doubleFactorials *= n;
+ else if (n != 0) // even and not zero
+ doubleFactorials /= n;
+
+ assert(compare(std::legendre(n, Real{0}), (n & 1) ? Real{0} : std::pow(-1, n / 2) * doubleFactorials));
+ }
}
}
+struct TestFloat {
+ template <class Real>
+ void operator()() {
+ test<Real>();
+ }
+};
+
+struct TestInt {
+ template <class Integer>
+ void operator()() {
+ // checks that std::legendre(unsigned, Integer) actually wraps std::legendre(unsigned, double)
+ for (unsigned n = 0; n < g_max_n; ++n)
+ for (Integer x : {-1, 0, 1})
+ assert(std::legendre(n, x) == std::legendre(n, static_cast<double>(x)));
+ }
+};
+
int main() {
- testLegendre<float>(1e-6f, 1e-6f);
- testLegendre<double>(1e-9, 1e-9);
- testLegendre<long double>(1e-12, 1e-12);
+ types::for_each(types::floating_point_types(), TestFloat());
+ types::for_each(types::type_list<short, int, long, long long>(), TestInt());
}
More information about the libcxx-commits
mailing list