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,8 @@ #include <__config> #include #include +#include +#include #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER) #pragma GCC system_header @@ -606,6 +618,214 @@ return isfinite(__lcpp_x); } + +#if _LIBCPP_STD_VER > 14 + +template struct __promote_floating_point {}; +template <> struct __promote_floating_point { typedef double type; }; +template <> struct __promote_floating_point { typedef long double type; }; +template <> struct __promote_floating_point { typedef long double type; }; + +/// \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(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(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 for float 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. +template +inline _LIBCPP_INLINE_VISIBILITY +typename __lazy_enable_if::value, __promote<_A1>>::type +assoc_legendre(unsigned __lcpp_l, unsigned __lcpp_m, _A1 __lcpp_x) { + typedef typename __promote<_A1>::type __result_type; + typedef typename __promote_floating_point<__result_type>::type __internal_type; + + return (__result_type) __libcpp_assoc_legendre<__internal_type>( + __lcpp_l, __lcpp_m, (__internal_type) __lcpp_x); +} + +inline _LIBCPP_INLINE_VISIBILITY +float assoc_legendref(unsigned __lcpp_l, unsigned __lcpp_m, + float __lcpp_x) { + return assoc_legendre(__lcpp_l, __lcpp_m, __lcpp_x); +} + +inline _LIBCPP_INLINE_VISIBILITY +long double assoc_legendrel(unsigned __lcpp_l, unsigned __lcpp_m, + long double __lcpp_x) { + return assoc_legendre(__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). +template +inline _LIBCPP_INLINE_VISIBILITY +typename __lazy_enable_if::value, __promote<_A1>>::type +legendre(unsigned __lcpp_n, _A1 __lcpp_x) { + typedef typename __promote<_A1>::type __result_type; + typedef typename __promote_floating_point<__result_type>::type __internal_type; + return (__result_type)__libcpp_legendre<__internal_type>(__lcpp_n, + (__internal_type)__lcpp_x); +} + +inline _LIBCPP_INLINE_VISIBILITY +float legendref(unsigned __lcpp_n, float __lcpp_x) { + return legendre(__lcpp_n, __lcpp_x); +} + +inline _LIBCPP_INLINE_VISIBILITY +long double legendrel(unsigned __lcpp_n, long double __lcpp_x) { + return legendre(__lcpp_n, __lcpp_x); +} + +#endif + _LIBCPP_END_NAMESPACE_STD #endif // _LIBCPP_CMATH diff --git a/libcxx/test/std/numerics/c.math/cmath.pass.cpp b/libcxx/test/std/numerics/c.math/cmath.pass.cpp --- a/libcxx/test/std/numerics/c.math/cmath.pass.cpp +++ b/libcxx/test/std/numerics/c.math/cmath.pass.cpp @@ -100,6 +100,11 @@ Ambiguous tgamma(Ambiguous){ return Ambiguous(); } Ambiguous trunc(Ambiguous){ return Ambiguous(); } +#if _LIBCPP_STD_VER > 14 +Ambiguous assoc_legendre(unsigned, unsigned, Ambiguous){ return Ambiguous(); } +Ambiguous legendre(unsigned, Ambiguous){ return Ambiguous(); } +#endif + void test_abs() { static_assert((std::is_same::value), ""); @@ -1514,6 +1519,47 @@ assert(std::trunc(1) == 1); } +#if _LIBCPP_STD_VER > 14 + +void test_assoc_legendre() +{ + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); +} + +void test_legendre() +{ + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); +} + +#endif + + int main(int, char**) { test_abs(); @@ -1586,6 +1632,11 @@ test_scalbn(); test_tgamma(); test_trunc(); + +#if _LIBCPP_STD_VER > 14 + test_assoc_legendre(); + test_legendre(); +#endif return 0; } 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,128 @@ +//===----------------------------------------------------------------------===// +// +// 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 long double Pi = 3.14159265358979323846264338327950288L; + const long double LogUpperBound = + 0.5L*std::lgammal(l + m + 1u) - 0.5L*std::lgammal(l - m + 1u) - + 0.25L*std::logl(m) + 1.25L * std::logl(2) - 0.75L * std::logl(Pi); + const T UpperBound = static_cast(std::expl(LogUpperBound)); + if (std::isinf(UpperBound)) + return; + + const T UpperBoundWithTolerance = + UpperBound + UpperBound * RelTolerance + AbsTolerance; + const T LegendreAbs = std::abs(std::assoc_legendre(l, m, x)); + assert(LegendreAbs <= 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,356 @@ +//===----------------------------------------------------------------------===// +// +// 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); +} + +/// \brief Compares std::legendre with analytic expressions +/// for n = 0..6. +/// +/// \tparam ArgContainer a container with a floating-point value_type. +/// \tparam TestFunction a callable that accepts the arguments +/// (unsigned, unsigned, T, T) where T is the value_type of ArgContainer. +/// \param [in] XValues arguments x of std::assoc_legendre(l,m,x) +/// \param [in] ValueSink will be invoked with the arguments of +/// legendre and its expected outcome, i.e., (l, m, x, expected). +template +void compareWithAnalyticExpressions(const ArgContainer& XValues, + TestFunction&& ValueSink) { + std::vector > > TableOfCoefficients; + + // m = 1 and l = 1..4 + TableOfCoefficients.emplace_back(); + TableOfCoefficients.back().push_back({1.0L}); + TableOfCoefficients.back().push_back({3.0L, 0.0L}); + TableOfCoefficients.back().push_back({7.5L, 0.0L, -1.5L}); + TableOfCoefficients.back().push_back({17.5L, 0.0L, -7.5L, 0.0L}); + + // 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({3.0L}); + TableOfCoefficients.back().push_back({15.0L, 0.0L}); + TableOfCoefficients.back().push_back({52.5L, 0.0L, -7.5L}); + + // 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({15.0L}); + TableOfCoefficients.back().push_back({105.0L, 0.0L}); + + // 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({105.0L}); + + for (auto x : XValues) { + for (size_t i = 0; i < TableOfCoefficients.size(); ++i) { + const unsigned m = static_cast(i) + 1; + const long double Factor = std::powl((1.0L - x) * (1.0L + x), 0.5L * m); + for (size_t j = 0; j < TableOfCoefficients[i].size(); ++j) { + const unsigned int l = static_cast(j) + 1; + const long double ExpectedValue = + Factor * evaluatePolynomial(TableOfCoefficients[i][j], + static_cast(x)); + ValueSink(l, m, x, static_cast(ExpectedValue)); + } + } + } +} + +template +void compareWithAnalyticExpressions(const T AbsTolerance, + const T RelTolerance) { + 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); + }; + + const T Samples[] = {T(0), T(0.1), T(0.2), T(0.5), T(0.9), T(1.0)}; + compareWithAnalyticExpressions(Samples, 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 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) { + 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 long double UpperBound = + 1.0L / + std::sqrtl(std::sqrtl(1.0L + n * (n + 1u) * (1.0L - x) * (1.0L + x))); + const T LegendreAbs = std::abs(std::legendre(n, x)); + const T UpperBoundWithTolerance = + static_cast(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,285 @@ +//===----------------------------------------------------------------------===// +// +// 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 = { + {}, + {0.0}, + {0.577350269189626}, + {0.0, 0.774596669241483}, + {0.339981043584856, 0.86113631159405}, + {0.0, 0.538469310105683, 0.90617984593866}, + {0.238619186083197, 0.661209386466265, 0.93246951420315}, + {0.0, 0.405845151377397, 0.741531185599394, 0.949107912342759}, + {0.183434642495650, 0.525532409916329, 0.796666477413627, + 0.960289856497536}, + {0.0, 0.324253423403809, 0.613371432700590, 0.836031107326636, + 0.968160239507626}, + {0.148874338981631, 0.433395394129247, 0.679409568299024, + 0.865063366688985, 0.973906528517172}, + {0.0, 0.269543155952345, 0.51909612920681, 0.730152005574049, + 0.887062599768095, 0.978228658146057}, + {0.125233408511469, 0.367831498998180, 0.587317954286617, + 0.769902674194305, 0.904117256370475, 0.981560634246719}, + {0.0, 0.230458315955135, 0.448492751036447, 0.642349339440340, + 0.801578090733310, 0.917598399222978, 0.984183054718588}, + {0.108054948707344, 0.319112368927890, 0.515248636358154, + 0.687292904811685, 0.827201315069765, 0.928434883663574, + 0.986283808696812}, + {0.0, 0.201194093997435, 0.394151347077563, 0.570972172608539, + 0.724417731360170, 0.848206583410427, 0.937273392400706, + 0.987992518020485}, + {0.095012509837637, 0.281603550779259, 0.458016777657227, + 0.617876244402644, 0.755404408355003, 0.865631202387832, + 0.944575023073233, 0.989400934991650}}; + for (size_t n = 0; n < LegendreZeros.size(); ++n) { + for (double Root : LegendreZeros[n]) { + const T HalfWidthInterval = + T(0.5) * (RelTolerance * T(Root) + AbsTolerance); + const T LeftOfRoot = T(Root) - HalfWidthInterval; + const T RightOfRoot = T(Root) + HalfWidthInterval; + assert(LeftOfRoot != T(Root)); + assert(RightOfRoot != T(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 ArgContainer a container with a floating-point value_type. +/// \tparam TestFunction a callable that accepts the arguments +/// (unsigned, unsigned, T) where T is the value_type of ArgContainer. +/// \param [in] XValues arguments x of std::legendre(n,x) +/// \param [in] ValueSink will be invoked with the arguments of +/// legendre and its expected outcome, i.e., (n, x, expected). +template +void compareWithAnalyticExpressions(const ArgContainer& XValues, + TestFunction&& ValueSink) { + std::vector > VectorOfCoeffients; + VectorOfCoeffients.push_back({1.0L}); + VectorOfCoeffients.push_back({1.0L, 0.0L}); + VectorOfCoeffients.push_back({1.5L, 0.0L, -0.5L}); + VectorOfCoeffients.push_back({2.5L, 0.0L, -1.5L, 0.0L}); + VectorOfCoeffients.push_back({4.375L, 0.0L, -3.75L, 0.0L, 0.375L}); + VectorOfCoeffients.push_back({7.875L, 0.0L, -8.75L, 0.0L, 1.875L, 0.0L}); + VectorOfCoeffients.push_back( + {14.4375L, 0.0L, -19.6875L, 0.0L, 6.5625L, 0.0L, -0.3125L}); + + for (auto x : XValues) { + for (const auto& Coefficients : VectorOfCoeffients) { + const unsigned int n = static_cast(Coefficients.size()) - 1; + const long double ExpectedValue = + ::evaluatePolynomial(Coefficients, static_cast(x)); + ValueSink(n, x, static_cast(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 T Samples[] = {T(0), T(0.1), T(0.2), T(0.5), T(0.9), T(1.0)}; + compareWithAnalyticExpressions(Samples, 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 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) { + 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(); +}