From 9c6716e910664ae0c3348f400a7728ef8d2109ff 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/__math/special_functions.h | 49 +++++++- libcxx/include/cmath | 8 ++ libcxx/modules/std/cmath.inc | 2 + .../std/numerics/c.math/legendre.pass.cpp | 111 ++++++++++++++++++ 4 files changed, 169 insertions(+), 1 deletion(-) create mode 100644 libcxx/test/std/numerics/c.math/legendre.pass.cpp 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..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 + +// + +// 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 + +template +void testLegendreNaNPropagation() { + const unsigned MaxN = 127; + const T x = std::numeric_limits::quiet_NaN(); + for (unsigned n = 0; n <= MaxN; ++n) { + assert(std::isnan(std::legendre(n, x))); + } +} + +template +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 +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 +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 +void testLegendre(const T AbsTolerance, const T RelTolerance) { + testLegendreNaNPropagation(); + testLegendreThrows(T(-5)); + testLegendreThrows(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(1e-6f, 1e-6f); + testLegendre(1e-9, 1e-9); + testLegendre(1e-12, 1e-12); +} From 06b211375d094dfbd798178f10a9f0afc09bffa0 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 --- .../std/numerics/c.math/legendre.pass.cpp | 206 ++++++++++++------ 1 file changed, 137 insertions(+), 69 deletions(-) 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 // double legendre(unsigned n, Integer x); +#include #include #include #include +#include "type_algorithms.h" + +inline constexpr unsigned g_max_n = 128; + template -void testLegendreNaNPropagation() { - const unsigned MaxN = 127; - const T x = std::numeric_limits::quiet_NaN(); - for (unsigned n = 0; n <= MaxN; ++n) { - assert(std::isnan(std::legendre(n, x))); - } +std::array sample_points() { + return {-1.0, -0.5, -0.1, 0.0, 0.1, 0.5, 1.0}; } -template -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 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; } -} -template -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 +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 -} + } -template -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()) { + 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))); + } + } - 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) + for (unsigned n = 0; n < g_max_n; ++n) + for (float x : sample_points()) + 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) + 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)); + } + } -template -void testLegendre(const T AbsTolerance, const T RelTolerance) { - testLegendreNaNPropagation(); - testLegendreThrows(T(-5)); - testLegendreThrows(T(5)); + { // 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)); + } - 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 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 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() { - testLegendre(1e-6f, 1e-6f); - testLegendre(1e-9, 1e-9); - testLegendre(1e-12, 1e-12); + types::for_each(types::floating_point_types(), TestFloat()); + types::for_each(types::type_list(), TestInt()); }