diff --git a/libcxx/include/__math/special_functions.h b/libcxx/include/__math/special_functions.h index 0b1c753a659ad..8b2b962a85a4d 100644 --- a/libcxx/include/__math/special_functions.h +++ b/libcxx/include/__math/special_functions.h @@ -11,11 +11,13 @@ #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 +#include #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER) # pragma GCC system_header @@ -49,7 +51,7 @@ _LIBCPP_HIDE_FROM_ABI _Real __hermite(unsigned __n, _Real __x) { } if (!__math::isfinite(__H_n)) { - // Overflow occured. Two possible cases: + // Overflow occurred. 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(); @@ -77,6 +79,51 @@ _LIBCPP_HIDE_FROM_ABI double hermite(unsigned __n, _Integer __x) { return std::hermite(__n, static_cast(__x)); } +template +_LIBCPP_HIDE_FROM_ABI _Real __legendre(unsigned __n, _Real __x) { + // The Legendre polynomial P_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) + std::__throw_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; + } + + 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 , int> = 0> +_LIBCPP_HIDE_FROM_ABI double legendre(unsigned __n, _Integer __x) { + return std::legendre(__n, static_cast(__x)); +} + #endif // _LIBCPP_STD_VER >= 17 _LIBCPP_END_NAMESPACE_STD diff --git a/libcxx/include/cmath b/libcxx/include/cmath index 3c22604a683c3..4473b356ca63c 100644 --- a/libcxx/include/cmath +++ b/libcxx/include/cmath @@ -224,6 +224,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 +double legendre(unsigned n, Integer x); // C++17 + floating_point lgamma (arithmetic x); float lgammaf(float x); long double lgammal(long double x); diff --git a/libcxx/modules/std/cmath.inc b/libcxx/modules/std/cmath.inc index fe8ac773c9d1c..8767a4ac0273c 100644 --- a/libcxx/modules/std/cmath.inc +++ b/libcxx/modules/std/cmath.inc @@ -346,12 +346,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..aee7850f0a5e8 --- /dev/null +++ b/libcxx/test/std/numerics/c.math/legendre.pass.cpp @@ -0,0 +1,179 @@ +//===----------------------------------------------------------------------===// +// +// 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 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 +// double legendre(unsigned n, Integer x); + +#include +#include +#include +#include + +#include "type_algorithms.h" + +inline constexpr unsigned g_max_n = 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 +void test() { + { // checks if NaNs are reported correctly (i.e. output == input for input == NaN) + using nl = std::numeric_limits; + 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()) + 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::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 + } + + { // check against analytic polynoms for order n=0..6 + for (Real x : sample_points()) { + 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 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))); + } + } + + { // checks std::legendref for bitwise equality with std::legendre(unsigned, float) + if constexpr (std::is_same_v) + for (unsigned n = 0; n < g_max_n; ++n) + for (float x : sample_points()) + assert(std::legendre(n, x) == std::legendref(n, x)); + } + + { // checks std::legendrel for bitwise equality with std::legendre(unsigned, long double) + if constexpr (std::is_same_v) + for (unsigned n = 0; n < g_max_n; ++n) + for (long double x : sample_points()) { + assert(std::legendre(n, x) == std::legendrel(n, x)); + assert(std::legendre(n, x) == std::legendrel(n, x)); + } + } + + { // evaluation at x=1: P_n(1) = 1 + const CompareFloatingValues compare; + for (unsigned n = 0; n < g_max_n; ++n) + assert(compare(std::legendre(n, Real{1}), 1)); + } + + { // evaluation at x=-1: P_n(-1) = (-1)^n + const CompareFloatingValues compare; + for (unsigned n = 0; n < g_max_n; ++n) + assert(compare(std::legendre(n, 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_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 + void operator()() { + test(); + } +}; + +struct TestInt { + template + 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(x))); + } +}; + +int main() { + types::for_each(types::floating_point_types(), TestFloat()); + types::for_each(types::type_list(), TestInt()); +}