diff --git a/libcxx/include/experimental/__hermite b/libcxx/include/experimental/__hermite new file mode 100644 --- /dev/null +++ b/libcxx/include/experimental/__hermite @@ -0,0 +1,54 @@ +// -*- C++ -*- +//===--------------------------- __hermite --------------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifndef _LIBCPP_HERMITE_H_ +#define _LIBCPP_HERMITE_H_ + +#include +#include +#include + + +/// @return the hermite polynomial \f$ H_{n}(x) \f$ +/// @note The implementation is based on the recurrence formula +/// @f[ +/// nH_{n+1}(x) = 2x H_{n}(x) - 2 n H_{n-1} +/// @f] +/// Press, William H., et al. Numerical recipes 3rd edition: The art of scientific computing. +/// Cambridge university press, 2007, p. 182. +template +_Real +__libcpp_hermite_recurrence(unsigned __n, _Real __x) +{ + if (__n == 0u) + return _Real(1); + + _Real __t2(1); + _Real __t1 = _Real(2)*__x; + for (unsigned __i = 1; __i < __n; ++__i) + { + const _Real __t0 = _Real(2)*(__x*__t1 - _Real(__i)*__t2); + __t2 = __t1; + __t1 = __t0; + } + return __t1; +} + +template +_Real +__libcpp_hermite(unsigned __n, _Real __x) +{ + if (std::isnan(__x)) + return std::numeric_limits<_Real>::quiet_NaN(); + + return __libcpp_hermite_recurrence(__n, __x); +} + + +#endif // _LIBCPP_HERMITE_H_ diff --git a/libcxx/include/experimental/__laguerre b/libcxx/include/experimental/__laguerre new file mode 100644 --- /dev/null +++ b/libcxx/include/experimental/__laguerre @@ -0,0 +1,75 @@ +// -*- C++ -*- +//===--------------------------- __assoc_laguerre --------------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifndef _LIBCPP_LAGUERRE_H_ +#define _LIBCPP_LAGUERRE_H_ + +#include +#include +#include +#include + + +/// @return the generalized laguerre polynomial \f$ L_{n}^{(\alpha)}(x) \f$ +/// @note The implementation is based on the recurrence formula +/// @f[ +/// nL_{n}^{(\alpha)}(x) = (-x + 2n + \alpha - 1) L_{n-1}^{(\alpha)}(x) - +/// (n + \alpha - 1) L_{n-2}^{(\alpha)}(x) +/// @f] +/// Press, William H., et al. Numerical recipes 3rd edition: The art of scientific computing. +/// Cambridge university press, 2007, p. 182. +template +_Real +__libcpp_generalized_laguerre_recurrence(unsigned __n, _Real __alpha, _Real __x) +{ + if (__n == 0u) + return _Real(1); + + _Real __delta = __alpha - __x; + _Real __li = _Real(1) + __delta; + const _Real __alpham1 = __alpha - _Real(1); + for (unsigned __i = 2; __i <= __n; ++__i) + { + __delta = (__delta*(_Real(__i) + __alpham1) - __x*__li)/_Real(__i); + __li += __delta; + } + return __li; +} + + +template +_Real +__libcpp_assoc_laguerre(unsigned __n, unsigned __m, _Real __x) +{ + if (std::isnan(__x)) + return std::numeric_limits<_Real>::quiet_NaN(); + + if (__x < _Real(0)) + _VSTD::__throw_domain_error("Argument of assoc_laguerre function is out of range"); + + + return __libcpp_generalized_laguerre_recurrence(__n, _Real(__m),__x); +} + +template +_Real +__libcpp_laguerre(unsigned __n, _Real __x) +{ + if (std::isnan(__x)) + return std::numeric_limits<_Real>::quiet_NaN(); + + if (__x < _Real(0)) + _VSTD::__throw_domain_error("Argument of laguerre function is out of range"); + + return __libcpp_generalized_laguerre_recurrence(__n, _Real(0),__x); +} + + + +#endif // _LIBCPP_LAGUERRE_H_ diff --git a/libcxx/include/experimental/__legendre b/libcxx/include/experimental/__legendre new file mode 100644 --- /dev/null +++ b/libcxx/include/experimental/__legendre @@ -0,0 +1,136 @@ +// -*- C++ -*- +//===--------------------------- __legendre --------------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifndef _LIBCPP_LEGENDRE_H_ +#define _LIBCPP_LEGENDRE_H_ + +#include +#include +#include +#include + + +/// @return the Legendre polynomial \f$ P_{n}(x) \f$ +/// @note 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] +/// Press, William H., et al. Numerical recipes 3rd edition: The art of scientific computing. +/// Cambridge university press, 2007, p. 182. +template +_Real +__libcpp_legendre_recurrence(unsigned __n, _Real __x) +{ + 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; +} + +template +_Real +__libcpp_legendre(unsigned __n, _Real __x) +{ + if (std::isnan(__x)) + return std::numeric_limits<_Real>::quiet_NaN(); + + if (std::abs(__x) > _Real(1)) + _VSTD::__throw_domain_error("Argument of legendre function is out of range"); + + return __libcpp_legendre_recurrence(__n, __x); +} + + +/// @return \f$ s^{-m} P_{l}^{m}(x) \f$ with an additonal scaling factor to prevent overflow. +/// @note 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] +/// @attention The starting point of the recursion grows exponentially with __m! +/// For large m, 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$ +template +_Real +__libcpp_assoc_legendre_recurrence(unsigned __l, unsigned __m, _Real __x, _Real __scale = _Real(1)) +{ + if (__m == 0u) + return __libcpp_legendre_recurrence(__l, __x); + + // libstdc++ throws a domain error but the return value as specified in + // ISO/IEC JTC 1/SC 22/WG 21 N3060 is actually well-defined for __m > __l + if (__l < __m) + return _Real(0); + + if (__l == 0u) + return _Real(1); + + + _Real __pmm = _Real(1); + // Alternatively, we could have This is the so-called Condon-Shortley phase term, + // which is omitted in the C++17 standard's definition of std::assoc_laguerre. + // This might be confusing because this term is often used in literature, e.g., + // Press, William H., et al. Numerical recipes 3rd edition: The art of scientific computing. + // Cambridge university press, 2007, p. 293. Even the current version of libstdc++ does that. + // msvc conforms to the C++17 standard in this respect. + // _Real __pmm = __m % 2 == 0 ? _Real(1) : _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)*__scale); + 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; +} + +template +_Real +__libcpp_assoc_legendre(unsigned __n, unsigned __m, _Real __x) +{ + if (std::isnan(__x)) + return std::numeric_limits<_Real>::quiet_NaN(); + + if (std::abs(__x) > _Real(1)) + _VSTD::__throw_domain_error("Argument of assoc_legendre function is out of range"); + + return __libcpp_assoc_legendre_recurrence(__n, __m, __x); +} + + +#endif // _LIBCPP_LEGENDRE_H_ diff --git a/libcxx/include/experimental/cmath b/libcxx/include/experimental/cmath new file mode 100644 --- /dev/null +++ b/libcxx/include/experimental/cmath @@ -0,0 +1,136 @@ +// -*- C++ -*- +//===---------------------------- cmath ----------------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifdef _LIBCPP_EXPERIMENTAL_CMATH + +#if defined(__STDCPP_WANT_MATH_SPEC_FUNCS__) && !((__STDCPP_WANT_MATH_SPEC_FUNCS__ + 0) == 0) +#define _LIBCPP_EXPECT_ENABLE_MATH_SPEC_FUNCS +#endif + +#if defined(_LIBCPP_EXPECT_ENABLE_MATH_SPEC_FUNCS) && !defined(_LIBCPP_ENABLE_MATH_SPEC_FUNCS) +#error "Detected inconsistent usage of __STDCPP_WANT_MATH_SPEC_FUNCS__!" \ + " has already been included previously" \ + "with disabled support for special mathematical functions." +#endif + +#if !defined(_LIBCPP_EXPECT_ENABLE_MATH_SPEC_FUNCS) && defined(_LIBCPP_ENABLE_MATH_SPEC_FUNCS) +#error "Detected inconsistent usage of __STDCPP_WANT_MATH_SPEC_FUNCS__!" \ + " has already been included previously" \ + "with enabled support for special mathematical functions." +#endif + +#ifdef _LIBCPP_EXPECT_ENABLE_MATH_SPEC_FUNCS +#undef _LIBCPP_EXPECT_ENABLE_MATH_SPEC_FUNCS +#endif + +#endif // _LIBCPP_EXPERIMENTAL_CMATH + + + +#ifndef _LIBCPP_EXPERIMENTAL_CMATH +#define _LIBCPP_EXPERIMENTAL_CMATH + + +// The following macro name shall be conditionally defined by the implementation +// to indicate conformance to the International Standard ISO/IEC JTC 1/SC 22/WG 21 N3060 +//#define __STDCPP_MATH_SPEC_FUNCS__ 201003L + +#include +#include + +#if defined(__STDCPP_WANT_MATH_SPEC_FUNCS__) && !((__STDCPP_WANT_MATH_SPEC_FUNCS__ + 0) == 0) +#define _LIBCPP_ENABLE_MATH_SPEC_FUNCS +#endif + +#ifdef _LIBCPP_ENABLE_MATH_SPEC_FUNCS +#if _LIBCPP_STD_VER > 14 + +#include +#include +#include + +_LIBCPP_BEGIN_NAMESPACE_EXPERIMENTAL + +inline _LIBCPP_INLINE_VISIBILITY double assoc_laguerre(unsigned __lcpp_n, unsigned __lcpp_m, double __lcpp_x) +{ + return __libcpp_assoc_laguerre(__lcpp_n , __lcpp_m, __lcpp_x); +} +inline _LIBCPP_INLINE_VISIBILITY float assoc_laguerref(unsigned __lcpp_n, unsigned __lcpp_m, float __lcpp_x) +{ + return static_cast(__libcpp_assoc_laguerre(__lcpp_n , __lcpp_m, static_cast(__lcpp_x))); +} +inline _LIBCPP_INLINE_VISIBILITY long double assoc_laguerrel(unsigned __lcpp_n, unsigned __lcpp_m, long double __lcpp_x) +{ + return __libcpp_assoc_laguerre(__lcpp_n , __lcpp_m, __lcpp_x); +} + + +inline _LIBCPP_INLINE_VISIBILITY double assoc_legendre(unsigned __lcpp_n, unsigned __lcpp_m, double __lcpp_x) +{ + return __libcpp_assoc_legendre(__lcpp_n , __lcpp_m, __lcpp_x); +} +inline _LIBCPP_INLINE_VISIBILITY float assoc_legendref(unsigned __lcpp_n, unsigned __lcpp_m, float __lcpp_x) +{ + // use double internally -- float is too prone to overflow! + return static_cast(__libcpp_assoc_legendre(__lcpp_n , __lcpp_m, static_cast(__lcpp_x))); +} +inline _LIBCPP_INLINE_VISIBILITY long double assoc_legendrel(unsigned __lcpp_n, unsigned __lcpp_m, long double __lcpp_x) +{ + return __libcpp_assoc_legendre(__lcpp_n , __lcpp_m, __lcpp_x); +} + + +inline _LIBCPP_INLINE_VISIBILITY double laguerre(unsigned __lcpp_n, double __lcpp_x) +{ + return __libcpp_laguerre(__lcpp_n ,__lcpp_x); +} +inline _LIBCPP_INLINE_VISIBILITY float laguerref(unsigned __lcpp_n, float __lcpp_x) +{ + return static_cast(__libcpp_laguerre(__lcpp_n, static_cast(__lcpp_x))); +} +inline _LIBCPP_INLINE_VISIBILITY long double laguerrel(unsigned __lcpp_n, long double __lcpp_x) +{ + return __libcpp_laguerre(__lcpp_n ,__lcpp_x); +} + + +inline _LIBCPP_INLINE_VISIBILITY double legendre(unsigned __lcpp_n, double __lcpp_x) +{ + return __libcpp_legendre(__lcpp_n ,__lcpp_x); +} +inline _LIBCPP_INLINE_VISIBILITY float legendref(unsigned __lcpp_n, float __lcpp_x) +{ + return static_cast(__libcpp_legendre(__lcpp_n, static_cast(__lcpp_x))); +} +inline _LIBCPP_INLINE_VISIBILITY long double legendrel(unsigned __lcpp_n, long double __lcpp_x) +{ + return __libcpp_legendre(__lcpp_n ,__lcpp_x); +} + +inline _LIBCPP_INLINE_VISIBILITY double hermite(unsigned __lcpp_n, double __lcpp_x) +{ + return __libcpp_hermite(__lcpp_n ,__lcpp_x); +} +inline _LIBCPP_INLINE_VISIBILITY float hermitef(unsigned __lcpp_n, float __lcpp_x) +{ + // use double internally -- float is too prone to overflow! + return static_cast(__libcpp_hermite(__lcpp_n , static_cast(__lcpp_x))); +} +inline _LIBCPP_INLINE_VISIBILITY long double hermitel(unsigned __lcpp_n, long double __lcpp_x) +{ + return __libcpp_hermite(__lcpp_n ,__lcpp_x); +} + + +_LIBCPP_END_NAMESPACE_EXPERIMENTAL + +#endif +#endif // _LIBCPP_ENABLE_MATH_SPEC_FUNCS + +#endif // _LIBCPP_EXPERIMENTAL_CMATH diff --git a/libcxx/test/libcxx/experimental/numerics/c.math/assoc_laguerre.pass.cpp b/libcxx/test/libcxx/experimental/numerics/c.math/assoc_laguerre.pass.cpp new file mode 100644 --- /dev/null +++ b/libcxx/test/libcxx/experimental/numerics/c.math/assoc_laguerre.pass.cpp @@ -0,0 +1,162 @@ +//===----------------------------------------------------------------------===// +// +// 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 +// +//===----------------------------------------------------------------------===// + +// + +#define __STDCPP_WANT_MATH_SPEC_FUNCS__ 1 +#include +#include +#include + +#if _LIBCPP_STD_VER > 14 + +template +void +testAssocLaguerreNaNPropagation() +{ + const unsigned maxN = 127; + const T x = std::numeric_limits::quiet_NaN(); + for (unsigned n = 0; n <= maxN; ++n) + { + for (unsigned m = 0; m <= maxN; ++m) + { + assert(std::isnan(std::experimental::assoc_laguerre(n, m, x))); + } + } +} + +template +void +testAssocLaguerreNotNaN(const T x) +{ + assert(!std::isnan(x)); + const unsigned maxN = 127; + for (unsigned n=0; n<=maxN; ++n) + { + for (unsigned m = 0; m <= maxN; ++m) + { + assert(!std::isnan(std::experimental::assoc_laguerre(n, m, x))); + } + } +} + +template +void +testAssocLaguerreThrows(const T x) +{ +#ifndef _LIBCPP_NO_EXCEPTIONS + const unsigned maxN = 127; + for (unsigned n=0; n<=maxN; ++n) + { + for (unsigned m = 0; m <= maxN; ++m) + { + bool throws = false; + try + { + std::experimental::assoc_laguerre(n, m, x); + } + catch(const std::domain_error&) + { + throws = true; + } + assert(throws); + } + } +#endif // _LIBCPP_NO_EXCEPTIONS +} + +template +void +testAssocLaguerreVsLaguerre(const T x, const T absTolerance, const T relTolerance) +{ + assert(!std::isnan(x)); + const unsigned maxN = 127; + for (unsigned n=0; n<=maxN; ++n) + { + for (unsigned m = 0; m <= maxN; ++m) + { + const T result = std::experimental::assoc_laguerre(n, 0, x); + const T expectedResult = std::experimental::laguerre(n, x); + const T tolerance = absTolerance + std::abs(expectedResult)*relTolerance; + const T difference = std::abs(result-expectedResult); + assert(difference <= tolerance); + } + } +} + +template +void +testAssocLaguerreAnalytic(const T x, const T absTolerance, const T relTolerance) +{ + assert(!std::isnan(x)); + const auto compareFloatingPoint = [absTolerance, relTolerance](const T result, const T expected) + { + if (std::isinf(expected) && std::isinf(result)) + return true; + + if (std::isnan(expected) || std::isnan(result)) + return false; + + const T tolerance = absTolerance + std::abs(expected)*relTolerance; + return std::abs(result-expected) < tolerance; + }; + + const auto l0 = [](T x, T a) { return T(1); }; + const auto l1 = [](T x, T a) { return -x+a+1; }; + const auto l2 = [](T x, T a) { return x*x/T(2)-(a+T(2))*x+(a+T(1))*(a+T(2))/T(2); }; + const auto l3 = [](T x, T a) { + T result = -x*x*x/T(6); + result += (a+T(3))*x*x/T(2); + result += -(a+T(2))*(a+T(3))*x/T(2); + result += (a+T(1))*(a+T(2))*(a+T(3))/T(6); + return result; + }; + + for (unsigned m=0;m<128;++m) + { + const T alpha = static_cast(m); + assert(compareFloatingPoint(std::experimental::assoc_laguerre(0,m,x), l0(x,alpha))); + assert(compareFloatingPoint(std::experimental::assoc_laguerre(1,m,x), l1(x,alpha))); + assert(compareFloatingPoint(std::experimental::assoc_laguerre(2,m,x), l2(x,alpha))); + assert(compareFloatingPoint(std::experimental::assoc_laguerre(3,m,x), l3(x,alpha))); + } + + + +} + +template +void +testAssocLaguerre(const T absTolerance, const T relTolerance) +{ + testAssocLaguerreNaNPropagation(); + testAssocLaguerreThrows(T(-5)); + + const T samples[] = {T(0.0),T(0.1),T(0.5),T(1.0),T(10.0)}; + + for (T x : samples) + { + testAssocLaguerreNotNaN(x); + testAssocLaguerreAnalytic(x, absTolerance, relTolerance); + testAssocLaguerreVsLaguerre(x, absTolerance, relTolerance); + } +} + +#endif + + + +int main(int, char**) +{ +#if _LIBCPP_STD_VER > 14 + testAssocLaguerre(1e-5f, 1e-5f); + testAssocLaguerre(1e-9, 1e-9); + testAssocLaguerre(1e-12, 1e-12); +#endif + return 0; +} diff --git a/libcxx/test/libcxx/experimental/numerics/c.math/assoc_legendre.pass.cpp b/libcxx/test/libcxx/experimental/numerics/c.math/assoc_legendre.pass.cpp new file mode 100644 --- /dev/null +++ b/libcxx/test/libcxx/experimental/numerics/c.math/assoc_legendre.pass.cpp @@ -0,0 +1,213 @@ +//===----------------------------------------------------------------------===// +// +// 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 +// +//===----------------------------------------------------------------------===// + +// + +#define __STDCPP_WANT_MATH_SPEC_FUNCS__ 1 +#include +#include +#include + +#if _LIBCPP_STD_VER > 14 + +template +void +testAssocLegendreNaNPropagation() +{ + const unsigned maxN = 127; + const T x = std::numeric_limits::quiet_NaN(); + for (unsigned n = 0; n <= maxN; ++n) + { + for (unsigned m = 0; m <= maxN; ++m) + { + assert(std::isnan(std::experimental::assoc_legendre(n, m, x))); + } + } +} + +template +void +testAssocLegendreNotNaN(const T x) +{ + assert(!std::isnan(x)); + const unsigned maxN = 127; + for (unsigned n=0; n<=maxN; ++n) + { + for (unsigned m = 0; m <= maxN; ++m) + { + assert(!std::isnan(std::experimental::assoc_legendre(n, m, x))); + } + } +} + + +template +void +testAssocLegendreThrows(const T x) +{ +#ifndef _LIBCPP_NO_EXCEPTIONS + const unsigned maxN = 127; + for (unsigned n=0; n<=maxN; ++n) + { + for (unsigned m = 0; m <= maxN; ++m) + { + bool throws = false; + try + { + std::experimental::assoc_legendre(n, m, x); + } + catch(const std::domain_error&) + { + throws = true; + } + assert(throws); + } + } +#endif // _LIBCPP_NO_EXCEPTIONS +} + + +template +void +testAssocLegendreVsLegendre(const T x, const T absTolerance, const T relTolerance) +{ + assert(!std::isnan(x)); + const unsigned maxN = 127; + for (unsigned n=0; n<=maxN; ++n) + { + for (unsigned m = 0; m <= maxN; ++m) + { + const T result = std::experimental::assoc_legendre(n, 0, x); + const T expectedResult = std::experimental::legendre(n, x); + const T tolerance = absTolerance + std::abs(expectedResult)*relTolerance; + const T difference = std::abs(result-expectedResult); + assert(difference <= tolerance); + } + } +} + +template +void +testAssocLegendreAnalytic(const T x, const T absTolerance, const T relTolerance) +{ + assert(!std::isnan(x)); + const auto compareFloatingPoint = [absTolerance, relTolerance](const T result, const T expectedResult) + { + if (std::isinf(expectedResult) && std::isinf(result)) + return true; + + if (std::isnan(expectedResult) || std::isnan(result)) + return false; + + const T tolerance = absTolerance + std::abs(expectedResult)*relTolerance; + return std::abs(result-expectedResult) < tolerance; + }; + + + const auto l00 = [](T x) { return T(1); }; + + const auto l10 = [](T x) { return x; }; + const auto l11 = [](T x) { return std::sqrt((T(1)-x)*(T(1)+x)); }; + + const auto l20 = [](T x) { return (T(3)*x*x-T(1))/T(2); }; + const auto l21 = [](T x) { return T(3)*x*std::sqrt((T(1)-x)*(T(1)+x)); }; + const auto l22 = [](T x) { return T(3)*(T(1)-x)*(T(1)+x); }; + + const auto l30 = [](T x) { return (T(5)*x*x-T(3))*x/T(2); }; + const auto l31 = [](T x) { return T(3)/T(2)*(T(5)*x*x-T(1))*std::sqrt((T(1)-x)*(T(1)+x)); }; + const auto l32 = [](T x) { return T(15)*x*(T(1)-x)*(T(1)+x); }; + const auto l33 = [](T x) + { + const T temp = (T(1)-x)*(T(1)+x); + return T(15)*temp*std::sqrt(temp); + }; + + const auto l40 = [](T x) { return (T(35)*x*x*x*x-T(30)*x*x + T(3))/T(8); }; + const auto l41 = [](T x) { return T(5)/T(2)*x*(T(7)*x*x-T(3))*std::sqrt((T(1)-x)*(T(1)+x)); }; + const auto l42 = [](T x) { return T(15)/T(2)*(T(7)*x*x-1)*(T(1)-x)*(T(1)+x); }; + const auto l43 = [](T x) + { + const T temp = (T(1)-x)*(T(1)+x); + return T(105)*x*temp*std::sqrt(temp); + }; + const auto l44 = [](T x) + { + const T temp = (T(1)-x)*(T(1)+x); + return T(105)*temp*temp; + }; + + assert(compareFloatingPoint(std::experimental::assoc_legendre(0,0,x), l00(x))); + + assert(compareFloatingPoint(std::experimental::assoc_legendre(1,0,x), l10(x))); + assert(compareFloatingPoint(std::experimental::assoc_legendre(1,1,x), l11(x))); + + assert(compareFloatingPoint(std::experimental::assoc_legendre(2,0,x), l20(x))); + assert(compareFloatingPoint(std::experimental::assoc_legendre(2,1,x), l21(x))); + assert(compareFloatingPoint(std::experimental::assoc_legendre(2,2,x), l22(x))); + + assert(compareFloatingPoint(std::experimental::assoc_legendre(3,0,x), l30(x))); + assert(compareFloatingPoint(std::experimental::assoc_legendre(3,1,x), l31(x))); + assert(compareFloatingPoint(std::experimental::assoc_legendre(3,2,x), l32(x))); + assert(compareFloatingPoint(std::experimental::assoc_legendre(3,3,x), l33(x))); + + assert(compareFloatingPoint(std::experimental::assoc_legendre(4,0,x), l40(x))); + assert(compareFloatingPoint(std::experimental::assoc_legendre(4,1,x), l41(x))); + assert(compareFloatingPoint(std::experimental::assoc_legendre(4,2,x), l42(x))); + assert(compareFloatingPoint(std::experimental::assoc_legendre(4,3,x), l43(x))); + assert(compareFloatingPoint(std::experimental::assoc_legendre(4,4,x), l44(x))); + + try + { + const unsigned maxN = 127; + for (unsigned n=0; n<=maxN; ++n) + { + for (unsigned m = m+1; m <= maxN; ++m) + { + assert(std::experimental::assoc_legendre(4,0,x) <= absTolerance); + } + } + } + catch (const std::domain_error&) + { + // Should not throw! The expression given in + // ISO/IEC JTC 1/SC 22/WG 21 N3060 is actually well-defined for m > n! + assert(false); + } +} + +template +void +testAssocLegendre(const T absTolerance, const T relTolerance) +{ + testAssocLegendreNaNPropagation(); + testAssocLegendreThrows(T(-5)); + testAssocLegendreThrows(T(5)); + + const T samples[] = {T(-1.0),T(-0.5),T(-0.1),T(0.0),T(0.1),T(0.5),T(1.0)}; + + for (T x : samples) + { + testAssocLegendreNotNaN(x); + testAssocLegendreVsLegendre(x, absTolerance, relTolerance); + testAssocLegendreAnalytic(x, absTolerance, relTolerance); + } +} + +#endif + + + +int main(int, char**) +{ +#if _LIBCPP_STD_VER > 14 + testAssocLegendre(1e-6f, 1e-6f); + testAssocLegendre(1e-9, 1e-9); + testAssocLegendre(1e-12, 1e-12); +#endif + return 0; +} diff --git a/libcxx/test/libcxx/experimental/numerics/c.math/hermite.pass.cpp b/libcxx/test/libcxx/experimental/numerics/c.math/hermite.pass.cpp new file mode 100644 --- /dev/null +++ b/libcxx/test/libcxx/experimental/numerics/c.math/hermite.pass.cpp @@ -0,0 +1,100 @@ +//===----------------------------------------------------------------------===// +// +// 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 +// +//===----------------------------------------------------------------------===// + +// + +#define __STDCPP_WANT_MATH_SPEC_FUNCS__ 1 +#include +#include +#include + +#if _LIBCPP_STD_VER > 14 + +template +void +testHermiteNaNPropagation() +{ + const unsigned maxN = 127; + const T x = std::numeric_limits::quiet_NaN(); + for (unsigned n=0; n<=maxN; ++n) + { + assert(std::isnan(std::experimental::hermite(n,x))); + } +} + +template +void +testHermiteNotNaN(const T x) +{ + assert(!std::isnan(x)); + const unsigned maxN = 127; + for (unsigned n=0; n<=maxN; ++n) + { + assert(!std::isnan(std::experimental::hermite(n,x))); + } +} + +template +void +testHermiteAnalytic(const T x, const T absTolerance, const T relTolerance) +{ + assert(!std::isnan(x)); + const auto compareFloatingPoint = [absTolerance, relTolerance](const T result, const T expected) + { + if (std::isinf(expected) && std::isinf(result)) + return true; + + if (std::isnan(expected) || std::isnan(result)) + return false; + + const T tolerance = absTolerance + std::abs(expected)*relTolerance; + return std::abs(result-expected) < tolerance; + }; + + const auto h0 = [](T) { return T(1); }; + const auto h1 = [](T x) { return T(2)*x; }; + const auto h2 = [](T x) { return T(4)*x*x - T(2);}; + const auto h3 = [](T x) { return x*(T(8)*x*x - T(12)); }; + const auto h4 = [](T x) { return (T(16)*x*x*x*x - T(48)*x*x + T(12)); }; + const auto h5 = [](T x) { return x*(T(32)*x*x*x*x - T(160)*x*x + T(120)); }; + + assert(compareFloatingPoint(std::experimental::hermite(0,x), h0(x))); + assert(compareFloatingPoint(std::experimental::hermite(1,x), h1(x))); + assert(compareFloatingPoint(std::experimental::hermite(2,x), h2(x))); + assert(compareFloatingPoint(std::experimental::hermite(3,x), h3(x))); + assert(compareFloatingPoint(std::experimental::hermite(4,x), h4(x))); + assert(compareFloatingPoint(std::experimental::hermite(5,x), h5(x))); +} + +template +void +testHermite(const T absTolerance, const T relTolerance) +{ + testHermiteNaNPropagation(); + const T samples[] = {T(-1.0),T(-0.5),T(-0.1),T(0.0),T(0.1),T(0.5),T(1.0)}; + + for (T x : samples) + { + testHermiteNotNaN(x); + testHermiteAnalytic(x, absTolerance, relTolerance); + } +} + +#endif + + + +int main(int, char**) +{ +#if _LIBCPP_STD_VER > 14 + testHermite(1e-6f, 1e-6f); + testHermite(1e-9, 1e-9); + testHermite(1e-12, 1e-12); +#endif + return 0; +} diff --git a/libcxx/test/libcxx/experimental/numerics/c.math/laguerre.pass.cpp b/libcxx/test/libcxx/experimental/numerics/c.math/laguerre.pass.cpp new file mode 100644 --- /dev/null +++ b/libcxx/test/libcxx/experimental/numerics/c.math/laguerre.pass.cpp @@ -0,0 +1,130 @@ +//===----------------------------------------------------------------------===// +// +// 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 +// +//===----------------------------------------------------------------------===// + +// + +#define __STDCPP_WANT_MATH_SPEC_FUNCS__ 1 +#include +#include +#include + +#if _LIBCPP_STD_VER > 14 + +template +void +testLaguerreNaNPropagation() +{ + const unsigned maxN = 127; + const T x = std::numeric_limits::quiet_NaN(); + for (unsigned n=0; n<=maxN; ++n) + { + assert(std::isnan(std::experimental::laguerre(n,x))); + } +} + +template +void +testLaguerreNotNaN(const T x) +{ + assert(!std::isnan(x)); + const unsigned maxN = 127; + for (unsigned n=0; n<=maxN; ++n) + { + assert(!std::isnan(std::experimental::laguerre(n,x))); + } +} + +template +void +testLaguerreThrows(const T x) +{ +#ifndef _LIBCPP_NO_EXCEPTIONS + const unsigned maxN = 127; + for (unsigned n=0; n<=maxN; ++n) + { + bool throws = false; + try + { + std::experimental::laguerre(n,x); + } + catch(const std::domain_error&) + { + throws = true; + } + assert(throws); + } +#endif // _LIBCPP_NO_EXCEPTIONS +} + +template +void +testLaguerreAnalytic(const T x, const T absTolerance, const T relTolerance) +{ + assert(!std::isnan(x)); + const auto compareFloatingPoint = [absTolerance, relTolerance](const T result, const T expected) + { + if (std::isinf(expected) && std::isinf(result)) + return true; + + if (std::isnan(expected) || std::isnan(result)) + return false; + + const T tolerance = absTolerance + std::abs(expected)*relTolerance; + return std::abs(result-expected) < tolerance; + }; + + + const auto l0 = [](T) { return T(1); }; + const auto l1 = [](T x) { return -x+1; }; + const auto l2 = [](T x) { return (x*x-T(4)*x+T(2))/T(2); }; + const auto l3 = [](T x) { return (-x*x*x+T(9)*x*x-T(18)*x+T(6))/T(6); }; + const auto l4 = [](T x) { return (x*x*x*x-T(16)*x*x*x+T(72)*x*x-T(96)*x+T(24))/T(24); }; + const auto l5 = [](T x) { return (-x*x*x*x*x+T(25)*x*x*x*x-T(200)*x*x*x+T(600)*x*x-T(600)*x+T(120))/T(120); }; + const auto l6 = [](T x) + { + return (x*x*x*x*x*x-T(36)*x*x*x*x*x+T(450)*x*x*x*x-T(2400)*x*x*x+T(5400)*x*x-T(4320)*x+T(720))/T(720); + }; + + assert(compareFloatingPoint(std::experimental::laguerre(0,x), l0(x))); + assert(compareFloatingPoint(std::experimental::laguerre(1,x), l1(x))); + assert(compareFloatingPoint(std::experimental::laguerre(2,x), l2(x))); + assert(compareFloatingPoint(std::experimental::laguerre(3,x), l3(x))); + assert(compareFloatingPoint(std::experimental::laguerre(4,x), l4(x))); + assert(compareFloatingPoint(std::experimental::laguerre(5,x), l5(x))); + assert(compareFloatingPoint(std::experimental::laguerre(6,x), l6(x))); +} + +template +void +testLaguerre(const T absTolerance, const T relTolerance) +{ + testLaguerreNaNPropagation(); + testLaguerreThrows(T(-5)); + + const T samples[] = {T(0.0),T(0.1),T(0.5),T(1.0),T(10.0)}; + + for (T x : samples) + { + testLaguerreNotNaN(x); + testLaguerreAnalytic(x, absTolerance, relTolerance); + } +} + +#endif + + + +int main(int, char**) +{ +#if _LIBCPP_STD_VER > 14 + testLaguerre(1e-6f, 1e-6f); + testLaguerre(1e-9, 1e-9); + testLaguerre(1e-12, 1e-12); +#endif + return 0; +} diff --git a/libcxx/test/libcxx/experimental/numerics/c.math/legendre.pass.cpp b/libcxx/test/libcxx/experimental/numerics/c.math/legendre.pass.cpp new file mode 100644 --- /dev/null +++ b/libcxx/test/libcxx/experimental/numerics/c.math/legendre.pass.cpp @@ -0,0 +1,132 @@ +//===----------------------------------------------------------------------===// +// +// 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 +// +//===----------------------------------------------------------------------===// + +// + +#define __STDCPP_WANT_MATH_SPEC_FUNCS__ 1 +#include +#include +#include + +#if _LIBCPP_STD_VER > 14 + +template +void +testLegendreNaNPropagation() +{ + const unsigned maxN = 127; + const T x = std::numeric_limits::quiet_NaN(); + for (unsigned n=0; n<=maxN; ++n) + { + assert(std::isnan(std::experimental::legendre(n,x))); + } +} + +template +void +testLegendreNotNaN(const T x) +{ + assert(!std::isnan(x)); + const unsigned maxN = 127; + for (unsigned n=0; n<=maxN; ++n) + { + assert(!std::isnan(std::experimental::legendre(n,x))); + } +} + +template +void +testLegendreThrows(const T x) +{ +#ifndef _LIBCPP_NO_EXCEPTIONS + const unsigned maxN = 127; + for (unsigned n=0; n<=maxN; ++n) + { + bool throws = false; + try + { + std::experimental::legendre(n,x); + } + catch(const std::domain_error&) + { + throws = true; + } + assert(throws); + } +#endif // _LIBCPP_NO_EXCEPTIONS +} + +template +void +testLegendreAnalytic(const T x, const T absTolerance, const T relTolerance) +{ + assert(!std::isnan(x)); + const auto compareFloatingPoint = [absTolerance, relTolerance](const T result, const T expected) + { + if (std::isinf(expected) && std::isinf(result)) + return true; + + if (std::isnan(expected) || std::isnan(result)) + return false; + + const T tolerance = absTolerance + std::abs(expected)*relTolerance; + return std::abs(result-expected) < tolerance; + }; + + + const auto l0 = [](T) { return T(1); }; + const auto l1 = [](T x) { return x; }; + const auto l2 = [](T x) { return (T(3)*x*x-T(1))/T(2); }; + const auto l3 = [](T x) { return (T(5)*x*x-T(3))*x/T(2); }; + const auto l4 = [](T x) { return (T(35)*x*x*x*x-T(30)*x*x+T(3))/T(8); }; + const auto l5 = [](T x) { return (T(63)*x*x*x*x-T(70)*x*x+T(15))*x/T(8); }; + const auto l6 = [](T x) + { + const T x2 = x*x; + return (T(231)*x2*x2*x2-T(315)*x2*x2+T(105)*x2 - T(5))/T(16); + }; + + assert(compareFloatingPoint(std::experimental::legendre(0,x), l0(x))); + assert(compareFloatingPoint(std::experimental::legendre(1,x), l1(x))); + assert(compareFloatingPoint(std::experimental::legendre(2,x), l2(x))); + assert(compareFloatingPoint(std::experimental::legendre(3,x), l3(x))); + assert(compareFloatingPoint(std::experimental::legendre(4,x), l4(x))); + assert(compareFloatingPoint(std::experimental::legendre(5,x), l5(x))); + assert(compareFloatingPoint(std::experimental::legendre(6,x), l6(x))); +} + +template +void +testLegendre(const T absTolerance, const T relTolerance) +{ + testLegendreNaNPropagation(); + testLegendreThrows(T(-5)); + testLegendreThrows(T(5)); + + const T samples[] = {T(-1.0),T(-0.5),T(-0.1),T(0.0),T(0.1),T(0.5),T(1.0)}; + + for (T x : samples) + { + testLegendreNotNaN(x); + testLegendreAnalytic(x, absTolerance, relTolerance); + } +} + +#endif + + + +int main(int, char**) +{ +#if _LIBCPP_STD_VER > 14 + testLegendre(1e-6f, 1e-6f); + testLegendre(1e-9, 1e-9); + testLegendre(1e-12, 1e-12); +#endif + return 0; +}