| 1 | // -*- C++ -*- |
| 2 | //===----------------------------------------------------------------------===// |
| 3 | // |
| 4 | // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. |
| 5 | // See https://llvm.org/LICENSE.txt for license information. |
| 6 | // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception |
| 7 | // |
| 8 | //===----------------------------------------------------------------------===// |
| 9 | |
| 10 | #ifndef _LIBCPP___MATH_SPECIAL_FUNCTIONS_H |
| 11 | #define _LIBCPP___MATH_SPECIAL_FUNCTIONS_H |
| 12 | |
| 13 | #include <__config> |
| 14 | #include <__math/copysign.h> |
| 15 | #include <__math/traits.h> |
| 16 | #include <__type_traits/enable_if.h> |
| 17 | #include <__type_traits/is_integral.h> |
| 18 | #include <limits> |
| 19 | |
| 20 | #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER) |
| 21 | # pragma GCC system_header |
| 22 | #endif |
| 23 | |
| 24 | _LIBCPP_BEGIN_NAMESPACE_STD |
| 25 | |
| 26 | #if _LIBCPP_STD_VER >= 17 |
| 27 | |
| 28 | template <class _Real> |
| 29 | _LIBCPP_HIDE_FROM_ABI _Real __hermite(unsigned __n, _Real __x) { |
| 30 | // The Hermite polynomial H_n(x). |
| 31 | // The implementation is based on the recurrence formula: H_{n+1}(x) = 2x H_n(x) - 2n H_{n-1}. |
| 32 | // Press, William H., et al. Numerical recipes 3rd edition: The art of scientific computing. |
| 33 | // Cambridge university press, 2007, p. 183. |
| 34 | |
| 35 | // NOLINTBEGIN(readability-identifier-naming) |
| 36 | if (__math::isnan(__x)) |
| 37 | return __x; |
| 38 | |
| 39 | _Real __H_0{1}; |
| 40 | if (__n == 0) |
| 41 | return __H_0; |
| 42 | |
| 43 | _Real __H_n_prev = __H_0; |
| 44 | _Real __H_n = 2 * __x; |
| 45 | for (unsigned __i = 1; __i < __n; ++__i) { |
| 46 | _Real __H_n_next = 2 * (__x * __H_n - __i * __H_n_prev); |
| 47 | __H_n_prev = __H_n; |
| 48 | __H_n = __H_n_next; |
| 49 | } |
| 50 | |
| 51 | if (!__math::isfinite(__H_n)) { |
| 52 | // Overflow occured. Two possible cases: |
| 53 | // n is odd: return infinity of the same sign as x. |
| 54 | // n is even: return +Inf |
| 55 | _Real __inf = std::numeric_limits<_Real>::infinity(); |
| 56 | return (__n & 1) ? __math::copysign(__inf, __x) : __inf; |
| 57 | } |
| 58 | return __H_n; |
| 59 | // NOLINTEND(readability-identifier-naming) |
| 60 | } |
| 61 | |
| 62 | inline _LIBCPP_HIDE_FROM_ABI double hermite(unsigned __n, double __x) { return std::__hermite(__n, __x); } |
| 63 | |
| 64 | inline _LIBCPP_HIDE_FROM_ABI float hermite(unsigned __n, float __x) { |
| 65 | // use double internally -- float is too prone to overflow! |
| 66 | return static_cast<float>(std::hermite(__n, x: static_cast<double>(__x))); |
| 67 | } |
| 68 | |
| 69 | inline _LIBCPP_HIDE_FROM_ABI long double hermite(unsigned __n, long double __x) { return std::__hermite(__n, __x); } |
| 70 | |
| 71 | inline _LIBCPP_HIDE_FROM_ABI float hermitef(unsigned __n, float __x) { return std::hermite(__n, __x); } |
| 72 | |
| 73 | inline _LIBCPP_HIDE_FROM_ABI long double hermitel(unsigned __n, long double __x) { return std::hermite(__n, __x); } |
| 74 | |
| 75 | template <class _Integer, std::enable_if_t<std::is_integral_v<_Integer>, int> = 0> |
| 76 | _LIBCPP_HIDE_FROM_ABI double hermite(unsigned __n, _Integer __x) { |
| 77 | return std::hermite(__n, x: static_cast<double>(__x)); |
| 78 | } |
| 79 | |
| 80 | #endif // _LIBCPP_STD_VER >= 17 |
| 81 | |
| 82 | _LIBCPP_END_NAMESPACE_STD |
| 83 | |
| 84 | #endif // _LIBCPP___MATH_SPECIAL_FUNCTIONS_H |
| 85 | |