[libc++][math] Mathematical Special Functions: Hermite Polynomial (#89982)

Summary:
Implementing the Hermite polynomials which are part of C++17's
mathematical special functions. The goal is to get early feedback which
will make implementing the other functions easier. Integration of
functions in chunks (e.g. `std::hermite` at first, then `std::laguerre`,
etc.) might make sense as well (also see note on boost.math below).

I started out from this abandoned merge request:
https://reviews.llvm.org/D58876 .

The C++23 standard defines them in-terms of `/* floating-point type */`
arguments. I have not looked into that.

Note, there is still an ongoing discussion on discourse whether
importing boost.math is an option.

Test Plan: 

Reviewers: 

Subscribers: 

Tasks: 

Tags: 


Differential Revision: https://phabricator.intern.facebook.com/D60251212

Author: PaulXiCao

libcxx/docs/ImplementationDefinedBehavior.rst

@@ -51,6 +51,17 @@ Libc++ determines that a stream is Unicode-capable terminal by:
   <http://eel.is/c++draft/print.fun#7>`_. This function is used for other
   ``std::print`` overloads that don't take an ``ostream&`` argument.
 
+`[sf.cmath] <https://wg21.link/sf.cmath>`_ Mathematical Special Functions: Large indices
+----------------------------------------------------------------------------------------
+
+Most functions within the Mathematical Special Functions section contain integral indices.
+The Standard specifies the result for larger indices as implementation-defined.
+Libc++ pursuits reasonable results by choosing the same formulas as for indices below that threshold.
+E.g.
+
+- ``std::hermite(unsigned n, T x)`` for ``n >= 128``
+
+
 Listed in the index of implementation-defined behavior
 ======================================================
 

libcxx/docs/Status/Cxx17.rst

@@ -41,6 +41,7 @@ Paper Status
 .. note::
 
    .. [#note-P0067] P0067: ``std::(to|from)_chars`` for integrals has been available since version 7.0. ``std::to_chars`` for ``float`` and ``double`` since version 14.0 ``std::to_chars`` for ``long double`` uses the implementation for ``double``.
+   .. [#note-P0226] P0226: Progress is tracked `here <https://https://libcxx.llvm.org/Status/SpecialMath.html>`_.
    .. [#note-P0607] P0607: The parts of P0607 that are not done are the ``<regex>`` bits.
    .. [#note-P0154] P0154: The required macros are only implemented as of clang 19.
    .. [#note-P0452] P0452: The changes to ``std::transform_inclusive_scan`` and ``std::transform_exclusive_scan`` have not yet been implemented.

libcxx/docs/Status/Cxx17Papers.csv

@@ -26,7 +26,7 @@
 "`P0013R1 <https://wg21.link/p0013r1>`__","LWG","Logical type traits rev 2","Kona","|Complete|","3.8"
 "","","","","",""
 "`P0024R2 <https://wg21.link/P0024R2>`__","LWG","The Parallelism TS Should be Standardized","Jacksonville","|Partial|",""
-"`P0226R1 <https://wg21.link/P0226R1>`__","LWG","Mathematical Special Functions for C++17","Jacksonville","",""
+"`P0226R1 <https://wg21.link/P0226R1>`__","LWG","Mathematical Special Functions for C++17","Jacksonville","|In Progress| [#note-P0226]_",""
 "`P0220R1 <https://wg21.link/P0220R1>`__","LWG","Adopt Library Fundamentals V1 TS Components for C++17","Jacksonville","|Complete|","16.0"
 "`P0218R1 <https://wg21.link/P0218R1>`__","LWG","Adopt the File System TS for C++17","Jacksonville","|Complete|","7.0"
 "`P0033R1 <https://wg21.link/P0033R1>`__","LWG","Re-enabling shared_from_this","Jacksonville","|Complete|","3.9"

libcxx/docs/Status/SpecialMath.rst

@@ -0,0 +1,35 @@
+.. special-math-status:
+
+======================================================
+libc++ Mathematical Special Functions Status (P0226R1)
+======================================================
+
+.. include:: ../Helpers/Styles.rst
+
+.. contents::
+  :local:
+
+Overview
+========
+
+This document contains the status of the C++17 mathematical special functions implementation in libc++.
+It is used to track both the status of the sub-projects of the effort and who is assigned to these sub-projects.
+This avoids duplicating effort.
+
+If you are interested in contributing to this effort, please send a message
+to the #libcxx channel in the LLVM discord. Please *do not* start working
+on any items below that has already been assigned to someone else.
+
+Sub-projects in the Implementation Effort
+=========================================
+
+.. csv-table::
+  :file: SpecialMathProjects.csv
+  :header-rows: 1
+  :widths: auto
+
+Paper and Issue Status
+======================
+
+The underlying paper is `Mathematical Special Functions for C++17 (P0226) <https://wg21.link/P0226>`_ and is included in C++17.
+Implementation is *In Progress*.

libcxx/docs/Status/SpecialMathProjects.csv

@@ -0,0 +1,22 @@
+Section,Description,Assignee,Complete
+| `[sf.cmath.assoc.laguerre] <https://wg21.link/sf.cmath.assoc.laguerre>`_, std::assoc_laguerre, None, |Not Started|
+| `[sf.cmath.assoc.legendre] <https://wg21.link/sf.cmath.assoc.legendre>`_, std::assoc_legendre, None, |Not Started|
+| `[sf.cmath.beta] <https://wg21.link/sf.cmath.beta>`_, std::beta, None, |Not Started|
+| `[sf.cmath.comp.ellint.1] <https://wg21.link/sf.cmath.comp.ellint.1>`_, std::comp_ellint_1, None, |Not Started|
+| `[sf.cmath.comp.ellint.2] <https://wg21.link/sf.cmath.comp.ellint.2>`_, std::comp_ellint_2, None, |Not Started|
+| `[sf.cmath.comp.ellint.3] <https://wg21.link/sf.cmath.comp.ellint.3>`_, std::comp_ellint_3, None, |Not Started|
+| `[sf.cmath.cyl.bessel.i] <https://wg21.link/sf.cmath.cyl.bessel.i>`_, std::cyl_bessel_i, None, |Not Started|
+| `[sf.cmath.cyl.bessel.j] <https://wg21.link/sf.cmath.cyl.bessel.j>`_, std::cyl_bessel_j, None, |Not Started|
+| `[sf.cmath.cyl.bessel.k] <https://wg21.link/sf.cmath.cyl.bessel.k>`_, std::cyl_bessel_k, None, |Not Started|
+| `[sf.cmath.cyl.neumann] <https://wg21.link/sf.cmath.cyl.neumann>`_, std::cyl_neumann, None, |Not Started|
+| `[sf.cmath.ellint.1] <https://wg21.link/sf.cmath.ellint.1>`_, std::ellint_1, None, |Not Started|
+| `[sf.cmath.ellint.2] <https://wg21.link/sf.cmath.ellint.2>`_, std::ellint_2, None, |Not Started|
+| `[sf.cmath.ellint.3] <https://wg21.link/sf.cmath.ellint.3>`_, std::ellint_3, None, |Not Started|
+| `[sf.cmath.expint] <https://wg21.link/sf.cmath.expint>`_, std::expint, None, |Not Started|
+| `[sf.cmath.hermite] <https://wg21.link/sf.cmath.hermite>`_, std::hermite, Paul Xi Cao, |Complete|
+| `[sf.cmath.laguerre] <https://wg21.link/sf.cmath.laguerre>`_, std::laguerre, None, |Not Started|
+| `[sf.cmath.legendre] <https://wg21.link/sf.cmath.legendre>`_, std::legendre, None, |Not Started|
+| `[sf.cmath.riemann.zeta] <https://wg21.link/sf.cmath.riemann.zeta>`_, std::riemann_zeta, None, |Not Started|
+| `[sf.cmath.sph.bessel] <https://wg21.link/sf.cmath.sph.bessel>`_, std::sph_bessel, None, |Not Started|
+| `[sf.cmath.sph.legendre] <https://wg21.link/sf.cmath.sph.legendre>`_, std::sph_legendre, None, |Not Started|
+| `[sf.cmath.sph.neumann] <https://wg21.link/sf.cmath.sph.neumann>`_, std::sph_neumann, None, |Not Started|

libcxx/docs/index.rst

@@ -53,6 +53,7 @@ Getting Started with libc++
    Status/PSTL
    Status/Ranges
    Status/Spaceship
+   Status/SpecialMath
    Status/Zip
 
 

libcxx/include/CMakeLists.txt

@@ -509,6 +509,7 @@ set(files
   __math/remainder.h
   __math/roots.h
   __math/rounding_functions.h
+  __math/special_functions.h
   __math/traits.h
   __math/trigonometric_functions.h
   __mbstate_t.h

libcxx/include/__math/special_functions.h

@@ -0,0 +1,84 @@
+// -*- 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
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef _LIBCPP___MATH_SPECIAL_FUNCTIONS_H
+#define _LIBCPP___MATH_SPECIAL_FUNCTIONS_H
+
+#include <__config>
+#include <__math/copysign.h>
+#include <__math/traits.h>
+#include <__type_traits/enable_if.h>
+#include <__type_traits/is_integral.h>
+#include <limits>
+
+#if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
+#  pragma GCC system_header
+#endif
+
+_LIBCPP_BEGIN_NAMESPACE_STD
+
+#if _LIBCPP_STD_VER >= 17
+
+template <class _Real>
+_LIBCPP_HIDE_FROM_ABI _Real __hermite(unsigned __n, _Real __x) {
+  // The Hermite polynomial H_n(x).
+  // The implementation is based on the recurrence formula: H_{n+1}(x) = 2x H_n(x) - 2n H_{n-1}.
+  // Press, William H., et al. Numerical recipes 3rd edition: The art of scientific computing.
+  // Cambridge university press, 2007, p. 183.
+
+  // NOLINTBEGIN(readability-identifier-naming)
+  if (__math::isnan(__x))
+    return __x;
+
+  _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;
+  }
+
+  if (!__math::isfinite(__H_n)) {
+    // Overflow occured. Two possible cases:
+    //    n is odd:  return infinity of the same sign as x.
+    //    n is even: return +Inf
+    _Real __inf = std::numeric_limits<_Real>::infinity();
+    return (__n & 1) ? __math::copysign(__inf, __x) : __inf;
+  }
+  return __H_n;
+  // NOLINTEND(readability-identifier-naming)
+}
+
+inline _LIBCPP_HIDE_FROM_ABI double hermite(unsigned __n, double __x) { return std::__hermite(__n, __x); }
+
+inline _LIBCPP_HIDE_FROM_ABI float hermite(unsigned __n, float __x) {
+  // use double internally -- float is too prone to overflow!
+  return static_cast<float>(std::hermite(__n, static_cast<double>(__x)));
+}
+
+inline _LIBCPP_HIDE_FROM_ABI long double hermite(unsigned __n, long double __x) { return std::__hermite(__n, __x); }
+
+inline _LIBCPP_HIDE_FROM_ABI float hermitef(unsigned __n, float __x) { return std::hermite(__n, __x); }
+
+inline _LIBCPP_HIDE_FROM_ABI long double hermitel(unsigned __n, long double __x) { return std::hermite(__n, __x); }
+
+template <class _Integer, std::enable_if_t<std::is_integral_v<_Integer>, int> = 0>
+_LIBCPP_HIDE_FROM_ABI double hermite(unsigned __n, _Integer __x) {
+  return std::hermite(__n, static_cast<double>(__x));
+}
+
+#endif // _LIBCPP_STD_VER >= 17
+
+_LIBCPP_END_NAMESPACE_STD
+
+#endif // _LIBCPP___MATH_SPECIAL_FUNCTIONS_H

libcxx/include/cmath

@@ -204,6 +204,14 @@ floating_point fmin (arithmetic x, arithmetic y);
 float          fminf(float x, float y);
 long double    fminl(long double x, long double y);
 
+double         hermite(unsigned n, double x);                    // C++17
+float          hermite(unsigned n, float x);                     // C++17
+long double    hermite(unsigned n, long double x);               // C++17
+float          hermitef(unsigned n, float x);                    // C++17
+long double    hermitel(unsigned n, long double x);              // C++17
+template <class Integer>
+double         hermite(unsigned n, Integer x);                   // C++17
+
 floating_point hypot (arithmetic x, arithmetic y);
 float          hypotf(float x, float y);
 long double    hypotl(long double x, long double y);
@@ -315,6 +323,7 @@ constexpr long double lerp(long double a, long double b, long double t) noexcept
 #include <limits>
 #include <version>
 
+#include <__math/special_functions.h>
 #include <math.h>
 
 #ifndef _LIBCPP_MATH_H

libcxx/include/module.modulemap

@@ -1485,6 +1485,7 @@ module std_private_math_modulo                          [system] { header "__mat
 module std_private_math_remainder                       [system] { header "__math/remainder.h" }
 module std_private_math_roots                           [system] { header "__math/roots.h" }
 module std_private_math_rounding_functions              [system] { header "__math/rounding_functions.h" }
+module std_private_math_special_functions               [system] { header "__math/special_functions.h" }
 module std_private_math_traits                          [system] { header "__math/traits.h" }
 module std_private_math_trigonometric_functions         [system] { header "__math/trigonometric_functions.h" }
 

libcxx/modules/std/cmath.inc

@@ -334,12 +334,14 @@ export namespace std {
   using std::expint;
   using std::expintf;
   using std::expintl;
+#endif
 
   // [sf.cmath.hermite], Hermite polynomials
   using std::hermite;
   using std::hermitef;
   using std::hermitel;
 
+#if 0
   // [sf.cmath.laguerre], Laguerre polynomials
   using std::laguerre;
   using std::laguerref;