[libcxx-commits] [libcxx] af0d731 - [libc++][math] Mathematical Special Functions: Hermite Polynomial (#89982)
via libcxx-commits
libcxx-commits at lists.llvm.org
Sat Jul 20 08:50:09 PDT 2024
Author: PaulXiCao
Date: 2024-07-20T17:50:05+02:00
New Revision: af0d731b12983a649e07f25037ee0e16a68ca470
URL: https://github.com/llvm/llvm-project/commit/af0d731b12983a649e07f25037ee0e16a68ca470
DIFF: https://github.com/llvm/llvm-project/commit/af0d731b12983a649e07f25037ee0e16a68ca470.diff
LOG: [libc++][math] Mathematical Special Functions: Hermite Polynomial (#89982)
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.
Added:
libcxx/docs/Status/SpecialMath.rst
libcxx/docs/Status/SpecialMathProjects.csv
libcxx/include/__math/special_functions.h
libcxx/test/std/numerics/c.math/hermite.pass.cpp
Modified:
libcxx/docs/ImplementationDefinedBehavior.rst
libcxx/docs/Status/Cxx17.rst
libcxx/docs/Status/Cxx17Papers.csv
libcxx/docs/index.rst
libcxx/include/CMakeLists.txt
libcxx/include/cmath
libcxx/include/module.modulemap
libcxx/modules/std/cmath.inc
libcxx/utils/libcxx/test/modules.py
Removed:
################################################################################
diff --git a/libcxx/docs/ImplementationDefinedBehavior.rst b/libcxx/docs/ImplementationDefinedBehavior.rst
index 3000bb7cfa468..f0ef733fc2c55 100644
--- a/libcxx/docs/ImplementationDefinedBehavior.rst
+++ b/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
======================================================
diff --git a/libcxx/docs/Status/Cxx17.rst b/libcxx/docs/Status/Cxx17.rst
index d4426afa81638..ad4f8576f03db 100644
--- a/libcxx/docs/Status/Cxx17.rst
+++ b/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.
diff --git a/libcxx/docs/Status/Cxx17Papers.csv b/libcxx/docs/Status/Cxx17Papers.csv
index 2e560cfe0d576..6c657d51f5c7e 100644
--- a/libcxx/docs/Status/Cxx17Papers.csv
+++ b/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"
diff --git a/libcxx/docs/Status/SpecialMath.rst b/libcxx/docs/Status/SpecialMath.rst
new file mode 100644
index 0000000000000..fcc9f03e3ae64
--- /dev/null
+++ b/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*.
diff --git a/libcxx/docs/Status/SpecialMathProjects.csv b/libcxx/docs/Status/SpecialMathProjects.csv
new file mode 100644
index 0000000000000..f964e79de91d3
--- /dev/null
+++ b/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|
diff --git a/libcxx/docs/index.rst b/libcxx/docs/index.rst
index 69a9e575cfe7c..4bca3ccc8fa06 100644
--- a/libcxx/docs/index.rst
+++ b/libcxx/docs/index.rst
@@ -53,6 +53,7 @@ Getting Started with libc++
Status/PSTL
Status/Ranges
Status/Spaceship
+ Status/SpecialMath
Status/Zip
diff --git a/libcxx/include/CMakeLists.txt b/libcxx/include/CMakeLists.txt
index 1a4d9c7070f14..32579272858a8 100644
--- a/libcxx/include/CMakeLists.txt
+++ b/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
diff --git a/libcxx/include/__math/special_functions.h b/libcxx/include/__math/special_functions.h
new file mode 100644
index 0000000000000..0b1c753a659ad
--- /dev/null
+++ b/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
diff --git a/libcxx/include/cmath b/libcxx/include/cmath
index 7a87e35c84603..3c22604a683c3 100644
--- a/libcxx/include/cmath
+++ b/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
diff --git a/libcxx/include/module.modulemap b/libcxx/include/module.modulemap
index 3443fbc3347a3..13d0dce34d97e 100644
--- a/libcxx/include/module.modulemap
+++ b/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" }
diff --git a/libcxx/modules/std/cmath.inc b/libcxx/modules/std/cmath.inc
index a463c1e3ccf86..fe8ac773c9d1c 100644
--- a/libcxx/modules/std/cmath.inc
+++ b/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;
diff --git a/libcxx/test/std/numerics/c.math/hermite.pass.cpp b/libcxx/test/std/numerics/c.math/hermite.pass.cpp
new file mode 100644
index 0000000000000..08fbd5c3283c1
--- /dev/null
+++ b/libcxx/test/std/numerics/c.math/hermite.pass.cpp
@@ -0,0 +1,341 @@
+//===----------------------------------------------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+// UNSUPPORTED: c++03, c++11, c++14
+
+// <cmath>
+
+// double hermite(unsigned n, double x);
+// float hermite(unsigned n, float x);
+// long double hermite(unsigned n, long double x);
+// float hermitef(unsigned n, float x);
+// long double hermitel(unsigned n, long double x);
+// template <class Integer>
+// double hermite(unsigned n, Integer x);
+
+#include <array>
+#include <cassert>
+#include <cmath>
+#include <limits>
+#include <vector>
+
+#include "type_algorithms.h"
+
+inline constexpr unsigned g_max_n = 128;
+
+template <class T>
+std::array<T, 11> sample_points() {
+ return {-12.34, -7.42, -1.0, -0.5, -0.1, 0.0, 0.1, 0.5, 1.0, 5.67, 15.67};
+}
+
+template <class Real>
+class CompareFloatingValues {
+private:
+ Real abs_tol;
+ Real rel_tol;
+
+public:
+ CompareFloatingValues() {
+ abs_tol = []() -> Real {
+ if (std::is_same_v<Real, float>)
+ return 1e-5f;
+ else if (std::is_same_v<Real, double>)
+ return 1e-11;
+ else
+ return 1e-12l;
+ }();
+
+ rel_tol = abs_tol;
+ }
+
+ bool operator()(Real result, Real expected) const {
+ if (std::isinf(expected) && std::isinf(result))
+ return result == expected;
+
+ if (std::isnan(expected) || std::isnan(result))
+ return false;
+
+ Real tol = abs_tol + std::abs(expected) * rel_tol;
+ return std::abs(result - expected) < tol;
+ }
+};
+
+// 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 <class T>
+std::vector<T> get_roots(unsigned n) {
+ switch (n) {
+ case 0:
+ return {};
+ case 1:
+ return {T(0)};
+ case 2:
+ return {T(0.707106781186548)};
+ case 3:
+ return {T(0), T(1.224744871391589)};
+ case 4:
+ return {T(0.524647623275290), T(1.650680123885785)};
+ case 5:
+ return {T(0), T(0.958572464613819), T(2.020182870456086)};
+ case 6:
+ return {T(0.436077411927617), T(1.335849074013697), T(2.350604973674492)};
+ case 7:
+ return {T(0), T(0.816287882858965), T(1.673551628767471), T(2.651961356835233)};
+ case 8:
+ return {T(0.381186990207322), T(1.157193712446780), T(1.981656756695843), T(2.930637420257244)};
+ case 9:
+ return {T(0), T(0.723551018752838), T(1.468553289216668), T(2.266580584531843), T(3.190993201781528)};
+ case 10:
+ return {
+ T(0.342901327223705), T(1.036610829789514), T(1.756683649299882), T(2.532731674232790), T(3.436159118837738)};
+ case 11:
+ return {T(0),
+ T(0.65680956682100),
+ T(1.326557084494933),
+ T(2.025948015825755),
+ T(2.783290099781652),
+ T(3.668470846559583)};
+
+ case 12:
+ return {T(0.314240376254359),
+ T(0.947788391240164),
+ T(1.597682635152605),
+ T(2.279507080501060),
+ T(3.020637025120890),
+ T(3.889724897869782)};
+
+ case 13:
+ return {T(0),
+ T(0.605763879171060),
+ T(1.220055036590748),
+ T(1.853107651601512),
+ T(2.519735685678238),
+ T(3.246608978372410),
+ T(4.101337596178640)};
+
+ case 14:
+ return {T(0.29174551067256),
+ T(0.87871378732940),
+ T(1.47668273114114),
+ T(2.09518325850772),
+ T(2.74847072498540),
+ T(3.46265693360227),
+ T(4.30444857047363)};
+
+ case 15:
+ 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)};
+
+ case 16:
+ 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)};
+
+ case 17:
+ 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)};
+
+ case 18:
+ 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)};
+
+ case 19:
+ 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)};
+
+ case 20:
+ 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)};
+
+ default: // polynom degree n>20 is unsupported
+ assert(false);
+ return {T(-42)};
+ }
+}
+
+template <class Real>
+void test() {
+ { // checks if NaNs are reported correctly (i.e. output == input for input == NaN)
+ using nl = std::numeric_limits<Real>;
+ for (Real NaN : {nl::quiet_NaN(), nl::signaling_NaN()})
+ for (unsigned n = 0; n < g_max_n; ++n)
+ assert(std::isnan(std::hermite(n, NaN)));
+ }
+
+ { // simple sample points for n=0..127 should not produce NaNs.
+ for (Real x : sample_points<Real>())
+ for (unsigned n = 0; n < g_max_n; ++n)
+ assert(!std::isnan(std::hermite(n, x)));
+ }
+
+ { // checks std::hermite(n, x) for n=0..5 against analytic polynoms
+ const auto h0 = [](Real) -> Real { return 1; };
+ const auto h1 = [](Real y) -> Real { return 2 * y; };
+ const auto h2 = [](Real y) -> Real { return 4 * y * y - 2; };
+ const auto h3 = [](Real y) -> Real { return y * (8 * y * y - 12); };
+ const auto h4 = [](Real y) -> Real { return (16 * std::pow(y, 4) - 48 * y * y + 12); };
+ const auto h5 = [](Real y) -> Real { return y * (32 * std::pow(y, 4) - 160 * y * y + 120); };
+
+ for (Real x : sample_points<Real>()) {
+ const CompareFloatingValues<Real> compare;
+ assert(compare(std::hermite(0, x), h0(x)));
+ assert(compare(std::hermite(1, x), h1(x)));
+ assert(compare(std::hermite(2, x), h2(x)));
+ assert(compare(std::hermite(3, x), h3(x)));
+ assert(compare(std::hermite(4, x), h4(x)));
+ assert(compare(std::hermite(5, x), h5(x)));
+ }
+ }
+
+ { // checks std::hermitef for bitwise equality with std::hermite(unsigned, float)
+ if constexpr (std::is_same_v<Real, float>)
+ for (unsigned n = 0; n < g_max_n; ++n)
+ for (float x : sample_points<float>())
+ assert(std::hermite(n, x) == std::hermitef(n, x));
+ }
+
+ { // checks std::hermitel for bitwise equality with std::hermite(unsigned, long double)
+ if constexpr (std::is_same_v<Real, long double>)
+ for (unsigned n = 0; n < g_max_n; ++n)
+ for (long double x : sample_points<long double>())
+ assert(std::hermite(n, x) == std::hermitel(n, x));
+ }
+
+ { // Checks if the characteristic recurrence relation holds: H_{n+1}(x) = 2x H_n(x) - 2n H_{n-1}(x)
+ for (Real x : sample_points<Real>()) {
+ for (unsigned n = 1; n < g_max_n - 1; ++n) {
+ Real H_next = std::hermite(n + 1, x);
+ Real H_next_recurrence = 2 * (x * std::hermite(n, x) - n * std::hermite(n - 1, x));
+
+ if (std::isinf(H_next))
+ break;
+ const CompareFloatingValues<Real> compare;
+ assert(compare(H_next, H_next_recurrence));
+ }
+ }
+ }
+
+ { // sanity checks: hermite polynoms need to change signs at (simple) roots. checked upto order n<=20.
+
+ // root tolerance: must be smaller than the smallest
diff erence between adjacent roots
+ Real tol = []() -> Real {
+ if (std::is_same_v<Real, float>)
+ return 1e-5f;
+ else if (std::is_same_v<Real, double>)
+ return 1e-9;
+ else
+ return 1e-10l;
+ }();
+
+ const auto is_sign_change = [tol](unsigned n, Real x) -> bool {
+ return std::hermite(n, x - tol) * std::hermite(n, x + tol) < 0;
+ };
+
+ for (unsigned n = 0; n <= 20u; ++n) {
+ for (Real x : get_roots<Real>(n)) {
+ // the roots are symmetric: if x is a root, so is -x
+ if (x > 0)
+ assert(is_sign_change(n, -x));
+ assert(is_sign_change(n, x));
+ }
+ }
+ }
+
+ { // check input infinity is handled correctly
+ Real inf = std::numeric_limits<Real>::infinity();
+ for (unsigned n = 1; n < g_max_n; ++n) {
+ assert(std::hermite(n, +inf) == inf);
+ assert(std::hermite(n, -inf) == ((n & 1) ? -inf : inf));
+ }
+ }
+
+ { // check: if overflow occurs that it is mapped to the correct infinity
+ if constexpr (std::is_same_v<Real, double>) {
+ // Q: Why only double?
+ // A: The numeric values (e.g. overflow threshold `n`) below are
diff erent for other types.
+ static_assert(sizeof(double) == 8);
+ for (unsigned n = 0; n < g_max_n; ++n) {
+ // Q: Why n=111 and x=300?
+ // A: Both are chosen s.t. the first overlow occurs for some `n<g_max_n`.
+ if (n < 111) {
+ assert(std::isfinite(std::hermite(n, +300.0)));
+ assert(std::isfinite(std::hermite(n, -300.0)));
+ } else {
+ double inf = std::numeric_limits<double>::infinity();
+ assert(std::hermite(n, +300.0) == inf);
+ assert(std::hermite(n, -300.0) == ((n & 1) ? -inf : inf));
+ }
+ }
+ }
+ }
+}
+
+struct TestFloat {
+ template <class Real>
+ void operator()() {
+ test<Real>();
+ }
+};
+
+struct TestInt {
+ template <class Integer>
+ void operator()() {
+ // checks that std::hermite(unsigned, Integer) actually wraps std::hermite(unsigned, double)
+ for (unsigned n = 0; n < g_max_n; ++n)
+ for (Integer x : {-42, -7, -5, -1, 0, 1, 5, 7, 42})
+ assert(std::hermite(n, x) == std::hermite(n, static_cast<double>(x)));
+ }
+};
+
+int main() {
+ types::for_each(types::floating_point_types(), TestFloat());
+ types::for_each(types::type_list<short, int, long, long long>(), TestInt());
+}
diff --git a/libcxx/utils/libcxx/test/modules.py b/libcxx/utils/libcxx/test/modules.py
index aab7651c7bb03..b7758dc9a41ee 100644
--- a/libcxx/utils/libcxx/test/modules.py
+++ b/libcxx/utils/libcxx/test/modules.py
@@ -76,6 +76,13 @@
# This declaration is in the ostream header.
ExtraDeclarations["system_error"] = ["std::operator<<"]
+# TODO MODULES avoid this work-around
+# This is a work-around for the special math functions. They are declared in
+# __math/special_functions.h. Adding this as an ExtraHeader works for the std
+# module. However these functions are special; they are not available in the
+# global namespace.
+ExtraDeclarations["cmath"] = ["std::hermite", "std::hermitef", "std::hermitel"]
+
### ExtraHeader
# Adds extra headers file to scan
More information about the libcxx-commits
mailing list