Index: libcxx/include/cmath =================================================================== --- libcxx/include/cmath +++ 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 #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER) #pragma GCC system_header @@ -606,6 +617,220 @@ return isfinite(__lcpp_x); } +_LIBCPP_FUNC_VIS +void __libcpp_cmath_domain_error() _NOEXCEPT; + + +#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_{l}(x) \f$. +/// +/// The implementation is based on the recurrence formula +/// \f[ +/// (l+1)P_{l+1}(x) = (2l+1)xP_{l}(x) - lP_{l-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] __l corresponds to the degree \f$ l \f$. +/// \param [in] __x corresponds to the argument \f$ x \f$. +/// \return the Legendre polynomial \f$ P_{l}(x) \f$ +template +inline _LIBCPP_INLINE_VISIBILITY _Real +__libcpp_legendre_recurrence(unsigned __l, _Real __x) _NOEXCEPT { + if (__l == 0u) + return _Real(1); + + _Real __t2(1); + _Real __t1 = __x; + for (unsigned __i = 1; __i < __l; ++__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_{l}(x) \f$. +/// +/// \tparam T one of float, double, or long double +/// \param [in] __l corresponds to the degree \f$ l \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(unsigned __l, _Real __x) { + if (_VSTD::isnan(__x)) + return _VSTD::numeric_limits<_Real>::quiet_NaN(); + + if (_VSTD::abs(__x) > _Real(1)) { + _VSTD::__libcpp_cmath_domain_error(); + return _VSTD::numeric_limits<_Real>::quiet_NaN(); + } + + return __libcpp_legendre_recurrence(__l, __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) + // See p. 38 of Goldberg, D. (1991). + // What every computer scientist should know about floating-point arithmetic. + // ACM Computing Surveys (CSUR), 23(1), 5-48. + const _Real __t = _VSTD::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$ +template +inline _LIBCPP_INLINE_VISIBILITY _Real +__libcpp_assoc_legendre(unsigned __l, unsigned __m, _Real __x) { + if (_VSTD::isnan(__x)) + return _VSTD::numeric_limits<_Real>::quiet_NaN(); + + if (_VSTD::abs(__x) > _Real(1)) { + _VSTD::__libcpp_cmath_domain_error(); + return _VSTD::numeric_limits<_Real>::quiet_NaN(); + } + + 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 + #if _LIBCPP_STD_VER > 17 template constexpr Index: libcxx/src/CMakeLists.txt =================================================================== --- libcxx/src/CMakeLists.txt +++ libcxx/src/CMakeLists.txt @@ -7,6 +7,7 @@ bind.cpp charconv.cpp chrono.cpp + cmath.cpp condition_variable.cpp debug.cpp exception.cpp Index: libcxx/src/cmath.cpp =================================================================== --- /dev/null +++ libcxx/src/cmath.cpp @@ -0,0 +1,36 @@ +//===------------------------------ cmath.cpp -----------------------------===// +// +// 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 +// +//===----------------------------------------------------------------------===// + +#include "cmath" +#include "cerrno" +#include "cfenv" + +_LIBCPP_BEGIN_NAMESPACE_STD + +/// It is unspecified whether math_errhandling is a macro or an +/// identifier with external linkage (see J.1 of the C11 standard). +inline _LIBCPP_INLINE_VISIBILITY +bool __libcpp_cmath_use_errno() _NOEXCEPT { + return (math_errhandling & MATH_ERRNO) != 0; +} + +inline _LIBCPP_INLINE_VISIBILITY +bool __libcpp_cmath_use_floating_point_exceptions() _NOEXCEPT { + return (math_errhandling & MATH_ERREXCEPT) != 0; +} + +/// \brief Signal errors from mathematics functions according to +/// section 7.12.1 of the C11 standard. +void __libcpp_cmath_domain_error() _NOEXCEPT { + if (__libcpp_cmath_use_floating_point_exceptions()) + _VSTD::feraiseexcept(FE_INVALID); + if (__libcpp_cmath_use_errno()) + errno = EDOM; +} + +_LIBCPP_END_NAMESPACE_STD Index: libcxx/test/std/numerics/c.math/c.math.legendre/assoc_legendre_basic.pass.cpp =================================================================== --- /dev/null +++ libcxx/test/std/numerics/c.math/c.math.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); +} Index: libcxx/test/std/numerics/c.math/c.math.legendre/assoc_legendre_comparisons.pass.cpp =================================================================== --- /dev/null +++ libcxx/test/std/numerics/c.math/c.math.legendre/assoc_legendre_comparisons.pass.cpp @@ -0,0 +1,244 @@ +//===----------------------------------------------------------------------===// +// +// 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 = static_cast(xy.second); + if (std::isinf(ExpectedValue)) + continue; + const T Tolerance = std::abs(ExpectedValue) * RelTolerance + AbsTolerance; + const T Error = std::abs( + std::assoc_legendre(l, m, static_cast(xy.first)) - ExpectedValue); + assert(Error <= Tolerance); + } + } +} + +template +void compareWithBoost(T RelTolerance, T AbsTolerance) { + // Generated with gcc 8.3.1 and boost 1.66.0 + // using boost::multiprecision::cpp_bin_float_oct + std::map, std::map > Table; + Table[{10, 10}] = {{0.000000000e+00f, 6.547290750000000000000e+08L}, + {1.000000015e-01f, 6.226408347680519008579e+08L}, + {2.000000030e-01f, 5.338482087653367613489e+08L}, + {3.000000119e-01f, 4.085719730748052363114e+08L}, + {4.000000060e-01f, 2.738155104297545025258e+08L}, + {5.000000000e-01f, 1.553702785400390625000e+08L}, + {6.000000238e-01f, 7.030098340813657269702e+07L}, + {6.999999881e-01f, 2.258981004004501399553e+07L}, + {8.000000119e-01f, 3.958895299377443170423e+06L}, + {8.999999762e-01f, 1.621175838753980941588e+05L}}; + Table[{20, 10}] = {{0.000000000e+00f, -1.612052956674316406250e+12L}, + {1.000000015e-01f, 3.560417986699956011359e+11L}, + {2.000000030e-01f, 1.466609229293489889458e+12L}, + {3.000000119e-01f, -1.092420733564380088043e+12L}, + {4.000000060e-01f, -8.999733248634433357581e+11L}, + {5.000000000e-01f, 1.744502295946207410190e+12L}, + {6.000000238e-01f, -4.244728166109593568832e+11L}, + {6.999999881e-01f, -1.548364851025354883254e+12L}, + {8.000000119e-01f, 2.392758198486276572460e+12L}, + {8.999999762e-01f, 9.900181633268510553219e+11L}}; + Table[{20, 20}] = {{0.000000000e+00f, 3.198309867728777708156e+23L}, + {1.000000015e-01f, 2.892494105990308964133e+23L}, + {2.000000030e-01f, 2.126340753675630014195e+23L}, + {3.000000119e-01f, 1.245473315336152979943e+23L}, + {4.000000060e-01f, 5.593882940857043104312e+22L}, + {5.000000000e-01f, 1.801080697817960690393e+22L}, + {6.000000238e-01f, 3.687398576504732856657e+21L}, + {6.999999881e-01f, 3.807346833982522757930e+20L}, + {8.000000119e-01f, 1.169352142138083731714e+19L}, + {8.999999762e-01f, 1.960909400307456079701e+16L}}; + Table[{50, 10}] = {{0.000000000e+00f, 1.145223824924634340402e+16L}, + {1.000000015e-01f, 2.794368988457766056445e+15L}, + {2.000000030e-01f, -9.926438554299240031969e+15L}, + {3.000000119e-01f, -9.451388348844008469839e+15L}, + {4.000000060e-01f, 8.891978701142798866308e+14L}, + {5.000000000e-01f, 9.180035388465756390148e+15L}, + {6.000000238e-01f, 1.217668478629230978049e+16L}, + {6.999999881e-01f, 1.215354000371287394087e+16L}, + {8.000000119e-01f, 9.322971451346198071300e+14L}, + {8.999999762e-01f, -8.994700004228369487463e+15L}}; + Table[{50, 20}] = {{0.000000000e+00f, -7.873593569733218615651e+32L}, + {1.000000015e-01f, 5.415130362958155159222e+31L}, + {2.000000030e-01f, 7.930926719080271686628e+32L}, + {3.000000119e-01f, -4.005360288314231186108e+31L}, + {4.000000060e-01f, -8.232432232229365316108e+32L}, + {5.000000000e-01f, -4.016986081427444525172e+32L}, + {6.000000238e-01f, 4.086002375842282237800e+32L}, + {6.999999881e-01f, 8.602255288380065959535e+32L}, + {8.000000119e-01f, 1.062267415526937607175e+33L}, + {8.999999762e-01f, 1.663092579695240408210e+33L}}; + Table[{50, 50}] = {{0.000000000e+00f, 2.725392139750729502981e+78L}, + {1.000000015e-01f, 2.119868203082874867369e+78L}, + {2.000000030e-01f, 9.822223488246167325621e+77L}, + {3.000000119e-01f, 2.579073515630015723133e+77L}, + {4.000000060e-01f, 3.486671958530380886198e+76L}, + {5.000000000e-01f, 2.050976025703723881955e+75L}, + {6.000000238e-01f, 3.889805295979395319202e+73L}, + {6.999999881e-01f, 1.332550746119903553704e+71L}, + {8.000000119e-01f, 2.202880522487869978317e+67L}, + {8.999999762e-01f, 2.536713821604446050049e+60L}}; + Table[{100, 10}] = {{0.000000000e+00f, -8.249459718167038283208e+18L}, + {1.000000015e-01f, 6.862784831827366756776e+18L}, + {2.000000030e-01f, -2.346878604334024285454e+18L}, + {3.000000119e-01f, -4.908944436938120307561e+18L}, + {4.000000060e-01f, 8.237081820519236786637e+18L}, + {5.000000000e-01f, 4.238712302196339607619e+18L}, + {6.000000238e-01f, -9.593338266778539462151e+17L}, + {6.999999881e-01f, 4.435490321308495360216e+18L}, + {8.000000119e-01f, 1.577846540648397376697e+18L}, + {8.999999762e-01f, 2.565015383247689253549e+17L}}; + Table[{100, 20}] = {{0.000000000e+00f, 7.772743983615954054163e+38L}, + {1.000000015e-01f, -7.050326156048125707974e+38L}, + {2.000000030e-01f, 4.396623183489620262706e+38L}, + {3.000000119e-01f, 1.147459374439165019899e+38L}, + {4.000000060e-01f, -7.612060446652667077043e+38L}, + {5.000000000e-01f, 3.075049703999956744883e+38L}, + {6.000000238e-01f, 8.267992424776928215129e+38L}, + {6.999999881e-01f, 7.902231488613681556252e+38L}, + {8.000000119e-01f, -8.396400636583794923594e+38L}, + {8.999999762e-01f, 7.161037327097612617070e+37L}}; + Table[{100, 50}] = {{0.000000000e+00f, -1.171214503157837849026e+98L}, + {1.000000015e-01f, 9.011934838994649829986e+97L}, + {2.000000030e-01f, -2.779117999991753228251e+97L}, + {3.000000119e-01f, -3.352047766735719610538e+97L}, + {4.000000060e-01f, 7.143132937401246198217e+97L}, + {5.000000000e-01f, -8.225262033972737473325e+97L}, + {6.000000238e-01f, 5.936375526121234287759e+97L}, + {6.999999881e-01f, 1.808915146385161650542e+97L}, + {8.000000119e-01f, -1.102981549212707410800e+98L}, + {8.999999762e-01f, 9.323157549433173993957e+96L}}; + Table[{100, 100}] = {{0.000000000e+00f, 6.666308670072953744411e+186L}, + {1.000000015e-01f, 4.033157130099916855894e+186L}, + {2.000000030e-01f, 8.658587377156835360803e+185L}, + {3.000000119e-01f, 5.969732406590800712476e+184L}, + {4.000000060e-01f, 1.091062423367808056023e+183L}, + {5.000000000e-01f, 3.775274968288970794232e+180L}, + {6.000000238e-01f, 1.357948023775031303522e+177L}, + {6.999999881e-01f, 1.593657292608815830992e+172L}, + {8.000000119e-01f, 4.355212066771058773438e+164L}, + {8.999999762e-01f, 5.775244447989821167677e+150L}}; + compareWithTable(Table, RelTolerance, AbsTolerance); +} + +int main(int, char**) { + compareWithAnalyticExpressions(1e-6f, 1e-6f); + compareWithAnalyticExpressions(1e-12, 1e-12); + + compareWithBoost(1e-6f, 0.0f); + compareWithBoost(1e-12, 0.0); +} Index: libcxx/test/std/numerics/c.math/c.math.legendre/assoc_legendre_exceptions.pass.cpp =================================================================== --- /dev/null +++ libcxx/test/std/numerics/c.math/c.math.legendre/assoc_legendre_exceptions.pass.cpp @@ -0,0 +1,53 @@ +//===----------------------------------------------------------------------===// +// +// 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 Tests if std::assoc_legendre(l, m, x) raises a domain error via errno +/// if x is not in the range [-1,1] and not NaN. +/// +/// \tparam T one of float, double, or long double +template +void testErrnoOnInvalidArgument() { + 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) { + errno = 0; + const T Value = std::assoc_legendre(l, m, x); + assert(errno == EDOM); + // The values returned by the mathematics functions on domain errors + // are implementation-defined according to J.3.12 of the C11 standard. + // We are going for NaN here. + assert(std::isnan(Value)); + } + } + } +} + +bool doMathematicsFunctionSetErrno() { + return (math_errhandling & MATH_ERRNO) != 0; +} + +int main(int, char**) { + if (doMathematicsFunctionSetErrno()) { + testErrnoOnInvalidArgument(); + testErrnoOnInvalidArgument(); + testErrnoOnInvalidArgument(); + } +} Index: libcxx/test/std/numerics/c.math/c.math.legendre/legendre_basic.pass.cpp =================================================================== --- /dev/null +++ libcxx/test/std/numerics/c.math/c.math.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(); +} Index: libcxx/test/std/numerics/c.math/c.math.legendre/legendre_comparisons.pass.cpp =================================================================== --- /dev/null +++ libcxx/test/std/numerics/c.math/c.math.legendre/legendre_comparisons.pass.cpp @@ -0,0 +1,232 @@ +//===----------------------------------------------------------------------===// +// +// 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, static_cast(xy.first)); + const T ExpectedValue = static_cast(xy.second); + const T Tolerance = std::abs(ExpectedValue) * RelTolerance + AbsTolerance; + const T Error = std::abs(Value - ExpectedValue); + assert(Error <= Tolerance); + } + } +} + +template +void compareWithBoost(T RelTolerance, T AbsTolerance) { + // Generated with gcc 8.3.1 and boost 1.66.0 + // using boost::multiprecision::cpp_bin_float_oct + std::map > Table; + Table[10u] = {{0.000000000e+00f, -2.460937500000000000000e-01L}, + {1.000000015e-01f, -1.221249940215520335990e-01L}, + {2.000000030e-01f, 1.290720324341292720607e-01L}, + {3.000000119e-01f, 2.514763479777521008693e-01L}, + {4.000000060e-01f, 9.683904825713287901691e-02L}, + {5.000000000e-01f, -1.882286071777343750000e-01L}, + {6.000000238e-01f, -2.436627083118402910094e-01L}, + {6.999999881e-01f, 8.580574597893920865941e-02L}, + {8.000000119e-01f, 3.005297781926033732426e-01L}, + {8.999999762e-01f, -2.631454513650919607382e-01L}}; + Table[20u] = {{0.000000000e+00f, 1.761970520019531250000e-01L}, + {1.000000015e-01f, -8.207713575422256062399e-02L}, + {2.000000030e-01f, -9.804218510815285071285e-02L}, + {3.000000119e-01f, 1.802871614549846022334e-01L}, + {4.000000060e-01f, -1.015926361969871907174e-01L}, + {5.000000000e-01f, -4.835838106737355701625e-02L}, + {6.000000238e-01f, 1.591674599490453944508e-01L}, + {6.999999881e-01f, -2.045739567029130982825e-01L}, + {8.000000119e-01f, 2.242045928729420716539e-01L}, + {8.999999762e-01f, -1.493084750921506729662e-01L}}; + Table[50u] = {{0.000000000e+00f, -1.122751726592170484764e-01L}, + {1.000000015e-01f, -3.820582067191225865450e-02L}, + {2.000000030e-01f, 8.343226042617237258518e-02L}, + {3.000000119e-01f, 1.091105387969200243060e-01L}, + {4.000000060e-01f, 4.156906945843224360418e-02L}, + {5.000000000e-01f, -3.105909923960982286895e-02L}, + {6.000000238e-01f, -5.886063262789513678078e-02L}, + {6.999999881e-01f, -1.457284285628031113290e-02L}, + {8.000000119e-01f, 1.387974171832264305537e-01L}, + {8.999999762e-01f, -1.700376526005852748269e-01L}}; + Table[100u] = {{0.000000000e+00f, 7.958923738717876149813e-02L}, + {1.000000015e-01f, -6.389509124671504719492e-02L}, + {2.000000030e-01f, 1.468181119397740238999e-02L}, + {3.000000119e-01f, 5.712746529625854606684e-02L}, + {4.000000060e-01f, -7.225817535745479257209e-02L}, + {5.000000000e-01f, -6.051802596186118687465e-02L}, + {6.000000238e-01f, -2.374728102587777222763e-02L}, + {6.999999881e-01f, -7.713241590611060241231e-02L}, + {8.000000119e-01f, 5.086134685690222334621e-02L}, + {8.999999762e-01f, 1.022654639320334207095e-01L}}; + + compareWithTable(Table, 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); + + compareWithBoost(1e-6f, 0.0f); + compareWithBoost(1e-12, 0.0); +} Index: libcxx/test/std/numerics/c.math/c.math.legendre/legendre_exceptions.pass.cpp =================================================================== --- /dev/null +++ libcxx/test/std/numerics/c.math/c.math.legendre/legendre_exceptions.pass.cpp @@ -0,0 +1,51 @@ +//===----------------------------------------------------------------------===// +// +// 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 std::legendre(l, x) raises a domain error via errno +/// if x is not in the range [-1,1] and not NaN. +/// +/// \tparam T one of float, double, or long double +template +void testErrnoOnInvalidArgument() { + const T InvalidArguments[] = {T(-5), T(2)}; + for (unsigned n = 0; n < 128; ++n) { + for (T x : InvalidArguments) { + errno = 0; + const T Value = std::legendre(n, x); + assert(errno == EDOM); + // The values returned by the mathematics functions on domain errors + // are implementation-defined according to J.3.12 of the C11 standard. + // We are going for NaN here. + assert(std::isnan(Value)); + } + } +} + +bool doMathematicsFunctionSetErrno() { + return (math_errhandling & MATH_ERRNO) != 0; +} + +int main(int, char**) { + if (doMathematicsFunctionSetErrno()) { + testErrnoOnInvalidArgument(); + testErrnoOnInvalidArgument(); + testErrnoOnInvalidArgument(); + } +} Index: libcxx/test/std/numerics/c.math/cmath.pass.cpp =================================================================== --- libcxx/test/std/numerics/c.math/cmath.pass.cpp +++ 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 + template ()))> std::true_type has_abs_imp(int); template @@ -1553,6 +1558,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(); @@ -1626,5 +1672,10 @@ test_tgamma(); test_trunc(); +#if _LIBCPP_STD_VER > 14 + test_assoc_legendre(); + test_legendre(); +#endif + return 0; }