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,50 @@ +//===------------------------ __hermite -------------------------*- C++ -*-===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// +/// +/// \file +/// This file contains the internal implementations of std::hermite. +/// +//===----------------------------------------------------------------------===// + +#ifndef _LIBCPP_EXPERIMENTAL___HERMITE +#define _LIBCPP_EXPERIMENTAL___HERMITE + +#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_EXPERIMENTAL___HERMITE 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,71 @@ +//===------------------------ __laguerre ------------------------*- C++ -*-===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// +/// +/// \file +/// This file contains the internal implementations of std::laguerre +/// and std::assoc_laguerre. +/// +//===----------------------------------------------------------------------===// + + +#ifndef _LIBCPP_EXPERIMENTAL___LAGUERRE +#define _LIBCPP_EXPERIMENTAL___LAGUERRE + +#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_EXPERIMENTAL___LAGUERRE 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,123 @@ +//===------------------------ __legendre ------------------------*- C++ -*-===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// +/// +/// \file +/// This file contains the internal implementations of std::legendre +/// and std::assoc_legendre. +/// +//===----------------------------------------------------------------------===// + + +#ifndef _LIBCPP_EXPERIMENTAL___LEGENDRE +#define _LIBCPP_EXPERIMENTAL___LEGENDRE + +#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$ +/// \attention The so-called Condon-Shortley phase term is omitted in the C++17 +/// standard's definition of std::assoc_laguerre. +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); + + 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) * __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_EXPERIMENTAL___LEGENDRE 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,103 @@ +// -*- C++ -*- +//===---------------------------- cmath ----------------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifndef _LIBCPP_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 _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 // _LIBCPP_STD_VER > 14 + +#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,134 @@ +//===----------------------------------------------------------------------===// +// +// 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 +#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 ExpectedResult) { + if (std::isinf(ExpectedResult) && std::isinf(Result)) + return true; + + if (std::isnan(ExpectedResult) || std::isnan(Result)) + return false; + + const T Tolerance = + AbsTolerance + std::abs(ExpectedResult) * RelTolerance; + return std::abs(Result - ExpectedResult) < Tolerance; + }; + + const auto l0 = [](T, unsigned) { return T(1); }; + const auto l1 = [](T x, unsigned m) { return -x + T(m + 1); }; + const auto l2 = [](T x, unsigned m) { + return x * x / T(2) - T(m + 2) * x + T(m + 1) * T(m + 2) / T(2); + }; + const auto l3 = [](T x, unsigned m) { + return -x * x * x / T(6) + T(m + 3) * x * x / T(2) - + T(m + 2) * T(m + 3) * x / T(2) + + T(m + 1) * T(m + 2) * T(m + 3) / T(6); + }; + + for (unsigned m = 0; m < 128; ++m) { + assert(compareFloatingPoint(std::experimental::assoc_laguerre(0, m, x), + l0(x, m))); + assert(compareFloatingPoint(std::experimental::assoc_laguerre(1, m, x), + l1(x, m))); + assert(compareFloatingPoint(std::experimental::assoc_laguerre(2, m, x), + l2(x, m))); + assert(compareFloatingPoint(std::experimental::assoc_laguerre(3, m, x), + l3(x, m))); + } +} + +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,203 @@ +//===----------------------------------------------------------------------===// +// +// 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 +#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) { 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 = n + 1; m <= MaxN; ++m) { + assert(std::experimental::assoc_legendre(n, m, 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,290 @@ +//===----------------------------------------------------------------------===// +// +// 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 +#include +#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))); +} + +/// \details This method checks if the following recurrence relation holds: +/// \f[ +/// H_{n+1}(x) = 2x H_{n}(x) - 2n H_{n-1}(x) +/// \f] +template +void testRecurrenceRelation(T x, T RelTolerance, T AbsTolerance) { + const unsigned MaxN = 127; + for (unsigned n = 1; n < MaxN; ++n) { + const T HermiteNext = std::experimental::hermite(n + 1, x); + const T HermiteNextRecurrence = + T(2) * x * std::experimental::hermite(n, x) - + T(2) * T(n) * std::experimental::hermite(n - 1, x); + const T Tolerance = AbsTolerance + std::abs(HermiteNext) * RelTolerance; + const T Error = std::abs(HermiteNextRecurrence - HermiteNext); + + if (std::isinf(HermiteNext)) + break; + assert(Error < Tolerance); + } +} + +template void testRecurrenceRelation(T RelTolerance, T AbsTolerance) { + 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) + testRecurrenceRelation(x, RelTolerance, AbsTolerance); +} + +/// \note Roots are taken from +/// Salzer, Herbert E., Ruth Zucker, and Ruth Capuano. +/// Table of the zeros and weight factors of the first twenty Hermite +/// polynomials. US Government Printing Office, 1952. +template std::vector getHermiteRoots(unsigned n) { + if (n == 0u) + return {}; + if (n == 1u) + return {T(0)}; + if (n == 2u) + return {T(0.707106781186548)}; + if (n == 3u) + return {T(0), + T(1.224744871391589)}; + if (n == 4u) + return {T(0.524647623275290), + T(1.650680123885785)}; + if (n == 5u) + return {T(0), T(0.958572464613819), + T(2.020182870456086)}; + if (n == 6u) + return {T(0.436077411927617), + T(1.335849074013697), + T(2.350604973674492)}; + if (n == 7u) + return {T(0), + T(0.816287882858965), + T(1.673551628767471), + T(2.651961356835233)}; + if (n == 8u) + return {T(0.381186990207322), + T(1.157193712446780), + T(1.981656756695843), + T(2.930637420257244)}; + if (n == 9u) + return {T(0), + T(0.723551018752838), + T(1.468553289216668), + T(2.266580584531843), + T(3.190993201781528)}; + if (n == 10u) + return {T(0.342901327223705), + T(1.036610829789514), + T(1.756683649299882), + T(2.532731674232790), + T(3.436159118837738)}; + if (n == 11u) + return {T(0), + T(0.65680956682100), + T(1.326557084494933), + T(2.025948015825755), + T(2.783290099781652), + T(3.668470846559583)}; + + if (n == 12u) + return {T(0.314240376254359), + T(0.947788391240164), + T(1.597682635152605), + T(2.279507080501060), + T(3.020637025120890), + T(3.889724897869782)}; + + if (n == 13u) + return {T(0), + T(0.605763879171060), + T(1.220055036590748), + T(1.853107651601512), + T(2.519735685678238), + T(3.246608978372410), + T(4.101337596178640)}; + + if (n == 14u) + return {T(0.29174551067256), + T(0.87871378732940), + T(1.47668273114114), + T(2.09518325850772), + T(2.74847072498540), + T(3.46265693360227), + T(4.30444857047363)}; + + if (n == 15u) + return {T(0.00000000000000), + T(0.56506958325558), + T(1.13611558521092), + T(1.71999257518649), + T(2.32573248617386), + T(2.96716692790560), + T(3.66995037340445), + T(4.49999070730939)}; + + if (n == 16u) + return {T(0.27348104613815), + T(0.82295144914466), + T(1.38025853919888), + T(1.95178799091625), + T(2.54620215784748), + T(3.17699916197996), + T(3.86944790486012), + T(4.68873893930582)}; + + if (n == 17u) + return {T(0), + T(0.5316330013427), + T(1.0676487257435), + T(1.6129243142212), + T(2.1735028266666), + T(2.7577629157039), + T(3.3789320911415), + T(4.0619466758755), + T(4.8713451936744)}; + if (n == 18u) + return {T(0.2582677505191), + T(0.7766829192674), + T(1.3009208583896), + T(1.8355316042616), + T(2.3862990891667), + T(2.9613775055316), + T(3.5737690684863), + T(4.2481178735681), + T(5.0483640088745)}; + if (n == 19u) + return {T(0), + T(0.5035201634239), + T(1.0103683871343), + T(1.5241706193935), + T(2.0492317098506), + T(2.5911337897945), + T(3.1578488183476), + T(3.7621873519640), + T(4.4285328066038), + T(5.2202716905375)}; + if (n == 20u) + return {T(0.2453407083009), + T(0.7374737285454), + T(1.2340762153953), + T(1.7385377121166), + T(2.2549740020893), + T(2.7888060584281), + T(3.347854567332), + T(3.9447640401156), + T(4.6036824495507), + T(5.3874808900112)}; + + return {}; +} + +/// \param [in] Tolerance of the root. This value must be smaller than +/// the smallest difference between adjacent roots in the given range +/// with n <= 20. +template void testHermiteRoots(T Tolerance) { + for (unsigned n = 0; n <= 20u; ++n) { + const auto Roots = getHermiteRoots(n); + for (T x : Roots) { + // the roots are symmetric: if x is a root, so is -x + if (x > T(0)) + assert(std::signbit(std::experimental::hermite(n, -x + Tolerance)) != + std::signbit(std::experimental::hermite(n, -x - Tolerance))); + assert(std::signbit(std::experimental::hermite(n, x + Tolerance)) != + std::signbit(std::experimental::hermite(n, x - Tolerance))); + } + } +} + +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-12l, 1e-12l); + + testRecurrenceRelation(1e-6f, 1e-6f); + testRecurrenceRelation(1e-9, 1e-9); + testRecurrenceRelation(1e-12l, 1e-12l); + + testHermiteRoots(1e-6f); + testHermiteRoots(1e-9); + testHermiteRoots(1e-10l); +#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,119 @@ +//===----------------------------------------------------------------------===// +// +// 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 +#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 ExpectedResult) { + if (std::isinf(ExpectedResult) && std::isinf(Result)) + return true; + + if (std::isnan(ExpectedResult) || std::isnan(Result)) + return false; + + const T Tolerance = + AbsTolerance + std::abs(ExpectedResult) * RelTolerance; + return std::abs(Result - ExpectedResult) < Tolerance; + }; + + const auto l0 = [](T) { return T(1); }; + const auto l1 = [](T 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,114 @@ +//===----------------------------------------------------------------------===// +// +// 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 +#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 ExpectedResult) { + if (std::isinf(ExpectedResult) && std::isinf(Result)) + return true; + + if (std::isnan(ExpectedResult) || std::isnan(Result)) + return false; + + const T Tolerance = + AbsTolerance + std::abs(ExpectedResult) * RelTolerance; + return std::abs(Result - ExpectedResult) < Tolerance; + }; + + const auto l0 = [](T) { return T(1); }; + const auto l1 = [](T 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; +}