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 +#include + +#if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER) +# pragma GCC system_header +#endif + +_LIBCPP_BEGIN_NAMESPACE_STD + +#if _LIBCPP_STD_VER >= 17 + +template +_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(std::__assoc_legendre(__l, __m, static_cast(__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 , int> = 0> +_LIBCPP_HIDE_FROM_ABI double assoc_legendre(unsigned __l, unsigned __m, _Integer __x) { + return std::assoc_legendre(__l, __m, static_cast(__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 +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 #include +#include <__math/special_functions.h> #include #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 + +// + +// 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 +// double assoc_legendre( unsigned l, unsigned m, Integer x); + +#include +#include +#include +#include +#include +#include +#include + +#include "type_algorithms.h" + +inline constexpr unsigned g_max_lm = 128; + +template +std::array sample_points() { + return {-1.0, -0.5, -0.1, 0.0, 0.1, 0.5, 1.0}; +} + +template +class CompareFloatingValues { +private: + Real tol; + +public: + CompareFloatingValues() { + if (std::is_same_v) + tol = 1e-6f; + else if (std::is_same_v) + 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 +std::array range() { + std::array arr; + std::iota(arr.begin(), arr.end(), T{0}); + return arr; +} + +template +void test() { +#if 1 + { // checks if NaNs are reported correctly (i.e. output == input for input == NaN) + using nl = std::numeric_limits; + 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()) + 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::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; + 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()) { + const CompareFloatingValues 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) + for (unsigned l = 0; l < g_max_lm; ++l) + for (unsigned m = 0; m < g_max_lm; ++m) + for (float x : sample_points()) + 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) + for (unsigned l = 0; l < g_max_lm; ++l) + for (unsigned m = 0; m < g_max_lm; ++m) + for (long double x : sample_points()) + 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 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 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 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 + void operator()() { + test(); + } +}; + +struct TestInt { + template + 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(x))); + } +}; + +int main() { + types::for_each(types::floating_point_types(), TestFloat()); + types::for_each(types::type_list(), TestInt()); +}