diff --git a/libcxx/include/__cmath b/libcxx/include/__cmath new file mode 100644 --- /dev/null +++ b/libcxx/include/__cmath @@ -0,0 +1,50 @@ +// -*- C++ -*- +//===------------------------------ __cmath -------------------------------===// +// +// 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__CMATH +#define _LIBCPP__CMATH + +#include <__config> + +#if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER) +#pragma GCC system_header +#endif + +_LIBCPP_BEGIN_NAMESPACE_STD + +// __libcpp_assoc_legendre + +_LIBCPP_FUNC_VIS +float __libcpp_assoc_legendref(unsigned int __lcpp_l, unsigned int __lcpp_m, + float __lcpp_x); + +_LIBCPP_FUNC_VIS +double __libcpp_assoc_legendre(unsigned int __lcpp_l, unsigned int __lcpp_m, + double __lcpp_x); + +_LIBCPP_FUNC_VIS +long double __libcpp_assoc_legendrel(unsigned int __lcpp_l, + unsigned int __lcpp_m, + long double __lcpp_x); + + +// __libcpp_legendre + +_LIBCPP_FUNC_VIS +float __libcpp_legendref(unsigned int __lcpp_n, float __lcpp_x); + +_LIBCPP_FUNC_VIS +double __libcpp_legendre(unsigned int __lcpp_n, double __lcpp_x); + +_LIBCPP_FUNC_VIS +long double __libcpp_legendrel(unsigned int __lcpp_n, long double __lcpp_x); + +_LIBCPP_END_NAMESPACE_STD + +#endif // _LIBCPP__CMATH diff --git a/libcxx/include/cmath b/libcxx/include/cmath --- a/libcxx/include/cmath +++ b/libcxx/include/cmath @@ -296,6 +296,16 @@ float truncf(float x); long double truncl(long double x); +// C++17 Mathematical special functions [sf.cmath] + +floating_point assoc_legendre (unsigned l, unsigned m, arithmetic x); +float assoc_legendref(unsigned l, unsigned m, float x); +long double assoc_legendrel(unsigned l, unsigned m, long double x); + +floating_point legendre (unsigned n, arithmetic x); +float legendref(unsigned n, float x); +long double legendrel(unsigned n, long double x); + } // std */ @@ -303,6 +313,7 @@ #include <__config> #include #include +#include <__cmath> #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER) #pragma GCC system_header @@ -523,6 +534,87 @@ using ::truncl; #if _LIBCPP_STD_VER > 14 + +// assoc_legendre + +inline _LIBCPP_INLINE_VISIBILITY +float +assoc_legendref(unsigned int __lcpp_l, unsigned int __lcpp_m, float __lcpp_x) { + return _VSTD::__libcpp_assoc_legendref(__lcpp_l, __lcpp_m, __lcpp_x); +} + +inline _LIBCPP_INLINE_VISIBILITY +double +assoc_legendre(unsigned int __lcpp_l, unsigned int __lcpp_m, double __lcpp_x) { + return _VSTD::__libcpp_assoc_legendre(__lcpp_l, __lcpp_m, __lcpp_x); +} + +inline _LIBCPP_INLINE_VISIBILITY +long double +assoc_legendrel(unsigned int __lcpp_l, unsigned int __lcpp_m, long double __lcpp_x) { + return _VSTD::__libcpp_assoc_legendrel(__lcpp_l, __lcpp_m, __lcpp_x); +} + +inline _LIBCPP_INLINE_VISIBILITY +float +assoc_legendre(unsigned int __lcpp_l, unsigned int __lcpp_m, float __lcpp_x) { + return _VSTD::assoc_legendref(__lcpp_l, __lcpp_m, __lcpp_x); +} + +inline _LIBCPP_INLINE_VISIBILITY +long double +assoc_legendre(unsigned int __lcpp_l, unsigned int __lcpp_m, long double __lcpp_x) { + return _VSTD::assoc_legendrel(__lcpp_l, __lcpp_m, __lcpp_x); +} + +template +inline _LIBCPP_INLINE_VISIBILITY +typename std::enable_if::value, double>::type +assoc_legendre(unsigned int __lcpp_l, unsigned int __lcpp_m, _A1 __lcpp_x) { + return _VSTD::assoc_legendre(__lcpp_l, __lcpp_m, (double)__lcpp_x); +} + +// legendre + +inline _LIBCPP_INLINE_VISIBILITY +float +legendref(unsigned int __lcpp_n, float __lcpp_x) { + return _VSTD::__libcpp_legendref(__lcpp_n, __lcpp_x); +} + +inline _LIBCPP_INLINE_VISIBILITY +double +legendre(unsigned int __lcpp_n, double __lcpp_x) { + return _VSTD::__libcpp_legendre(__lcpp_n, __lcpp_x); +} + +inline _LIBCPP_INLINE_VISIBILITY +long double +legendrel(unsigned int __lcpp_n, long double __lcpp_x) { + return _VSTD::__libcpp_legendrel(__lcpp_n, __lcpp_x); +} + +inline _LIBCPP_INLINE_VISIBILITY +float +legendre(unsigned int __lcpp_n, float __lcpp_x) { + return _VSTD::legendref(__lcpp_n, __lcpp_x); +} + +inline _LIBCPP_INLINE_VISIBILITY +long double +legendre(unsigned int __lcpp_n, long double __lcpp_x) { + return _VSTD::legendrel(__lcpp_n, __lcpp_x); +} + +template +inline _LIBCPP_INLINE_VISIBILITY +typename std::enable_if::value, double>::type +legendre(unsigned int __lcpp_n, _A1 __lcpp_x) { + return _VSTD::legendre(__lcpp_n, (double)__lcpp_x); +} + +// hypot + inline _LIBCPP_INLINE_VISIBILITY float hypot( float x, float y, float z ) { return sqrt(x*x + y*y + z*z); } inline _LIBCPP_INLINE_VISIBILITY double hypot( double x, double y, double z ) { return sqrt(x*x + y*y + z*z); } inline _LIBCPP_INLINE_VISIBILITY long double hypot( long double x, long double y, long double z ) { return sqrt(x*x + y*y + z*z); } diff --git a/libcxx/lib/CMakeLists.txt b/libcxx/lib/CMakeLists.txt --- a/libcxx/lib/CMakeLists.txt +++ b/libcxx/lib/CMakeLists.txt @@ -23,6 +23,8 @@ endif() endif() +list(APPEND LIBCXX_SOURCES ../src/cmath/legendre.cpp) + # Add all the headers to the project for IDEs. if (LIBCXX_CONFIGURE_IDE) file(GLOB_RECURSE LIBCXX_HEADERS ${CMAKE_CURRENT_SOURCE_DIR}/../include/*) diff --git a/libcxx/src/cmath/legendre.cpp b/libcxx/src/cmath/legendre.cpp new file mode 100644 --- /dev/null +++ b/libcxx/src/cmath/legendre.cpp @@ -0,0 +1,209 @@ +//===----------------------- legendre.cpp -----------------------*- 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 +// +//===----------------------------------------------------------------------===// +/// +/// \file +/// This file contains the internal implementations of std::legendre and +/// and std::assoc_legendre (as of C++17) +/// std::assoc_legendre, std::assoc_legendref, and std::assoc_legendrel +/// correspond to __libcpp_assoc_legendre, __libcpp_assoc_legendref, and +/// __libcpp_assoc_legendrel. +/// std::legendre, std::legendref, and std::legendrel correspond to +/// __libcpp_legendre, __libcpp_legendref, and __libcpp_legendrel. +/// +//===----------------------------------------------------------------------===// + +#include +#include +#include + +_LIBCPP_BEGIN_NAMESPACE_STD + +/// \brief Evaluates the Legendre polynomial \f$ P_{n}(x) \f$. +/// +/// The implementation is based on the recurrence formula +/// \f[ +/// (n+1)P_{n+1}(x) = (2n+1)xP_{n}(x) - nP_{n-1}(x) +/// \f] +/// where \f$ P_{0}(x) = 1 \f$ and \f$ P_{1}(x) = x \f$. +/// Press, William H., et al. Numerical recipes 3rd edition: The art of +/// scientific computing. Cambridge university press, 2007, p. 182. +/// \tparam T one of float, double, or long double +/// \param [in] __n corresponds to the degree \f$ n \f$. +/// \param [in] __x corresponds to the argument \f$ x \f$. +/// \return the Legendre polynomial \f$ P_{n}(x) \f$ +template +inline _LIBCPP_INLINE_VISIBILITY _Real +__libcpp_legendre_recurrence(unsigned __n, _Real __x) _NOEXCEPT { + if (__n == 0u) + return _Real(1); + + _Real __t2(1); + _Real __t1 = __x; + for (unsigned __i = 1; __i < __n; ++__i) { + const _Real __k = _Real(__i); + _Real __t0 = ((_Real(2) * __k + _Real(1)) * __x * __t1 - __k * __t2) / + (__k + _Real(1)); + __t2 = __t1; + __t1 = __t0; + } + return __t1; +} + +/// \brief Evaluates the Legendre polynomial \f$ P_{n}(x) \f$. +/// +/// \tparam T one of float, double, or long double +/// \param [in] __n corresponds to the degree \f$ n \f$. +/// \param [in] __x corresponds to the argument \f$ x \f$. +/// \return the Legendre polynomial \f$ P_{n}(x) \f$ +/// \throw std::domain_error if __x is not inside the interval [-1,1]. +template +inline _LIBCPP_INLINE_VISIBILITY _Real __libcpp_legendre_impl(unsigned __n, + _Real __x) { + if (std::isnan(__x)) + return std::numeric_limits<_Real>::quiet_NaN(); + + if (std::abs(__x) > _Real(1)) + __throw_domain_error("Argument of legendre function is out of range"); + + return __libcpp_legendre_recurrence(__n, __x); +} + +/// \brief Evaluates the associated Legendre polynomial \f$ P_{l}^{m} \f$. +/// +/// The implementation is based on the recurrence formula +/// \f[ +/// (l-m+1)P_{l+1}^{m}(x) = (2l+1)xP_{l}^{m}(x) - +/// (l+m)P_{l-1}^{m}(x) +/// \f] +/// with +/// \f[ +/// P_{m}^{m}(x) = \sqrt{1 - x^2}^{m}\frac{(2m)!}{2^m m!} +/// \f] +/// and +/// \f[ +/// P_{m-1}^{m}(x) = 0 +/// \f] +/// +/// \note The starting point of the recursion grows exponentially with +/// \f$ m \f$. For large \f$ m \f$, we have the following relation: +/// \f[ +/// P_{m}^{m}(x) \approx \sqrt{1 - +/// x^2}^{m}\sqrt{2} 2^{n} \exp( n(\ln n - 1 )) +/// \f] +/// For example, for \f$ m = 40\f$, we already have +/// \f$ P_{40}^{40}(0) \approx 8 \cdot 10^{58} \f$. +/// \attention The so-called Condon-Shortley phase term is omitted in the C++17 +/// standard's definition of std::assoc_legendre. +/// \param [in] __l corresponds to \f$ l \f$. +/// \param [in] __m corresponds to \f$ m \f$. +/// \param [in] __x corresponds to the argument \f$ x \f$. +/// \return \f$ P_{l}^{m}(x) \f$ +template +inline _LIBCPP_INLINE_VISIBILITY _Real __libcpp_assoc_legendre_recurrence( + unsigned __l, unsigned __m, _Real __x) _NOEXCEPT { + if (__m == 0u) + return __libcpp_legendre_recurrence(__l, __x); + + if (__l < __m) + return _Real(0); + + if (__l == 0u) + return _Real(1); + + _Real __pmm = _Real(1); + // Note: (1-x)*(1+x) is more accurate than (1-x*x) + // "What Every Computer Scientist Should Know About Floating-Point + // Arithmetic", David Goldberg, p. 38 + const _Real __t = std::sqrt((_Real(1) - __x) * (_Real(1) + __x)) / (_Real(2)); + for (unsigned __i = 2u * __m; __i > __m; --__i) + __pmm *= __t * __i; + + if (__l == __m) + return __pmm; + + // Actually, we'd start with _pmm but it grows exponentially with __m. + // Luckily, the recursion scales. So we can start with 1 and multiply + // afterwards. + _Real __t2 = _Real(1); + _Real __t1 = _Real(2u * __m + 1u) * __x; // first iteration unfolded + for (unsigned __i = __m + 1u; __i < __l; ++__i) { + // As soon as one of the terms becomes inf, this will quickly lead to NaNs. + // float just doesn't do it for the whole range up to l==127. + const _Real __t0 = + (_Real(2u * __i + 1u) * __x * __t1 - _Real(__i + __m) * __t2) / + _Real(__i - __m + 1u); + __t2 = __t1; + __t1 = __t0; + } + return __t1 * __pmm; +} + +/// \brief Evaluates the associated Legendre polynomial \f$ P_{l}^{m} \f$. +/// +/// \param [in] __l corresponds to \f$ l \f$. +/// \param [in] __m corresponds to \f$ m \f$. +/// \param [in] __x corresponds to the argument \f$ x \f$. +/// \return \f$ P_{l}^{m}(x) \f$ +/// \throw std::domain_error if __x is not inside the interval [-1,1]. +template +inline _LIBCPP_INLINE_VISIBILITY _Real +__libcpp_assoc_legendre_impl(unsigned __l, unsigned __m, _Real __x) { + if (std::isnan(__x)) + return std::numeric_limits<_Real>::quiet_NaN(); + + if (std::abs(__x) > _Real(1)) + __throw_domain_error("Argument of assoc_legendre function is out of range"); + + return __libcpp_assoc_legendre_recurrence(__l, __m, __x); +} + +/// \brief Evaluates the associated Legendre polynomial \f$ P_{l}^{m} \f$. +/// +/// We use the double version internally because assoc_legendre +/// is prone to overflow. It can be shown that assoc_legendre is bound +/// by a value less than std::numeric_limits::max() for the whole +/// range required by the standard, __lcpp_x in [-1,1}, __lcpp_l < 128, +/// and __lcpp_m < 128. +/// This is elaborated in more detail in the unit tests +/// See Eq. (22) on p. 183 of Lohoefer, G. (1998). +/// Inequalities for the associated Legendre functions. +float __libcpp_assoc_legendref(unsigned int __lcpp_l, unsigned int __lcpp_m, + float __lcpp_x) { + return static_cast(__libcpp_assoc_legendre_impl(__lcpp_l, __lcpp_m, static_cast(__lcpp_x))); +} + +double __libcpp_assoc_legendre(unsigned int __lcpp_l, unsigned int __lcpp_m, + double __lcpp_x) { + return __libcpp_assoc_legendre_impl(__lcpp_l, __lcpp_m, __lcpp_x); +} + +long double __libcpp_assoc_legendrel(unsigned int __lcpp_l, + unsigned int __lcpp_m, + long double __lcpp_x) { + return __libcpp_assoc_legendre_impl(__lcpp_l, __lcpp_m, + __lcpp_x); +} + +/// \brief Evaluates the Legendre polynomial \f$ P_{n}(x) \f$. +/// +/// In contrast to assoc_legendre, legendre won't overflow in common scenarios. +/// Nevertheless, we also use the double implementation to make sure that +/// legendre(n,x) returns assoc_legendre(n,0,x). +float __libcpp_legendref(unsigned int __lcpp_n, float __lcpp_x) { + return static_cast(__libcpp_legendre_impl(__lcpp_n, static_cast(__lcpp_x))); +} + +double __libcpp_legendre(unsigned int __lcpp_n, double __lcpp_x) { + return __libcpp_legendre_impl(__lcpp_n, __lcpp_x); +} + +long double __libcpp_legendrel(unsigned int __lcpp_n, long double __lcpp_x) { + return __libcpp_legendre_impl(__lcpp_n, __lcpp_x); +} + +_LIBCPP_END_NAMESPACE_STD diff --git a/libcxx/test/std/numerics/legendre/assoc_legendre_basic.pass.cpp b/libcxx/test/std/numerics/legendre/assoc_legendre_basic.pass.cpp new file mode 100644 --- /dev/null +++ b/libcxx/test/std/numerics/legendre/assoc_legendre_basic.pass.cpp @@ -0,0 +1,127 @@ +//===----------------------------------------------------------------------===// +// +// 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++98, c++03, c++11, c++14 + +// + +// floating_point assoc_legendre (unsigned l, unsigned m, arithmetic x); +// float assoc_legendref(unsigned l, unsigned m, float x); +// long double assoc_legendrel(unsigned l, unsigned m, long double x); + +#include +#include +#include + +/// \brief std::assoc_legendre(n,m,x) must silently return NaN if x is NaN. +/// +/// \tparam T one of float, double, or long double +template +void testNaNPropagation() { + for (unsigned l = 0; l < 128; ++l) { + for (unsigned m = 0; m < 128; ++m) { + const T ExpectNaN = + std::assoc_legendre(l, m, std::numeric_limits::quiet_NaN()); + assert(std::isnan(ExpectNaN)); + } + } +} + +/// \brief Tests if the associated Legendre polynomial behaves +/// correctly under the transformation x -> -x +/// +/// The expected property is +/// \f$ P_{l} ^{m} (-x) = (-1)^{l + m} P_{l}^{m}(x) \f$. +/// For odd \f$ l + m \f$, this test checks for a root at \f$ x = 0\f$. +/// \tparam T one of float, double, or long double +/// \param [in] l first argument of assoc_legendre +/// \param [in] m second argument of assoc_legendre +/// \param [in] x argument of the Legendre polynomial +/// \param [in] RelTolerance relative tolerance for the value of legendre +/// \param [in] AbsTolerance absolute tolerance for the value of legendre +template +void testParity(unsigned l, unsigned m, T x, T RelTolerance, T AbsTolerance) { + const T Parity = (l + m) % 2 == 0 ? T(1) : T(-1); + const T ExpectedValue = Parity * std::assoc_legendre(l, m, x); + const T Tolerance = RelTolerance * std::abs(ExpectedValue) + AbsTolerance; + const T Error = std::abs(std::assoc_legendre(l, m, -x) - ExpectedValue); + assert(Error <= Tolerance); +} + +template +void testParity(T RelTolerance, T AbsTolerance) { + const T Samples[] = {T(0), T(0.1), T(0.2), T(0.5), T(0.9), T(1.0)}; + for (unsigned l = 0; l < 128; ++l) { + for (unsigned m = 0; m < 128; ++m) { + for (T x : Samples) + testParity(l, m, x, RelTolerance, AbsTolerance); + } + } +} + +/// \brief Tests if the associated Legendre polynomial adheres +/// to an upper bound from the literature. +/// +/// For $ 1 \leq m \leq l $, we have: +/// \f[ +/// \frac{1}{\sqrt{2.22(m+1)}} < +/// \sqrt{\frac{(l-m)!}{(l+m)!}} \cdot +/// \max\limits_{x\in\left[-1,1\right]}\left| P_{l}^{m} (x) \right| +/// < \sqrt[4]{\frac{2^5}{\pi^3 m}} +/// \f] +/// See Eq. (22) on p. 183 of Lohoefer, G. (1998). +/// Inequalities for the associated Legendre functions. +/// Journal of Approximation Theory, 95(2), 178-193. +/// \note For \f$ l \leq 127, the maximum upper bound is +/// 3.44177e+250, which is well below std::numeric_limits::max() +/// of 1.79769e+308 but larger than the maximum value of float (3.40282e+38). +/// The lower bound for the maximum is 6.80027e+249. So the float version +/// will have to deal with overflows. +/// \tparam T one of float, double, or long double +/// \param [in] l first argument of assoc_legendre +/// \param [in] m second argument of assoc_legendre +/// \param [in] x argument of the Legendre polynomial +/// \param [in] RelTolerance relative tolerance for the value of legendre +/// \param [in] AbsTolerance absolute tolerance for the value of legendre +template +void testUpperBound(unsigned l, unsigned m, T x, T RelTolerance, + T AbsTolerance) { + const T Pi = T(3.14159265358979323846264338327950288L); + const T LogUpperBound = + std::lgamma(T(l + m + 1u)) / T(2) - std::lgamma(T(l - m + 1u)) / T(2) - + std::log(m) / T(4) + T(1.25) * std::log(2) - T(0.75) * std::log(Pi); + const T UpperBound = std::exp(LogUpperBound); + if (std::isinf(UpperBound)) + return; + + const T UpperBoundWithTolerance = + UpperBound + UpperBound * RelTolerance + AbsTolerance; + assert(std::assoc_legendre(l, m, x) <= UpperBoundWithTolerance); +} + +template +void testUpperBound(T RelTolerance, T AbsTolerance) { + const T Samples[] = {T(0), T(0.1), T(0.2), T(0.5), T(0.9), T(1.0)}; + for (unsigned l = 1; l < 128; ++l) { + for (unsigned m = 1; m <= l; ++m) { + for (T x : Samples) + testUpperBound(l, m, x, RelTolerance, AbsTolerance); + } + } +} + +int main(int, char**) { + testNaNPropagation(); + testNaNPropagation(); + testNaNPropagation(); + + // We only test double here because the float version of std::assoc_legendre + // often returns inf for large l and m. + testParity(1e-12, 0.0); + testUpperBound(1e-12, 0.0); +} diff --git a/libcxx/test/std/numerics/legendre/assoc_legendre_comparisons.pass.cpp b/libcxx/test/std/numerics/legendre/assoc_legendre_comparisons.pass.cpp new file mode 100644 --- /dev/null +++ b/libcxx/test/std/numerics/legendre/assoc_legendre_comparisons.pass.cpp @@ -0,0 +1,413 @@ +//===----------------------------------------------------------------------===// +// +// 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++98, c++03, c++11, c++14 + +// + +// floating_point assoc_legendre (unsigned l, unsigned m, arithmetic x); +// float assoc_legendref(unsigned l, unsigned m, float x); +// long double assoc_legendrel(unsigned l, unsigned m, long double x); + +#include +#include +#include +#include + +template +Real evaluatePolynomial(Iterator CoefficientsBegin, Iterator CoefficientsEnd, + Real x) { + Real Result = Real(0); + for (Iterator Iter = CoefficientsBegin; Iter != CoefficientsEnd; ++Iter) { + Result = Result * x + *Iter; + } + return Result; +} + +template +Real evaluatePolynomial(const Container& Coefficients, Real x) { + return evaluatePolynomial(std::begin(Coefficients), std::end(Coefficients), + x); +} + +template +void compareWithAnalyticExpressions(const T x, TestFunction&& ValueSink) { + std::vector > > TableOfCoefficients; + + // m = 1 and l = 1..4 + TableOfCoefficients.emplace_back(); + TableOfCoefficients.back().push_back({T(1)}); + TableOfCoefficients.back().push_back({T(3), T(0)}); + TableOfCoefficients.back().push_back({T(15) / T(2), T(0), -T(3) / T(2)}); + TableOfCoefficients.back().push_back( + {T(35) / T(2), T(0), -T(15) / T(2), T(0)}); + + // m = 2 and l = 1..4 without the factor (1-x^2)^{m/2} + TableOfCoefficients.emplace_back(); + TableOfCoefficients.back().emplace_back(); + TableOfCoefficients.back().push_back({T(3)}); + TableOfCoefficients.back().push_back({T(15), T(0)}); + TableOfCoefficients.back().push_back({T(105) / T(2), T(0), -T(15) / T(2)}); + + // m = 3 and l = 1..4 without the factor (1-x^2)^{m/2} + TableOfCoefficients.emplace_back(); + TableOfCoefficients.back().emplace_back(); + TableOfCoefficients.back().emplace_back(); + TableOfCoefficients.back().push_back({T(15)}); + TableOfCoefficients.back().push_back({T(105), T(0)}); + + // m = 4 and l = 1..4 without the factor (1-x^2)^{m/2} + TableOfCoefficients.emplace_back(); + TableOfCoefficients.back().emplace_back(); + TableOfCoefficients.back().emplace_back(); + TableOfCoefficients.back().emplace_back(); + TableOfCoefficients.back().push_back({T(105)}); + + for (size_t i = 0; i < TableOfCoefficients.size(); ++i) { + const unsigned m = static_cast(i) + 1; + const T Factor = std::pow((T(1) - x) * (T(1) + x), T(m) / T(2)); + for (size_t j = 0; j < TableOfCoefficients[i].size(); ++j) { + const unsigned int l = static_cast(j) + 1; + const T ExpectedValue = + Factor * evaluatePolynomial(TableOfCoefficients[i][j], x); + ValueSink(l, m, x, ExpectedValue); + } + } +} + +template +void compareWithAnalyticExpressions(const T AbsTolerance, + const T RelTolerance) { + const auto Samples = {T(0), T(0.1), T(0.2), T(0.5), T(0.9), T(1.0)}; + auto testFunction = [RelTolerance, AbsTolerance](unsigned l, unsigned m, T x, + T ExpectedValue) { + const T Result = std::assoc_legendre(l, m, x); + const T Tolerance = std::abs(ExpectedValue) * RelTolerance + AbsTolerance; + const T Error = std::abs(Result - ExpectedValue); + assert(Error <= Tolerance); + }; + for (T x : Samples) + compareWithAnalyticExpressions(x, testFunction); +} + +template +void compareWithTable( + std::map, std::map > + Table, + T RelTolerance, T AbsTolerance) { + for (const auto& Entry : Table) { + const unsigned l = Entry.first.first; + const unsigned m = Entry.first.second; + for (const auto& xy : Entry.second) { + const T ExpectedValue = xy.second; + const T Tolerance = std::abs(ExpectedValue) * RelTolerance + AbsTolerance; + const T Error = + std::abs(std::assoc_legendre(l, m, xy.first) - ExpectedValue); + assert(Error <= Tolerance); + } + } +} + +template +void compareWithMSVC(T RelTolerance, T AbsTolerance) { + // Generated with msvc 191627027 + std::map, std::map > + Table; + Table[{10, 10}] = { + {0.0L, 6.54729075000000000e+08L}, {0.1L, 6.22640835705230832e+08L}, + {0.2L, 5.33848212079902649e+08L}, {0.3L, 4.08571989131588876e+08L}, + {0.4L, 2.73815518201505244e+08L}, {0.5L, 1.55370278540039063e+08L}, + {0.6L, 7.03009991216332912e+07L}, {0.7L, 2.25898063438872844e+07L}, + {0.8L, 3.95889634812671319e+06L}, {0.9L, 1.62117400787842256e+05L}}; + Table[{20, 10}] = { + {0.0L, -1.61205295667431641e+12L}, {0.1L, 3.56041756390717590e+11L}, + {0.2L, 1.46660926765988135e+12L}, {0.3L, -1.09242045429771387e+12L}, + {0.4L, -8.99973487310552368e+11L}, {0.5L, 1.74450229594620776e+12L}, + {0.6L, -4.24471915195055054e+11L}, {0.7L, -1.54836452494957959e+12L}, + {0.8L, 2.39275793328105029e+12L}, {0.9L, 9.90017476694990845e+11L}}; + Table[{20, 20}] = { + {0.0L, 3.19830986772877752e+23L}, {0.1L, 2.89249411469768665e+23L}, + {0.2L, 2.12634078007975140e+23L}, {0.3L, 1.24547341322977873e+23L}, + {0.4L, 5.59388325840124914e+22L}, {0.5L, 1.80108069781796049e+22L}, + {0.6L, 3.68740022490078904e+21L}, {0.7L, 3.80734558806206054e+20L}, + {0.8L, 1.16935276168332288e+19L}, {0.9L, 1.96090497120238360e+16L}}; + Table[{50, 10}] = { + {0.0L, 1.14522382492463460e+16L}, {0.1L, 2.79436816274998850e+15L}, + {0.2L, -9.92643944667355200e+15L}, {0.3L, -9.45138403285157200e+15L}, + {0.4L, 8.89201701866905625e+14L}, {0.5L, 9.18003538846575400e+15L}, + {0.6L, 1.21766907555422340e+16L}, {0.7L, 1.21535350115069900e+16L}, + {0.8L, 9.32311375306574375e+14L}, {0.9L, -8.99466171009315500e+15L}}; + Table[{50, 20}] = { + {0.0L, -7.87359356973322099e+32L}, {0.1L, 5.41513582908981306e+31L}, + {0.2L, 7.93092660564343958e+32L}, {0.3L, -4.00540672362456500e+31L}, + {0.4L, -8.23243252797740372e+32L}, {0.5L, -4.01698608142744523e+32L}, + {0.6L, 4.08601287568084080e+32L}, {0.7L, 8.60225211647169535e+32L}, + {0.8L, 1.06226766578920635e+33L}, {0.9L, 1.66309251586455094e+33L}}; + Table[{50, 50}] = { + {0.0L, 2.72539213975072988e+78L}, {0.1L, 2.11986821903666169e+78L}, + {0.2L, 9.82222379316811538e+77L}, {0.3L, 2.57907402241500578e+77L}, + {0.4L, 3.48667245334433154e+76L}, {0.5L, 2.05097602570372398e+75L}, + {0.6L, 3.88980964317819002e+73L}, {0.7L, 1.33254965595665771e+71L}, + {0.8L, 2.20288344031012791e+67L}, {0.9L, 2.53669949744314049e+60L}}; + Table[{100, 10}] = { + {0.0L, -8.24945971816703795e+18L}, {0.1L, 6.86278552250344038e+18L}, + {0.2L, -2.34688103580907827e+18L}, {0.3L, -4.90893583880523878e+18L}, + {0.4L, 8.23708346184265728e+18L}, {0.5L, 4.23871230219633664e+18L}, + {0.6L, -9.59361116591237120e+17L}, {0.7L, 4.43550484874968269e+18L}, + {0.8L, 1.57786735456734694e+18L}, {0.9L, 2.56433959576586016e+17L}}; + Table[{100, 20}] = { + {0.0L, 7.77274398361595813e+38L}, {0.1L, -7.05032664517025152e+38L}, + {0.2L, 4.39662513074438558e+38L}, {0.3L, 1.14744968919013248e+38L}, + {0.4L, -7.61205860438438888e+38L}, {0.5L, 3.07504970399995784e+38L}, + {0.6L, 8.26800055741207868e+38L}, {0.7L, 7.90222368531099936e+38L}, + {0.8L, -8.39638951169623065e+38L}, {0.9L, 7.16043448787788025e+37L}}; + Table[{100, 50}] = { + {0.0L, -1.17121450315783826e+98L}, {0.1L, 9.01193385501805844e+97L}, + {0.2L, -2.77911495881219591e+97L}, {0.3L, -3.35206020674796850e+97L}, + {0.4L, 7.14313850937400452e+97L}, {0.5L, -8.22526203397273908e+97L}, + {0.6L, 5.93640459256693667e+97L}, {0.7L, 1.80893339034652566e+97L}, + {0.8L, -1.10297979774543506e+98L}, {0.9L, 9.32312785168941904e+96L}}; + Table[{100, 100}] = { + {0.0L, 6.66630867007295431e+186L}, {0.1L, 4.03315719080569865e+186L}, + {0.2L, 8.65858791475270039e+185L}, {0.3L, 5.96973475268218189e+184L}, + {0.4L, 1.09106273304589198e+183L}, {0.5L, 3.77527496828897112e+180L}, + {0.6L, 1.35795105902891242e+177L}, {0.7L, 1.59365468505949483e+172L}, + {0.8L, 4.35522360415857506e+164L}, {0.9L, 5.77517922557585923e+150L}}; + compareWithTable(Table, RelTolerance, AbsTolerance); +} + +template +void compareWithGCC(T RelTolerance, T AbsTolerance) { + // Generated with gcc 8.3.1 + std::map, std::map > + Table; + Table[{10, 10}] = {{0.0L, 6.547290750000000000000e+08L}, + {0.1L, 6.226408357052308421698e+08L}, + {0.2L, 5.338482120799027200555e+08L}, + {0.3L, 4.085719891315887827659e+08L}, + {0.4L, 2.738155182015052800416e+08L}, + {0.5L, 1.553702785400390624272e+08L}, + {0.6L, 7.030099912163327995950e+07L}, + {0.7L, 2.258980634388728250269e+07L}, + {0.8L, 3.958896348126719998163e+06L}, + {0.9L, 1.621174007878425001223e+05L}}; + Table[{20, 10}] = {{0.0L, -1.612052956674316406250e+12L}, + {0.1L, 3.560417563907174735665e+11L}, + {0.2L, 1.466609267659881798506e+12L}, + {0.3L, -1.092420454297714176655e+12L}, + {0.4L, -8.999734873105529630780e+11L}, + {0.5L, 1.744502295946207409143e+12L}, + {0.6L, -4.244719151950563148558e+11L}, + {0.7L, -1.548364524949579275370e+12L}, + {0.8L, 2.392757933281052241802e+12L}, + {0.9L, 9.900174766949915345311e+11L}}; + Table[{20, 20}] = {{0.0L, 3.198309867728777708175e+23L}, + {0.1L, 2.892494114697687005594e+23L}, + {0.2L, 2.126340780079752328315e+23L}, + {0.3L, 1.245473413229778328863e+23L}, + {0.4L, 5.593883258401251200614e+22L}, + {0.5L, 1.801080697817960689050e+22L}, + {0.6L, 3.687400224900788126976e+21L}, + {0.7L, 3.807345588062060330880e+20L}, + {0.8L, 1.169352761683327100000e+19L}, + {0.9L, 1.960904971202388847266e+16L}}; + Table[{50, 10}] = {{0.0L, 1.145223824924634340332e+16L}, + {0.1L, 2.794368162749986069824e+15L}, + {0.2L, -9.926439446673560716797e+15L}, + {0.3L, -9.451384032851574123047e+15L}, + {0.4L, 8.892017018669203554688e+14L}, + {0.5L, 9.180035388465756383789e+15L}, + {0.6L, 1.217669075554222969727e+16L}, + {0.7L, 1.215353501150698340625e+16L}, + {0.8L, 9.323113753066150476685e+14L}, + {0.9L, -8.994661710093201598633e+15L}}; + Table[{50, 20}] = {{0.0L, -7.873593569733218616572e+32L}, + {0.1L, 5.415135829089829040302e+31L}, + {0.2L, 7.930926605643441888714e+32L}, + {0.3L, -4.005406723624505754186e+31L}, + {0.4L, -8.232432527977407811105e+32L}, + {0.5L, -4.016986081427444523960e+32L}, + {0.6L, 4.086012875680830742676e+32L}, + {0.7L, 8.602252116471695076665e+32L}, + {0.8L, 1.062267665789207148447e+33L}, + {0.9L, 1.663092515864551144796e+33L}}; + Table[{50, 50}] = {{0.0L, 2.725392139750729502972e+78L}, + {0.1L, 2.119868219036661812206e+78L}, + {0.2L, 9.822223793168123441126e+77L}, + {0.3L, 2.579074022415003126022e+77L}, + {0.4L, 3.486672453344334274446e+76L}, + {0.5L, 2.050976025703723878091e+75L}, + {0.6L, 3.889809643178187639203e+73L}, + {0.7L, 1.332549655956656855416e+71L}, + {0.8L, 2.202883440310146416308e+67L}, + {0.9L, 2.536699497443158134445e+60L}}; + Table[{100, 10}] = {{0.0L, -8.249459718167038282000e+18L}, + {0.1L, 6.862785522503444437500e+18L}, + {0.2L, -2.346881035809095210750e+18L}, + {0.3L, -4.908935838805249369500e+18L}, + {0.4L, 8.237083461842661207000e+18L}, + {0.5L, 4.238712302196339605250e+18L}, + {0.6L, -9.593611165912185080625e+17L}, + {0.7L, 4.435504848749747270500e+18L}, + {0.8L, 1.577867354567400466500e+18L}, + {0.9L, 2.564339595766638741406e+17L}}; + Table[{100, 20}] = {{0.0L, 7.772743983615954055061e+38L}, + {0.1L, -7.050326645170251933206e+38L}, + {0.2L, 4.396625130744397503342e+38L}, + {0.3L, 1.147449689190144464586e+38L}, + {0.4L, -7.612058604384389981878e+38L}, + {0.5L, 3.075049703999956739724e+38L}, + {0.6L, 8.268000557412077024506e+38L}, + {0.7L, 7.902223685310973476875e+38L}, + {0.8L, -8.396389511696198700493e+38L}, + {0.9L, 7.160434487878827530346e+37L}}; + Table[{100, 50}] = {{0.0L, -1.171214503157837848611e+98L}, + {0.1L, 9.011933855018056106092e+97L}, + {0.2L, -2.779114958812187282420e+97L}, + {0.3L, -3.352060206747950450847e+97L}, + {0.4L, 7.143138509374017832078e+97L}, + {0.5L, -8.225262033972737461155e+97L}, + {0.6L, 5.936404592566909963016e+97L}, + {0.7L, 1.808933390346613713076e+97L}, + {0.8L, -1.102979797745434149687e+98L}, + {0.9L, 9.323127851689438561522e+96L}}; + Table[{100, 100}] = {{0.0L, 6.666308670072953743181e+186L}, + {0.1L, 4.033157190805700092363e+186L}, + {0.2L, 8.658587914752715633363e+185L}, + {0.3L, 5.969734752682171415182e+184L}, + {0.4L, 1.091062733045894109826e+183L}, + {0.5L, 3.775274968288970779608e+180L}, + {0.6L, 1.357951059028910874245e+177L}, + {0.7L, 1.593654685059493405648e+172L}, + {0.8L, 4.355223604158650443018e+164L}, + {0.9L, 5.775179225575939991155e+150L}}; + compareWithTable(Table, RelTolerance, AbsTolerance); +} + +template +void compareWithBoost(T RelTolerance, T AbsTolerance) { + // Generated with gcc 8.3.1 and boost 1.66.0 + std::map, std::map > + Table; + Table[{10, 10}] = {{0.0L, 6.547290750000000000000e+08L}, + {0.1L, 6.226408357052308425191e+08L}, + {0.2L, 5.338482120799027199391e+08L}, + {0.3L, 4.085719891315887824458e+08L}, + {0.4L, 2.738155182015052799543e+08L}, + {0.5L, 1.553702785400390625000e+08L}, + {0.6L, 7.030099912163327999588e+07L}, + {0.7L, 2.258980634388728249723e+07L}, + {0.8L, 3.958896348126719997708e+06L}, + {0.9L, 1.621174007878425002218e+05L}}; + Table[{20, 10}] = {{0.0L, -1.612052956674316406250e+12L}, + {0.1L, 3.560417563907174737155e+11L}, + {0.2L, 1.466609267659881798267e+12L}, + {0.3L, -1.092420454297714175999e+12L}, + {0.4L, -8.999734873105529626012e+11L}, + {0.5L, 1.744502295946207410216e+12L}, + {0.6L, -4.244719151950563156009e+11L}, + {0.7L, -1.548364524949579275131e+12L}, + {0.8L, 2.392757933281052242041e+12L}, + {0.9L, 9.900174766949915348291e+11L}}; + Table[{20, 20}] = {{0.0L, 3.198309867728777708175e+23L}, + {0.1L, 2.892494114697687008379e+23L}, + {0.2L, 2.126340780079752327168e+23L}, + {0.3L, 1.245473413229778327224e+23L}, + {0.4L, 5.593883258401251196928e+22L}, + {0.5L, 1.801080697817960690381e+22L}, + {0.6L, 3.687400224900788130304e+21L}, + {0.7L, 3.807345588062060328960e+20L}, + {0.8L, 1.169352761683327099800e+19L}, + {0.9L, 1.960904971202388849219e+16L}}; + Table[{50, 10}] = {{0.0L, 1.145223824924634340332e+16L}, + {0.1L, 2.794368162749986070557e+15L}, + {0.2L, -9.926439446673560716797e+15L}, + {0.3L, -9.451384032851574119141e+15L}, + {0.4L, 8.892017018669203523560e+14L}, + {0.5L, 9.180035388465756391602e+15L}, + {0.6L, 1.217669075554222970703e+16L}, + {0.7L, 1.215353501150698341113e+16L}, + {0.8L, 9.323113753066150515747e+14L}, + {0.9L, -8.994661710093201604004e+15L}}; + Table[{50, 20}] = {{0.0L, -7.873593569733218616572e+32L}, + {0.1L, 5.415135829089829046900e+31L}, + {0.2L, 7.930926605643441882380e+32L}, + {0.3L, -4.005406723624505746709e+31L}, + {0.4L, -8.232432527977407807587e+32L}, + {0.5L, -4.016986081427444519738e+32L}, + {0.6L, 4.086012875680830746546e+32L}, + {0.7L, 8.602252116471695071036e+32L}, + {0.8L, 1.062267665789207148447e+33L}, + {0.9L, 1.663092515864551146625e+33L}}; + Table[{50, 50}] = {{0.0L, 2.725392139750729502972e+78L}, + {0.1L, 2.119868219036661818834e+78L}, + {0.2L, 9.822223793168123432087e+77L}, + {0.3L, 2.579074022415003118238e+77L}, + {0.4L, 3.486672453344334267855e+76L}, + {0.5L, 2.050976025703723882014e+75L}, + {0.6L, 3.889809643178187647171e+73L}, + {0.7L, 1.332549655956656853859e+71L}, + {0.8L, 2.202883440310146415724e+67L}, + {0.9L, 2.536699497443158141239e+60L}}; + Table[{100, 10}] = {{0.0L, -8.249459718167038282000e+18L}, + {0.1L, 6.862785522503444443500e+18L}, + {0.2L, -2.346881035809095210750e+18L}, + {0.3L, -4.908935838805249365500e+18L}, + {0.4L, 8.237083461842661209000e+18L}, + {0.5L, 4.238712302196339609250e+18L}, + {0.6L, -9.593611165912185088125e+17L}, + {0.7L, 4.435504848749747265500e+18L}, + {0.8L, 1.577867354567400459375e+18L}, + {0.9L, 2.564339595766638627500e+17L}}; + Table[{100, 20}] = {{0.0L, 7.772743983615954055061e+38L}, + {0.1L, -7.050326645170251938371e+38L}, + {0.2L, 4.396625130744397503342e+38L}, + {0.3L, 1.147449689190144463479e+38L}, + {0.4L, -7.612058604384389980402e+38L}, + {0.5L, 3.075049703999956748025e+38L}, + {0.6L, 8.268000557412077031146e+38L}, + {0.7L, 7.902223685310973468759e+38L}, + {0.8L, -8.396389511696198700493e+38L}, + {0.9L, 7.160434487878827394302e+37L}}; + Table[{100, 50}] = {{0.0L, -1.171214503157837848611e+98L}, + {0.1L, 9.011933855018056135735e+97L}, + {0.2L, -2.779114958812187282420e+97L}, + {0.3L, -3.352060206747950445289e+97L}, + {0.4L, 7.143138509374017815033e+97L}, + {0.5L, -8.225262033972737473753e+97L}, + {0.6L, 5.936404592566909988953e+97L}, + {0.7L, 1.808933390346613708259e+97L}, + {0.8L, -1.102979797745434147612e+98L}, + {0.9L, 9.323127851689438586533e+96L}}; + Table[{100, 100}] = {{0.0L, 6.666308670072953744596e+186L}, + {0.1L, 4.033157190805700115242e+186L}, + {0.2L, 8.658587914752715622749e+185L}, + {0.3L, 5.969734752682171379064e+184L}, + {0.4L, 1.091062733045894106141e+183L}, + {0.5L, 3.775274968288970794454e+180L}, + {0.6L, 1.357951059028910880506e+177L}, + {0.7L, 1.593654685059493402464e+172L}, + {0.8L, 4.355223604158650442269e+164L}, + {0.9L, 5.775179225575940023096e+150L}}; + compareWithTable(Table, RelTolerance, AbsTolerance); +} + +template +void testReferenceValues(T RelTolerance, T AbsTolerance) { + compareWithMSVC(RelTolerance, AbsTolerance); + compareWithGCC(RelTolerance, AbsTolerance); + compareWithBoost(RelTolerance, AbsTolerance); +} + +int main(int, char**) { + compareWithAnalyticExpressions(1e-6f, 1e-6f); + compareWithAnalyticExpressions(1e-12, 1e-12); + + testReferenceValues(1e-6f, 0.0f); + testReferenceValues(1e-12, 0.0); +} diff --git a/libcxx/test/std/numerics/legendre/assoc_legendre_exceptions.pass.cpp b/libcxx/test/std/numerics/legendre/assoc_legendre_exceptions.pass.cpp new file mode 100644 --- /dev/null +++ b/libcxx/test/std/numerics/legendre/assoc_legendre_exceptions.pass.cpp @@ -0,0 +1,48 @@ +//===----------------------------------------------------------------------===// +// +// 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++98, c++03, c++11, c++14, libcpp-no-exceptions + +// + +// floating_point assoc_legendre (unsigned l, unsigned m, arithmetic x); +// float assoc_legendref(unsigned l, unsigned m, float x); +// long double assoc_legendrel(unsigned l, unsigned m, long double x); + +#include +#include +#include + +/// \brief std::legendre(n,x) throw if x is not in the range [-1,1] and not NaN. +/// +/// \tparam T one of float, double, or long double +template +void testInvalidArguments() { + const T InvalidArguments[] = {T(-5), T(2)}; + for (unsigned l = 0; l < 128; ++l) { + for (unsigned m = 0; m < 128; ++m) { + for (T x : InvalidArguments) { + bool ThrowsDomainError = false; + try { + const T Value = std::assoc_legendre(l, m, x); + (void)Value; + } catch (const std::domain_error&) { + ThrowsDomainError = true; + } catch (...) { + } + assert(ThrowsDomainError); + } + } + } +} + +int main(int, char**) { + testInvalidArguments(); + testInvalidArguments(); + testInvalidArguments(); +} diff --git a/libcxx/test/std/numerics/legendre/legendre_basic.pass.cpp b/libcxx/test/std/numerics/legendre/legendre_basic.pass.cpp new file mode 100644 --- /dev/null +++ b/libcxx/test/std/numerics/legendre/legendre_basic.pass.cpp @@ -0,0 +1,136 @@ +//===----------------------------------------------------------------------===// +// +// 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++98, c++03, c++11, c++14 + +// + +// floating_point legendre (unsigned n, arithmetic x); +// float legendref(unsigned n, float x); +// long double legendrel(unsigned n, long double x); + +#include +#include +#include + +/// \brief Tests if the Legendre polynomial behaves correctly +/// under the transformation x -> -x +/// +/// The expected property is +/// \f$ P_{n}(-x) = (-1)^{n} P_{n}(x) \f$. +/// For odd \f$ n \f$, this test checks for a root at \f$ x = 0\f$. +/// \tparam T one of float, double, or long double +/// \param [in] n degree of the Legendre polynomial +/// \param [in] x argument of the Legendre polynomial +/// \param [in] RelTolerance relative tolerance for the value of legendre +/// \param [in] AbsTolerance absolute tolerance for the value of legendre +template +void testParity(unsigned n, T x, T RelTolerance, T AbsTolerance) { + const T Parity = n % 2 == 0 ? T(1) : T(-1); + const T ExpectedValue = Parity * std::legendre(n, x); + const T Tolerance = RelTolerance * std::abs(ExpectedValue) + AbsTolerance; + const T Error = std::abs(std::legendre(n, -x) - ExpectedValue); + assert(Error <= Tolerance); +} + +template +void testParity(const T RelTolerance, const T AbsTolerance) { + const T Samples[] = {T(0), T(0.1), T(0.2), T(0.5), T(0.9), T(1.0)}; + for (unsigned n = 0; n < 128; ++n) { + for (T x : Samples) + testParity(n, x, RelTolerance, AbsTolerance); + } +} + +/// \brief Tests if the values of the Legendre polynomial at x = -1 +/// and x = -1 are correct. +/// +/// The expected values are: +/// \f$ P_{n}(1) = 1 \f$ and \f$ P_{n}(1) = (-1)^{n} \f$. +/// These properties can be easily proven by induction. +/// \tparam T one of float, double, or long double +/// \param [in] n degree of the legendre polynomial. +/// \param [in] RelTolerance relative tolerance for the value of legendre(x) +/// \param [in] AbsTolerance absolute tolerance for the value of legendre(x) +template +void testBoundaryValues(unsigned n, T RelTolerance, T AbsTolerance) { + const T Parity = n % 2 == 0 ? T(1) : T(-1); + const T Tolerance = RelTolerance + AbsTolerance; + const T ErrorRight = std::abs(std::legendre(n, T(1)) - T(1)); + const T ErrorLeft = std::abs(std::legendre(n, T(-1)) - Parity); + assert(ErrorRight <= Tolerance); + assert(ErrorLeft <= Tolerance); +} + +template +void testBoundaryValues(const T RelTolerance, const T AbsTolerance) { + for (unsigned n = 0; n < 128; ++n) + testBoundaryValues(n, RelTolerance, AbsTolerance); +} + +/// \brief Tests if an analytic upper bound holds +/// for the implementation of legendre. +/// +/// The upper bound is given by +/// \f[ +/// \left| P_n (x) \right| \leq +/// \left( 1 + n(n+1) (1-x^2)\right)^{-\frac{1}{4}}. +/// \f] +/// In particular, we check for \f$ \left| P_n (x) \right| \leq 1 \f$. +/// See Eq. A1, p. 1435, in Martin, A. (1963). +/// Unitarity and high-energy behavior of scattering amplitudes. +/// Physical Review, 129(3), 1432. +/// \tparam T one of float, double, or long double +/// \param [in] n degree of the Legendre polynomial +/// \param [in] x argument of the Legendre polynomial +/// \param [in] RelTolerance relative tolerance for the value of legendre +/// \param [in] AbsTolerance absolute tolerance for the value of legendre +template +void testUpperBounds(unsigned n, T x, T RelTolerance, T AbsTolerance) { + const T UpperBound = + T(1) / + std::sqrt(std::sqrt(T(1) + T(n) * T(n + 1u) * (T(1) - x) * (T(1) + x))); + const T LegendreAbs = std::abs(std::legendre(n, x)); + const T UpperBoundWithTolerance = + UpperBound + UpperBound * RelTolerance + AbsTolerance; + assert(LegendreAbs <= UpperBoundWithTolerance); +} + +template +void testUpperBounds(T RelTolerance, T AbsTolerance) { + const T Samples[] = {T(0), T(0.1), T(0.2), T(0.5), T(0.9), T(1.0)}; + for (unsigned n = 0; n < 128; ++n) { + for (T x : Samples) + testUpperBounds(n, x, RelTolerance, AbsTolerance); + } +} + +/// \brief std::legendre(n,x) must silently return NaN if x is NaN. +/// +/// \tparam T one of float, double, or long double +template +void testNaNPropagation() { + for (unsigned n = 0; n < 128; ++n) { + const T ShouldBeNaN = std::legendre(n, std::numeric_limits::quiet_NaN()); + assert(std::isnan(ShouldBeNaN)); + } +} + +int main(int, char**) { + testParity(1e-6f, 0.0f); + testParity(1e-12, 0.0); + + testBoundaryValues(1e-6f, 0.0f); + testBoundaryValues(1e-12, 0.0); + + testUpperBounds(1e-6f, 0.0f); + testUpperBounds(1e-12, 0.0); + + testNaNPropagation(); + testNaNPropagation(); +} diff --git a/libcxx/test/std/numerics/legendre/legendre_comparisons.pass.cpp b/libcxx/test/std/numerics/legendre/legendre_comparisons.pass.cpp new file mode 100644 --- /dev/null +++ b/libcxx/test/std/numerics/legendre/legendre_comparisons.pass.cpp @@ -0,0 +1,315 @@ +//===----------------------------------------------------------------------===// +// +// 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++98, c++03, c++11, c++14 + +// + +// floating_point legendre (unsigned n, arithmetic x); +// float legendref(unsigned n, float x); +// long double legendrel(unsigned n, long double x); + +#include +#include +#include +#include + +template +Real evaluatePolynomial(Iterator CoefficientsBegin, Iterator CoefficientsEnd, + Real x) { + Real Result = Real(0); + for (Iterator Iter = CoefficientsBegin; Iter != CoefficientsEnd; ++Iter) { + Result = Result * x + *Iter; + } + return Result; +} + +template +Real evaluatePolynomial(const Container& Coefficients, Real x) { + return evaluatePolynomial(std::begin(Coefficients), std::end(Coefficients), + x); +} + +/// \brief Checks if the Legendre polynomials up to n=16 +/// have roots at values from the literature. +/// +/// The test checks for a change of sign in the interval that is +/// centered around the expected root and has the width of the +/// provided tolerance. +/// The expected values were taken from: +/// Lowan, A. N., Davids, N., & Levenson, A. (1942). +/// Table of the zeros of the Legendre polynomials of order 1-16 and +/// the weight coefficients for Gauss� mechanical quadrature formula. +/// Bulletin of the American Mathematical Society, 48(10), 739-743. +/// For n=11 and n=12, I had to correct some misprints (in accordance +/// with wolframalpha). +/// \tparam T one of float, double, or long double +/// \param [in] RelTolerance Tolerance for x; contributes to the +/// width of the interval where legendre changes its sign. +/// \param [in] AbsTolerance Tolerance for x; contributes +/// to the +/// width of the interval where legendre changes its sign. +template +void compareZerosWithLowan1942(T RelTolerance, T AbsTolerance) { + std::vector > LegendreZeros = { + {}, + {T(0)}, + {T(0.577350269189626L)}, + {T(0), T(0.774596669241483)}, + {T(0.339981043584856), T(0.86113631159405)}, + {T(0), T(0.538469310105683), T(0.90617984593866)}, + {T(0.238619186083197), T(0.661209386466265), T(0.93246951420315)}, + {T(0), T(0.405845151377397), T(0.741531185599394), T(0.949107912342759)}, + {T(0.183434642495650), T(0.525532409916329), T(0.796666477413627), + T(0.960289856497536)}, + {T(0), T(0.324253423403809), T(0.613371432700590), T(0.836031107326636), + T(0.968160239507626)}, + {T(0.148874338981631), T(0.433395394129247), T(0.679409568299024), + T(0.865063366688985), T(0.973906528517172)}, + {T(0), T(0.269543155952345), T(0.51909612920681), T(0.730152005574049), + T(0.887062599768095), T(0.978228658146057)}, + {T(0.125233408511469), T(0.367831498998180), T(0.587317954286617), + T(0.769902674194305), T(0.904117256370475), T(0.981560634246719)}, + {T(0), T(0.230458315955135), T(0.448492751036447), T(0.642349339440340), + T(0.801578090733310), T(0.917598399222978), T(0.984183054718588)}, + {T(0.108054948707344), T(0.319112368927890), T(0.515248636358154), + T(0.687292904811685), T(0.827201315069765), T(0.928434883663574), + T(0.986283808696812)}, + {T(0), T(0.201194093997435), T(0.394151347077563), T(0.570972172608539), + T(0.724417731360170), T(0.848206583410427), T(0.937273392400706), + T(0.987992518020485)}, + {T(0.095012509837637), T(0.281603550779259), T(0.458016777657227), + T(0.617876244402644), T(0.755404408355003), T(0.865631202387832), + T(0.944575023073233), T(0.989400934991650)}}; + for (size_t n = 0; n < LegendreZeros.size(); ++n) { + for (T Root : LegendreZeros[n]) { + const T HalfWidthInterval = T(0.5) * (RelTolerance * Root + AbsTolerance); + const T LeftOfRoot = Root - HalfWidthInterval; + const T RightOfRoot = Root + HalfWidthInterval; + assert(LeftOfRoot != Root); + assert(RightOfRoot != Root); + const T LegendreAtLeft = + std::legendre(static_cast(n), LeftOfRoot); + const T LegendreAtRight = + std::legendre(static_cast(n), RightOfRoot); + assert(std::signbit(LegendreAtLeft) ^ std::signbit(LegendreAtRight)); + } + } +} + +/// \brief Compares std::legendre with analytic expressions +/// for n = 0..6. +/// +/// \tparam T one of float, double, or long double +/// \tparam TestFunction a callable with arguments (unsigned,T,T) +/// \param [in] x argument of legendre +/// \param [in] ValueSink will be invoked with the arguments of +/// legendre and its expected outcome, i.e., (n,x,expected). +template +void compareWithAnalyticExpressions(T x, TestFunction&& ValueSink) { + std::vector > VectorOfCoeffients; + VectorOfCoeffients.push_back({T(1)}); + VectorOfCoeffients.push_back({T(1), T(0)}); + VectorOfCoeffients.push_back({T(3) / T(2), T(0), -T(1) / T(2)}); + VectorOfCoeffients.push_back({T(5) / T(2), T(0), -T(3) / T(2), T(0)}); + VectorOfCoeffients.push_back( + {T(35) / T(8), T(0), -T(30) / T(8), T(0), T(3) / T(8)}); + VectorOfCoeffients.push_back( + {T(63) / T(8), T(0), -T(70) / T(8), T(0), T(15) / T(8), T(0)}); + VectorOfCoeffients.push_back({T(231) / T(16), T(0), -T(315) / T(16), T(0), + T(105) / T(16), T(0), -T(5) / T(16)}); + + for (const auto& Coefficients : VectorOfCoeffients) { + const unsigned int n = static_cast(Coefficients.size()) - 1; + const T ExpectedValue = ::evaluatePolynomial(Coefficients, x); + ValueSink(n, x, ExpectedValue); + } +} + +/// \brief Compares std::legendre with analytic expressions +/// for n = 0..6. +/// +/// \tparam T one of float, double, or long double +/// \param [in] RelTolerance relative tolerance for the value of legendre +/// \param [in] AbsTolerance absolute tolerance for the value of legendre +template +void compareWithAnalyticExpressions(T RelTolerance, T AbsTolerance) { + auto testFunction = [RelTolerance, AbsTolerance](unsigned n, T x, + T ExpectedValue) { + const T Value = std::legendre(n, x); + const T Tolerance = std::abs(ExpectedValue) * RelTolerance + AbsTolerance; + const T Error = std::abs(Value - ExpectedValue); + assert(Error <= Tolerance); + }; + + const auto Samples = {T(0), T(0.1), T(0.2), T(0.5), T(0.9), T(1.0)}; + for (T x : Samples) + compareWithAnalyticExpressions(x, testFunction); +} + +template +void compareWithTable( + const std::map >& Table, + T RelTolerance, T AbsTolerance) { + for (const auto& Entry : Table) { + const unsigned n = Entry.first; + for (const auto& xy : Entry.second) { + const T Value = std::legendre(n, xy.first); + const T ExpectedValue = xy.second; + const T Tolerance = std::abs(ExpectedValue) * RelTolerance + AbsTolerance; + const T Error = std::abs(Value - ExpectedValue); + assert(Error <= Tolerance); + } + } +} + +template +void compareWithMSVC(T RelTolerance, T AbsTolerance) { + // Generated with msvc 191627027 + std::map > Table; + Table[10u] = { + {0.0L, -2.46093750000000000e-01L}, {0.1L, -1.22124997387109330e-01L}, + {0.2L, 1.29072025600000029e-01L}, {0.3L, 2.51476349516015663e-01L}, + {0.4L, 9.68390643999999245e-02L}, {0.5L, -1.88228607177734375e-01L}, + {0.6L, -2.43662745600000008e-01L}, {0.7L, 8.58057955316403331e-02L}, + {0.8L, 3.00529795600000038e-01L}, {0.9L, -2.63145617855859604e-01L}}; + Table[20u] = { + {0.0L, 1.76197052001953125e-01L}, {0.1L, -8.20771309445277736e-02L}, + {0.2L, -9.80421943445946159e-02L}, {0.3L, 1.80287159479980419e-01L}, + {0.4L, -1.01592615586281562e-01L}, {0.5L, -4.83583810673735570e-02L}, + {0.6L, 1.59167529100981087e-01L}, {0.7L, -2.04573944638341720e-01L}, + {0.8L, 2.24204605417413466e-01L}, {0.9L, -1.49308235309848353e-01L}}; + Table[50u] = { + {0.0L, -1.12275172659217048e-01L}, {0.1L, -3.82058126613136831e-02L}, + {0.2L, 8.34322722041973136e-02L}, {0.3L, 1.09110515747147960e-01L}, + {0.4L, 4.15690333818253752e-02L}, {0.5L, -3.10590992396098109e-02L}, + {0.6L, -5.88607988440020963e-02L}, {0.7L, -1.45727316458928516e-02L}, + {0.8L, 1.38797373450931127e-01L}, {0.9L, -1.70037659943836711e-01L}}; + Table[100u] = { + {0.0L, 7.95892373871787268e-02L}, {0.1L, -6.38950984347501499e-02L}, + {0.2L, 1.46818353556592614e-02L}, {0.3L, 5.71273922028013650e-02L}, + {0.4L, -7.22582021256844703e-02L}, {0.5L, -6.05180259618611979e-02L}, + {0.6L, -2.37470239051331411e-02L}, {0.7L, -7.71325071997786410e-02L}, + {0.8L, 5.08611679135842279e-02L}, {0.9L, 1.02265820558718926e-01L}}; + + compareWithTable(Table, RelTolerance, AbsTolerance); +} + +template +void compareWithGCC(T RelTolerance, T AbsTolerance) { + // Generated with gcc 8.3.1 + std::map > Table; + Table[10u] = {{0.0L, -2.460937500000000000000e-01L}, + {0.1L, -1.221249973871093750063e-01L}, + {0.2L, 1.290720256000000000328e-01L}, + {0.3L, 2.514763495160156249876e-01L}, + {0.4L, 9.683906439999999996141e-02L}, + {0.5L, -1.882286071777343750000e-01L}, + {0.6L, -2.436627456000000000027e-01L}, + {0.7L, 8.580579553164062492864e-02L}, + {0.8L, 3.005297956000000000312e-01L}, + {0.9L, -2.631456178558593748918e-01L}}; + Table[20u] = {{0.0L, 1.761970520019531250000e-01L}, + {0.1L, -8.207713094452772853985e-02L}, + {0.2L, -9.804219434459463679125e-02L}, + {0.3L, 1.802871594799804655558e-01L}, + {0.4L, -1.015926155862814718285e-01L}, + {0.5L, -4.835838106737355701625e-02L}, + {0.6L, 1.591675291009810431662e-01L}, + {0.7L, -2.045739446383416052931e-01L}, + {0.8L, 2.242046054174133682042e-01L}, + {0.9L, -1.493082353098486828345e-01L}}; + Table[50u] = {{0.0L, -1.122751726592170484764e-01L}, + {0.1L, -3.820581266131368189813e-02L}, + {0.2L, 8.343227220419731034374e-02L}, + {0.3L, 1.091105157471479767954e-01L}, + {0.4L, 4.156903338182522686674e-02L}, + {0.5L, -3.105909923960982286266e-02L}, + {0.6L, -5.886079884400199673851e-02L}, + {0.7L, -1.457273164589237150986e-02L}, + {0.8L, 1.387973734509308074478e-01L}, + {0.9L, -1.700376599438367846763e-01L}}; + Table[100u] = {{0.0L, 7.958923738717876149887e-02L}, + {0.1L, -6.389509843475016413679e-02L}, + {0.2L, 1.468183535565932807048e-02L}, + {0.3L, 5.712739220280141857569e-02L}, + {0.4L, -7.225820212568462422759e-02L}, + {0.5L, -6.051802596186118686668e-02L}, + {0.6L, -2.374702390513330860656e-02L}, + {0.7L, -7.713250719977911342573e-02L}, + {0.8L, 5.086116791358337655272e-02L}, + {0.9L, 1.022658205587185505935e-01L}}; + + compareWithTable(Table, RelTolerance, AbsTolerance); +} + +template +void compareWithBoost(T RelTolerance, T AbsTolerance) { + // Generated with gcc 8.3.1 and boost 1.66.0 + std::map > Table; + Table[10u] = {{0.0L, -2.460937500000000000000e-01L}, + {0.1L, -1.221249973871093750063e-01L}, + {0.2L, 1.290720256000000000057e-01L}, + {0.3L, 2.514763495160156249876e-01L}, + {0.4L, 9.683906439999999996818e-02L}, + {0.5L, -1.882286071777343750000e-01L}, + {0.6L, -2.436627455999999999756e-01L}, + {0.7L, 8.580579553164062490831e-02L}, + {0.8L, 3.005297955999999999499e-01L}, + {0.9L, -2.631456178558593748647e-01L}}; + Table[20u] = {{0.0L, 1.761970520019531250000e-01L}, + {0.1L, -8.207713094452772851275e-02L}, + {0.2L, -9.804219434459463677092e-02L}, + {0.3L, 1.802871594799804655558e-01L}, + {0.4L, -1.015926155862814718082e-01L}, + {0.5L, -4.835838106737355701625e-02L}, + {0.6L, 1.591675291009810431391e-01L}, + {0.7L, -2.045739446383416053744e-01L}, + {0.8L, 2.242046054174133682177e-01L}, + {0.9L, -1.493082353098486828480e-01L}}; + Table[50u] = {{0.0L, -1.122751726592170484764e-01L}, + {0.1L, -3.820581266131368188457e-02L}, + {0.2L, 8.343227220419731035052e-02L}, + {0.3L, 1.091105157471479767683e-01L}, + {0.4L, 4.156903338182522687013e-02L}, + {0.5L, -3.105909923960982287282e-02L}, + {0.6L, -5.886079884400199674528e-02L}, + {0.7L, -1.457273164589237145396e-02L}, + {0.8L, 1.387973734509308073665e-01L}, + {0.9L, -1.700376599438367847305e-01L}}; + Table[100u] = {{0.0L, 7.958923738717876149209e-02L}, + {0.1L, -6.389509843475016413679e-02L}, + {0.2L, 1.468183535565932810267e-02L}, + {0.3L, 5.712739220280141854181e-02L}, + {0.4L, -7.225820212568462420048e-02L}, + {0.5L, -6.051802596186118688023e-02L}, + {0.6L, -2.374702390513330857607e-02L}, + {0.7L, -7.713250719977911350027e-02L}, + {0.8L, 5.086116791358337655272e-02L}, + {0.9L, 1.022658205587185507019e-01L}}; + + compareWithTable(Table, RelTolerance, AbsTolerance); +} + +template +void testReferenceValues(T RelTolerance, T AbsTolerance) { + compareWithMSVC(RelTolerance, AbsTolerance); + compareWithGCC(RelTolerance, AbsTolerance); + compareWithBoost(RelTolerance, AbsTolerance); +} + +int main(int, char**) { + compareZerosWithLowan1942(0.0f, 1e-6f); + compareZerosWithLowan1942(0.0f, 1e-12); + + compareWithAnalyticExpressions(1e-6f, 1e-6f); + compareWithAnalyticExpressions(1e-12, 1e-12); + + testReferenceValues(1e-6f, 0.0f); + testReferenceValues(1e-12, 0.0); +} diff --git a/libcxx/test/std/numerics/legendre/legendre_exceptions.pass.cpp b/libcxx/test/std/numerics/legendre/legendre_exceptions.pass.cpp new file mode 100644 --- /dev/null +++ b/libcxx/test/std/numerics/legendre/legendre_exceptions.pass.cpp @@ -0,0 +1,46 @@ +//===----------------------------------------------------------------------===// +// +// 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++98, c++03, c++11, c++14, libcpp-no-exceptions + +// + +// floating_point legendre (unsigned n, arithmetic x); +// float legendref(unsigned n, float x); +// long double legendrel(unsigned n, long double x); + +#include +#include +#include + +/// \brief std::legendre(n,x) must throw if x is not in the range [-1,1] and not NaN. +/// +/// \tparam T one of float, double, or long double +template +void testInvalidArguments() { + const T InvalidArguments[] = {T(-5), T(2)}; + for (unsigned n = 0; n < 128; ++n) { + for (T x : InvalidArguments) { + bool ThrowsDomainError = false; + try { + const T Value = std::legendre(n, x); + (void)Value; + } catch (const std::domain_error&) { + ThrowsDomainError = true; + } catch (...) { + } + assert(ThrowsDomainError); + } + } +} + +int main(int, char**) { + testInvalidArguments(); + testInvalidArguments(); + testInvalidArguments(); +}