[libcxx-commits] [libcxx] [libc++][math] Mathematical special functions: Implementing `std::assoc_legendre` (PR #92608)
via libcxx-commits
libcxx-commits at lists.llvm.org
Fri May 17 14:47:35 PDT 2024
https://github.com/PaulXiCao created https://github.com/llvm/llvm-project/pull/92608
Implementing the _associated Legendre polynomials_ `std::assoc_legendre`, `std::assoc_legendref`, `std::assoc_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.
The implementation dispatches to `std::legendre` which is waiting for integration in #92175.
>From e4411b2e7cf4cdbe582a363426178992cef3069e Mon Sep 17 00:00:00 2001
From: Paul <{ID}+{username}@users.noreply.github.com>
Date: Fri, 17 May 2024 20:20:28 +0200
Subject: [PATCH 1/2] implementation
---
libcxx/include/CMakeLists.txt | 1 +
libcxx/include/__math/special_functions.h | 105 ++++++++++++++++++++++
libcxx/include/cmath | 9 ++
libcxx/include/module.modulemap | 1 +
libcxx/modules/std/cmath.inc | 2 +
5 files changed, 118 insertions(+)
create mode 100644 libcxx/include/__math/special_functions.h
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..c282b9e33f879
--- /dev/null
+++ b/libcxx/include/__math/special_functions.h
@@ -0,0 +1,105 @@
+// -*- 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/roots.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 __assoc_legendre(unsigned __l, unsigned __m, _Real __x) {
+ // The associated Legendre polynomial P_l^m(x).
+ // The implementation is based on well-known recurrence formulas.
+ // Resources:
+ // - Book: Press - Numerical Recipes - 3rd Edition - p. 294
+ // - Excellect write-up: https://justinwillmert.com/articles/2020/calculating-legendre-polynomials/
+
+ // NOLINTBEGIN(readability-identifier-naming)
+ if (__l < __m)
+ return 0;
+ /*
+ // fixme: uncomment once `std::legendre` is integrated and git rebased.
+ if (__m == 0)
+ return std::legendre(__l, __x);
+ */
+ if (__math::isnan(__x))
+ return __x;
+ else if (__math::fabs(__x) > 1)
+ std::__throw_domain_error("Argument `x` of associated Legendre function is out of range: `|x| <= 1`.");
+
+ _Real __Pmm = 1; // init with P_0^0
+ // Compute P_m^m: Explicit loop unrolling to work around computing square root.
+ // Note: (1-x)*(1+x) is more accurate than (1-x*x)
+ for (unsigned __n = 0; __n < __m / 2; ++__n)
+ __Pmm *= (4 * __n + 1) * (4 * __n + 3) * (1 - __x) * (1 + __x);
+ // Odd m case: Cannot combine two iterations. Thus, needs to compute square root.
+ if (__m & 1)
+ __Pmm *= (2 * __m - 1) * __math::sqrt((1 - __x) * (1 + __x));
+
+ if (__l == __m)
+ return __Pmm;
+
+ // Factoring P^m_m out and multiplying it afterwards. Possible as recursion is linear.
+ _Real __Pml = __x * (2 * __m + 1) /* * __Pmm */; // init with P^m_{m+1}
+ _Real __Pml_prev = 1 /* * __Pmm */; // init with P^m_m
+ for (unsigned __n = __m + 2; __n <= __l; ++__n) {
+ _Real __Pml_prev2 = __Pml_prev; // P^m_{n-2}
+ __Pml_prev = __Pml; // P^m_{n-1}
+ __Pml = (__x * (2 * __n - 1) * __Pml_prev - (__n + __m - 1) * __Pml_prev2) / (__n - __m); // P^m_n
+ }
+ return __Pml * __Pmm;
+ // NOLINTEND(readability-identifier-naming)
+}
+
+inline _LIBCPP_HIDE_FROM_ABI double assoc_legendre(unsigned __l, unsigned __m, double __x) {
+ return std::__assoc_legendre(__l, __m, __x);
+}
+
+inline _LIBCPP_HIDE_FROM_ABI float assoc_legendre(unsigned __l, unsigned __m, float __x) {
+ return static_cast<float>(std::__assoc_legendre(__l, __m, static_cast<double>(__x)));
+}
+
+inline _LIBCPP_HIDE_FROM_ABI long double assoc_legendre(unsigned __l, unsigned __m, long double __x) {
+ return std::__assoc_legendre(__l, __m, __x);
+}
+
+inline _LIBCPP_HIDE_FROM_ABI float assoc_legendref(unsigned __l, unsigned __m, float __x) {
+ return std::assoc_legendre(__l, __m, __x);
+}
+
+inline _LIBCPP_HIDE_FROM_ABI long double assoc_legendrel(unsigned __l, unsigned __m, long double __x) {
+ return std::assoc_legendre(__l, __m, __x);
+}
+
+template <class _Integer, std::enable_if_t<std::is_integral_v<_Integer>, int> = 0>
+_LIBCPP_HIDE_FROM_ABI double assoc_legendre(unsigned __l, unsigned __m, _Integer __x) {
+ return std::assoc_legendre(__l, __m, 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..56c3add7bbac3 100644
--- a/libcxx/include/cmath
+++ b/libcxx/include/cmath
@@ -160,6 +160,14 @@ floating_point asinh (arithmetic x);
float asinhf(float x);
long double asinhl(long double x);
+double assoc_legendre(unsigned l, unsigned m, double x); // C++17
+float assoc_legendre(unsigned l, unsigned m, float x); // C++17
+long double assoc_legendre(unsigned l, unsigned m, long double x); // C++17
+float assoc_legendref(unsigned l, unsigned m, float x); // C++17
+long double assoc_legendrel(unsigned l, unsigned m, long double x); // C++17
+template <class Integer>
+double assoc_legendre(unsigned l, unsigned m, Integer x); // C++17
+
floating_point atanh (arithmetic x);
float atanhf(float x);
long double atanhl(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/include/module.modulemap b/libcxx/include/module.modulemap
index 70dac2f19846b..034796fcb6732 100644
--- a/libcxx/include/module.modulemap
+++ b/libcxx/include/module.modulemap
@@ -1488,6 +1488,7 @@ module std_private_math_modulo [system] { header "__mat
module std_private_math_remainder [system] { header "__math/remainder.h" }
module std_private_math_roots [system] { header "__math/roots.h" }
module std_private_math_rounding_functions [system] { header "__math/rounding_functions.h" }
+module std_private_math_special_functions [system] { header "__math/special_functions.h" }
module std_private_math_traits [system] { header "__math/traits.h" }
module std_private_math_trigonometric_functions [system] { header "__math/trigonometric_functions.h" }
diff --git a/libcxx/modules/std/cmath.inc b/libcxx/modules/std/cmath.inc
index a463c1e3ccf86..6800eebeb747d 100644
--- a/libcxx/modules/std/cmath.inc
+++ b/libcxx/modules/std/cmath.inc
@@ -268,12 +268,14 @@ export namespace std {
using std::assoc_laguerre;
using std::assoc_laguerref;
using std::assoc_laguerrel;
+#endif
// [sf.cmath.assoc.legendre], associated Legendre functions
using std::assoc_legendre;
using std::assoc_legendref;
using std::assoc_legendrel;
+#if 0
// [sf.cmath.beta], beta function
using std::beta;
using std::betaf;
>From 57946ec9515d421506dce2ee56c830f6641f162c Mon Sep 17 00:00:00 2001
From: Paul <{ID}+{username}@users.noreply.github.com>
Date: Fri, 17 May 2024 23:42:02 +0200
Subject: [PATCH 2/2] tests
---
.../numerics/c.math/assoc_legendre.pass.cpp | 221 ++++++++++++++++++
1 file changed, 221 insertions(+)
create mode 100644 libcxx/test/std/numerics/c.math/assoc_legendre.pass.cpp
diff --git a/libcxx/test/std/numerics/c.math/assoc_legendre.pass.cpp b/libcxx/test/std/numerics/c.math/assoc_legendre.pass.cpp
new file mode 100644
index 0000000000000..67dac84b414da
--- /dev/null
+++ b/libcxx/test/std/numerics/c.math/assoc_legendre.pass.cpp
@@ -0,0 +1,221 @@
+//===----------------------------------------------------------------------===//
+//
+// 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 assoc_legendre( unsigned l, unsigned m, double x);
+// float assoc_legendre( unsigned l, unsigned m, float x);
+// long double assoc_legendre( unsigned l, unsigned m, long double x);
+// float assoc_legendref(unsigned l, unsigned m, float x);
+// long double assoc_legendrel(unsigned l, unsigned m, long double x);
+// template <class Integer>
+// double assoc_legendre( unsigned l, unsigned m, Integer x);
+
+#include <array>
+#include <cassert>
+#include <functional>
+#include <cmath>
+#include <limits>
+#include <iostream>
+#include <numeric>
+
+#include "type_algorithms.h"
+
+inline constexpr unsigned g_max_lm = 128;
+
+template <class T>
+std::array<T, 7> sample_points() {
+ return {-1.0, -0.5, -0.1, 0.0, 0.1, 0.5, 1.0};
+}
+
+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;
+ }
+
+ 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 T, std::size_t N>
+std::array<T, N> range() {
+ std::array<T, N> arr;
+ std::iota(arr.begin(), arr.end(), T{0});
+ return arr;
+}
+
+template <class Real>
+void test() {
+#if 1
+ { // checks if NaNs are reported correctly (i.e. output == input for input == NaN)
+ using nl = std::numeric_limits<Real>;
+ for (unsigned l = 0; l < g_max_lm; ++l)
+ for (unsigned m = 0; m < g_max_lm; ++m)
+ for (Real NaN : {nl::quiet_NaN(), nl::signaling_NaN()})
+ assert(std::isnan(std::assoc_legendre(l, m, NaN)));
+ }
+
+ { // simple sample points for l,m=0..127 should not produce NaNs.
+ for (unsigned l = 0; l < g_max_lm; ++l)
+ for (unsigned m = 0; m < g_max_lm; ++m)
+ for (Real x : sample_points<Real>())
+ assert(!std::isnan(std::assoc_legendre(l, m, x)));
+ }
+
+ { // For x with abs(x) > 1 a domain_error exception should be thrown.
+# ifndef _LIBCPP_NO_EXCEPTIONS
+ for (unsigned l = 0; l < g_max_lm; ++l)
+ for (unsigned m = 0; m < g_max_lm; ++m) {
+ const auto is_domain_error = [l, m](Real x) -> bool {
+ try {
+ std::assoc_legendre(l, m, x);
+ } catch (const std::domain_error&) {
+ return true;
+ }
+ return false;
+ };
+
+ Real inf = std::numeric_limits<Real>::infinity();
+ for (Real absX : {std::nextafter(Real(1), inf), Real(2.0), Real(7.77), Real(42.42), inf})
+ for (Real x : {-absX, absX})
+ assert(is_domain_error(x));
+ }
+# endif // _LIBCPP_NO_EXCEPTIONS
+ }
+#endif
+
+#if 1
+ { // check against analytic polynoms for order n=0..6
+ using Func = std::function<Real(Real)>;
+ const Func P00 = [](Real) { return 1; };
+
+ const Func P10 = [](Real y) -> Real { return y; };
+ const Func P11 = [](Real y) -> Real { return std::sqrt((1 - y) * (1 + y)); };
+
+ const Func P20 = [](Real y) -> Real { return (3 * y * y - 1) / 2; };
+ const Func P21 = [](Real y) -> Real { return 3 * y * std::sqrt((1 - y) * (1 + y)); };
+ const Func P22 = [](Real y) -> Real { return 3 * (1 - y) * (1 + y); };
+
+ const Func P30 = [](Real y) -> Real { return (5 * y * y - 3) * y / 2; };
+ const Func P31 = [](Real y) -> Real { return Real(1.5) * (5 * y * y - 1) * std::sqrt((1 - y) * (1 + y)); };
+ const Func P32 = [](Real y) -> Real { return 15 * y * (1 - y) * (1 + y); };
+ const Func P33 = [](Real y) -> Real { return 15 * std::pow((1 - y) * (1 + y), Real(1.5)); };
+
+ const Func P40 = [](Real y) -> Real { return (35 * std::pow(y, 4) - 30 * y * y + 3) / 8; };
+ const Func P41 = [](Real y) -> Real { return Real(2.5) * y * (7 * y * y - 3) * std::sqrt((1 - y) * (1 + y)); };
+ const Func P42 = [](Real y) -> Real { return Real(7.5) * (7 * y * y - 1) * (1 - y) * (1 + y); };
+ const Func P43 = [](Real y) -> Real { return 105 * y * std::pow((1 - y) * (1 + y), Real(1.5)); };
+ const Func P44 = [](Real y) -> Real { return 105 * std::pow((1 - y) * (1 + y), Real(2)); };
+
+ unsigned l = 0;
+ unsigned m = 0;
+ for (auto P : {P00, P10, P11, P20, P21, P22, P30, P31, P32, P33, P40, P41, P42, P43, P44}) {
+ for (Real x : sample_points<Real>()) {
+ const CompareFloatingValues<Real> compare;
+ assert(compare(std::assoc_legendre(l, m, x), P(x)));
+ }
+
+ if (l == m) {
+ ++l;
+ m = 0;
+ } else
+ ++m;
+ }
+ }
+#endif
+
+#if 1
+ { // checks std::assoc_legendref for bitwise equality with std::assoc_legendre(unsigned, float)
+ if constexpr (std::is_same_v<Real, float>)
+ for (unsigned l = 0; l < g_max_lm; ++l)
+ for (unsigned m = 0; m < g_max_lm; ++m)
+ for (float x : sample_points<float>())
+ assert(std::assoc_legendre(l, m, x) == std::assoc_legendref(l, m, x));
+ }
+
+ { // checks std::assoc_legendrel for bitwise equality with std::assoc_legendre(unsigned, long double)
+ if constexpr (std::is_same_v<Real, long double>)
+ for (unsigned l = 0; l < g_max_lm; ++l)
+ for (unsigned m = 0; m < g_max_lm; ++m)
+ for (long double x : sample_points<long double>())
+ assert(std::assoc_legendre(l, m, x) == std::assoc_legendrel(l, m, x));
+ }
+#endif
+
+#if 0
+ { // evaluation at x=1: P_n(1) = 1
+ const CompareFloatingValues<Real> compare;
+ for (unsigned n = 0; n < g_max_lm; ++n)
+ assert(compare(std::assoc_legendre(l, m, Real{1}), 1));
+ }
+
+ { // evaluation at x=-1: P_n(-1) = (-1)^n
+ const CompareFloatingValues<Real> compare;
+ for (unsigned n = 0; n < g_max_lm; ++n)
+ assert(compare(std::assoc_legendre(l, m, Real{-1}), std::pow(-1, n)));
+ }
+
+ { // 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_lm; ++n) {
+ if (n & 1) // odd
+ doubleFactorials *= n;
+ else if (n != 0) // even and not zero
+ doubleFactorials /= n;
+
+ assert(compare(std::assoc_legendre(l, m, Real{0}), (n & 1) ? Real{0} : std::pow(-1, n / 2) * doubleFactorials));
+ }
+ }
+#endif
+}
+
+struct TestFloat {
+ template <class Real>
+ void operator()() {
+ test<Real>();
+ }
+};
+
+struct TestInt {
+ template <class Integer>
+ void operator()() {
+ // checks that std::assoc_legendre(unsigned, unsigned, Integer) actually wraps
+ // std::assoc_legendre(unsigned, unsigned, double)
+ for (unsigned l = 0; l < g_max_lm; ++l)
+ for (unsigned m = 0; m < g_max_lm; ++m)
+ for (Integer x : {-1, 0, 1})
+ assert(std::assoc_legendre(l, m, x) == std::assoc_legendre(l, m, static_cast<double>(x)));
+ }
+};
+
+int main() {
+ 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