|
| 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, 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, 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 |
0 commit comments