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: `_. This function is used for other ``std::print`` overloads that don't take an ``ostream&`` argument. +`[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 `_. .. [#note-P0607] P0607: The parts of P0607 that are not done are the ```` 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 `__","LWG","Logical type traits rev 2","Kona","|Complete|","3.8" "","","","","","" "`P0024R2 `__","LWG","The Parallelism TS Should be Standardized","Jacksonville","|Partial|","" -"`P0226R1 `__","LWG","Mathematical Special Functions for C++17","Jacksonville","","" +"`P0226R1 `__","LWG","Mathematical Special Functions for C++17","Jacksonville","|In Progress| [#note-P0226]_","" "`P0220R1 `__","LWG","Adopt Library Fundamentals V1 TS Components for C++17","Jacksonville","|Complete|","16.0" "`P0218R1 `__","LWG","Adopt the File System TS for C++17","Jacksonville","|Complete|","7.0" "`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) `_ 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] `_, std::assoc_laguerre, None, |Not Started| +| `[sf.cmath.assoc.legendre] `_, std::assoc_legendre, None, |Not Started| +| `[sf.cmath.beta] `_, std::beta, None, |Not Started| +| `[sf.cmath.comp.ellint.1] `_, std::comp_ellint_1, None, |Not Started| +| `[sf.cmath.comp.ellint.2] `_, std::comp_ellint_2, None, |Not Started| +| `[sf.cmath.comp.ellint.3] `_, std::comp_ellint_3, None, |Not Started| +| `[sf.cmath.cyl.bessel.i] `_, std::cyl_bessel_i, None, |Not Started| +| `[sf.cmath.cyl.bessel.j] `_, std::cyl_bessel_j, None, |Not Started| +| `[sf.cmath.cyl.bessel.k] `_, std::cyl_bessel_k, None, |Not Started| +| `[sf.cmath.cyl.neumann] `_, std::cyl_neumann, None, |Not Started| +| `[sf.cmath.ellint.1] `_, std::ellint_1, None, |Not Started| +| `[sf.cmath.ellint.2] `_, std::ellint_2, None, |Not Started| +| `[sf.cmath.ellint.3] `_, std::ellint_3, None, |Not Started| +| `[sf.cmath.expint] `_, std::expint, None, |Not Started| +| `[sf.cmath.hermite] `_, std::hermite, Paul Xi Cao, |Complete| +| `[sf.cmath.laguerre] `_, std::laguerre, None, |Not Started| +| `[sf.cmath.legendre] `_, std::legendre, None, |Not Started| +| `[sf.cmath.riemann.zeta] `_, std::riemann_zeta, None, |Not Started| +| `[sf.cmath.sph.bessel] `_, std::sph_bessel, None, |Not Started| +| `[sf.cmath.sph.legendre] `_, std::sph_legendre, None, |Not Started| +| `[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 07dd25604a9c7..55956b9ed2644 100644 --- a/libcxx/include/CMakeLists.txt +++ b/libcxx/include/CMakeLists.txt @@ -508,6 +508,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 + +#if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER) +# pragma GCC system_header +#endif + +_LIBCPP_BEGIN_NAMESPACE_STD + +#if _LIBCPP_STD_VER >= 17 + +template +_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(std::hermite(__n, static_cast(__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 , int> = 0> +_LIBCPP_HIDE_FROM_ABI double hermite(unsigned __n, _Integer __x) { + return std::hermite(__n, static_cast(__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 +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 #include +#include <__math/special_functions.h> #include #ifndef _LIBCPP_MATH_H diff --git a/libcxx/include/module.modulemap b/libcxx/include/module.modulemap index 4ad506781c489..249c138012a3e 100644 --- a/libcxx/include/module.modulemap +++ b/libcxx/include/module.modulemap @@ -1480,6 +1480,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 + +// + +// 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 +// double hermite(unsigned n, Integer x); + +#include +#include +#include +#include +#include + +#include "type_algorithms.h" + +inline constexpr unsigned g_max_n = 128; + +template +std::array 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 CompareFloatingValues { +private: + Real abs_tol; + Real rel_tol; + +public: + CompareFloatingValues() { + abs_tol = []() -> Real { + if (std::is_same_v) + return 1e-5f; + else if (std::is_same_v) + 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 +std::vector 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 +void test() { + { // checks if NaNs are reported correctly (i.e. output == input for input == NaN) + using nl = std::numeric_limits; + 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()) + 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()) { + const CompareFloatingValues 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) + for (unsigned n = 0; n < g_max_n; ++n) + for (float x : sample_points()) + 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) + for (unsigned n = 0; n < g_max_n; ++n) + for (long double x : sample_points()) + 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()) { + 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 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 difference between adjacent roots + Real tol = []() -> Real { + if (std::is_same_v) + return 1e-5f; + else if (std::is_same_v) + 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(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::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) { + // Q: Why only double? + // A: The numeric values (e.g. overflow threshold `n`) below are different 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::infinity(); + assert(std::hermite(n, +300.0) == inf); + assert(std::hermite(n, -300.0) == ((n & 1) ? -inf : inf)); + } + } + } + } +} + +struct TestFloat { + template + void operator()() { + test(); + } +}; + +struct TestInt { + template + 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(x))); + } +}; + +int main() { + types::for_each(types::floating_point_types(), TestFloat()); + types::for_each(types::type_list(), 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