[libcxx-commits] [libcxx] [libc++][math] Mathematical Special Functions: Hermite Polynomial (PR #89982)

[via libcxx-commits]

================
@@ -765,6 +766,53 @@ _LIBCPP_HIDE_FROM_ABI _LIBCPP_CONSTEXPR_SINCE_CXX20 _Tp __constexpr_scalbn(_Tp _
   return __builtin_scalbn(__x, __exp);
 }
 
+#if _LIBCPP_STD_VER >= 17
+
+/// \return the Hermite polynomial \f$ H_n(x) \f$
+/// \note The implementation is based on the recurrence formula
+/// \f[
+///   H_{n+1}(x) = 2x H_n(x) - 2n H_{n-1}
+/// \f]
+/// Press, William H., et al. Numerical recipes 3rd edition: The art of
+/// scientific computing. Cambridge university press, 2007, p. 183.
+template <class _Real>
+_LIBCPP_HIDE_FROM_ABI _Real __hermite(unsigned __n, _Real __x) {
+  if (std::isnan(__x))
+    return std::numeric_limits<_Real>::quiet_NaN();
+
+  _Real __H_0{1};
+  if (__n == 0)
+    return __H_0;
+
+  _Real __H_n_prev = __H_0;
+  _Real __H_n      = 2 * __x;
+  for (unsigned __i = 1; __i < __n; ++__i) {
+    _Real __H_n_next = 2 * (__x * __H_n - __i * __H_n_prev);
+    __H_n_prev       = __H_n;
+    __H_n            = __H_n_next;
+  }
+  return __H_n;
----------------
PaulXiCao wrote:

I looked some more into possible overflows. First thing to notice is that double overflows at ~1e308.
Hermite polynomial H_127 does not reach such large values inbetween any roots but only for larger x:
* largest root is at x=~16
* function values previous to the last root are significantly below the overflow threshold, i.e. |H_127(x)| < ~1e190 for |x| < ~16
* overflow threshold is reached at roughly: H_127(x=140)~1e308

For polynomials of smaller order it is even more extreme (e.g. largest root is smaller and overflow happens at even larger x).
I would not take polynomials of higher order into this discussion as their result is deemed implementation-defined by the Standard. Thus, possibly slowing down code for n<128 because of arguments for n>=128 seems strange.

My conclusion: I believe your initial suggestion is enough (e.g. checking the value of `__H_n` *after* the loop). Using `!std::isfinite()` also seems the approriate check here.

I will probably push a change similar to `hermite2` as in this [godbolt](https://godbolt.org/z/6xojME18W). (`hermite` is the current implementation, and `hermite1` the one you proposed in the previous godbolt.)

https://github.com/llvm/llvm-project/pull/89982