[libcxx-commits] [libcxx] [libc++][math] Mathematical special functions: Implementing `std::laguerre`, `std::assoc_laguerre` (PR #92301)
via libcxx-commits
libcxx-commits at lists.llvm.org
Wed May 15 11:29:00 PDT 2024
https://github.com/PaulXiCao updated https://github.com/llvm/llvm-project/pull/92301
>From cc0b8d5d868d2a8163fb744f1d628fed7c955599 Mon Sep 17 00:00:00 2001
From: Paul <{ID}+{username}@users.noreply.github.com>
Date: Wed, 15 May 2024 20:18:02 +0200
Subject: [PATCH 1/2] laguerre, assoc_laguerr: impl + tests
---
libcxx/include/CMakeLists.txt | 1 +
libcxx/include/__math/special_functions.h | 124 ++++++++++
libcxx/include/cmath | 17 ++
libcxx/modules/std/cmath.inc | 4 +-
.../numerics/c.math/assoc_laguerre.pass.cpp | 194 ++++++++++++++++
.../std/numerics/c.math/laguerre.pass.cpp | 217 ++++++++++++++++++
6 files changed, 556 insertions(+), 1 deletion(-)
create mode 100644 libcxx/include/__math/special_functions.h
create mode 100644 libcxx/test/std/numerics/c.math/assoc_laguerre.pass.cpp
create mode 100644 libcxx/test/std/numerics/c.math/laguerre.pass.cpp
diff --git a/libcxx/include/CMakeLists.txt b/libcxx/include/CMakeLists.txt
index 01e9c247560ca..02c46d87759d4 100644
--- a/libcxx/include/CMakeLists.txt
+++ b/libcxx/include/CMakeLists.txt
@@ -516,6 +516,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..a805089b56e6f
--- /dev/null
+++ b/libcxx/include/__math/special_functions.h
@@ -0,0 +1,124 @@
+// -*- 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/abs.h>
+#include <__math/copysign.h>
+#include <__math/traits.h>
+#include <__type_traits/enable_if.h>
+#include <__type_traits/is_integral.h>
+#include <limits>
+#include <stdexcept>
+
+#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 __assoc_laguerre(unsigned __n, unsigned __alpha, _Real __x) {
+ // The associated/generalized Laguerre polynomial L_n^\alpha(x).
+ // The implementation is based on the recurrence formula:
+ //
+ // (j+1) L_{j+1}^\alpha(x) = (-x + 2j + \alpha + 1) L_j^\alpha(x) - (j + \alpha) L_{j-1}^\alpha(x)
+ //
+ // 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;
+
+ if (__x < 0)
+ std::__throw_domain_error("Argument `x` of Laguerre function is out of range: `x >= 0`.");
+
+ _Real __L_0{1};
+ if (__n == 0)
+ return __L_0;
+
+ _Real __L_n_prev = __L_0;
+ _Real __L_n = 1 + __alpha - __x;
+ for (unsigned __i = 1; __i < __n; ++__i) {
+ _Real __L_n_next =
+ ((-__x + 2 * __i + __alpha + 1) * __L_n - (__i + __alpha) * __L_n_prev) / static_cast<_Real>(__i + 1);
+ __L_n_prev = __L_n;
+ __L_n = __L_n_next;
+ }
+
+ if (!__math::isfinite(__L_n)) {
+ // Overflow occured!
+ // Can only happen for $x >> 1$ as _Real is at least double, and $__n < 128$ can be assumed.
+ _Real __inf = std::numeric_limits<_Real>::infinity();
+ return (__n & 1) ? -__inf : __inf;
+ }
+
+ return __L_n;
+ // NOLINTEND(readability-identifier-naming)
+}
+
+template <class _Real>
+_LIBCPP_HIDE_FROM_ABI _Real __laguerre(unsigned __n, _Real __x) {
+ return std::__assoc_laguerre(__n, /*alpha=*/0, __x);
+}
+
+inline _LIBCPP_HIDE_FROM_ABI double laguerre(unsigned __n, double __x) { return std::__laguerre(__n, __x); }
+
+inline _LIBCPP_HIDE_FROM_ABI float laguerre(unsigned __n, float __x) {
+ // use double internally -- float is too prone to overflow!
+ return static_cast<float>(std::laguerre(__n, static_cast<double>(__x)));
+}
+
+inline _LIBCPP_HIDE_FROM_ABI long double laguerre(unsigned __n, long double __x) { return std::__laguerre(__n, __x); }
+
+inline _LIBCPP_HIDE_FROM_ABI float laguerref(unsigned __n, float __x) { return std::laguerre(__n, __x); }
+
+inline _LIBCPP_HIDE_FROM_ABI long double laguerrel(unsigned __n, long double __x) { return std::laguerre(__n, __x); }
+
+template <class _Integer, std::enable_if_t<std::is_integral_v<_Integer>, int> = 0>
+_LIBCPP_HIDE_FROM_ABI double laguerre(unsigned __n, _Integer __x) {
+ return std::laguerre(__n, static_cast<double>(__x));
+}
+
+inline _LIBCPP_HIDE_FROM_ABI double assoc_laguerre(unsigned __n, unsigned __m, double __x) {
+ return std::__assoc_laguerre(__n, __m, __x);
+}
+
+inline _LIBCPP_HIDE_FROM_ABI float assoc_laguerre(unsigned __n, unsigned __m, float __x) {
+ // use double internally -- float is too prone to overflow!
+ return static_cast<float>(std::assoc_laguerre(__n, __m, static_cast<double>(__x)));
+}
+
+inline _LIBCPP_HIDE_FROM_ABI long double assoc_laguerre(unsigned __n, unsigned __m, long double __x) {
+ return std::__assoc_laguerre(__n, __m, __x);
+}
+
+inline _LIBCPP_HIDE_FROM_ABI float assoc_laguerref(unsigned __n, unsigned __m, float __x) {
+ return std::assoc_laguerre(__n, __m, __x);
+}
+
+inline _LIBCPP_HIDE_FROM_ABI long double assoc_laguerrel(unsigned __n, unsigned __m, long double __x) {
+ return std::assoc_laguerre(__n, __m, __x);
+}
+
+template <class _Integer, std::enable_if_t<std::is_integral_v<_Integer>, int> = 0>
+_LIBCPP_HIDE_FROM_ABI double assoc_laguerre(unsigned __n, unsigned __m, _Integer __x) {
+ return std::assoc_laguerre(__n, __m, 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 dd194bbb55896..7097c8ac88f5a 100644
--- a/libcxx/include/cmath
+++ b/libcxx/include/cmath
@@ -160,6 +160,14 @@ floating_point asinh (arithmetic x);
float asinhf(float x);
long double asinhl(long double x);
+double assoc_laguerre(unsigned n, unsigned m, double x); // C++17
+float assoc_laguerre(unsigned n, unsigned m, float x); // C++17
+long double assoc_laguerre(unsigned n, unsigned m, long double x); // C++17
+float assoc_laguerref(unsigned n, unsigned m, float x); // C++17
+long double assoc_laguerrel(unsigned n, unsigned m, long double x); // C++17
+template <class Integer>
+double assoc_laguerre(unsigned n, unsigned m, Integer x); // C++17
+
floating_point atanh (arithmetic x);
float atanhf(float x);
long double atanhl(long double x);
@@ -216,6 +224,14 @@ int ilogb (arithmetic x);
int ilogbf(float x);
int ilogbl(long double x);
+double laguerre(unsigned n, double x); // C++17
+float laguerre(unsigned n, float x); // C++17
+long double laguerre(unsigned n, long double x); // C++17
+float laguerref(unsigned n, float x); // C++17
+long double laguerrel(unsigned n, long double x); // C++17
+template <class Integer>
+double laguerre(unsigned n, Integer x); // C++17
+
floating_point lgamma (arithmetic x);
float lgammaf(float x);
long double lgammal(long double x);
@@ -315,6 +331,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/modules/std/cmath.inc b/libcxx/modules/std/cmath.inc
index a463c1e3ccf86..0ed30e8d1eaa4 100644
--- a/libcxx/modules/std/cmath.inc
+++ b/libcxx/modules/std/cmath.inc
@@ -263,12 +263,12 @@ export namespace std {
using std::signbit _LIBCPP_USING_IF_EXISTS;
// [sf.cmath], mathematical special functions
-#if 0
// [sf.cmath.assoc.laguerre], associated Laguerre polynomials
using std::assoc_laguerre;
using std::assoc_laguerref;
using std::assoc_laguerrel;
+#if 0
// [sf.cmath.assoc.legendre], associated Legendre functions
using std::assoc_legendre;
using std::assoc_legendref;
@@ -339,12 +339,14 @@ export namespace std {
using std::hermite;
using std::hermitef;
using std::hermitel;
+#endif
// [sf.cmath.laguerre], Laguerre polynomials
using std::laguerre;
using std::laguerref;
using std::laguerrel;
+#if 0
// [sf.cmath.legendre], Legendre polynomials
using std::legendre;
using std::legendref;
diff --git a/libcxx/test/std/numerics/c.math/assoc_laguerre.pass.cpp b/libcxx/test/std/numerics/c.math/assoc_laguerre.pass.cpp
new file mode 100644
index 0000000000000..268630f730da9
--- /dev/null
+++ b/libcxx/test/std/numerics/c.math/assoc_laguerre.pass.cpp
@@ -0,0 +1,194 @@
+//===----------------------------------------------------------------------===//
+//
+// 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 assoc_laguerre( unsigned n, unsigned m, double x);
+// float assoc_laguerre( unsigned n, unsigned m, float x);
+// long double assoc_laguerre( unsigned n, unsigned m, long double x);
+// float assoc_laguerref(unsigned n, unsigned m, float x);
+// long double assoc_laguerrel(unsigned n, unsigned m, long double x);
+// template <class Integer>
+// double assoc_laguerre( unsigned n, unsigned m, Integer x);
+
+#include <array>
+#include <cassert>
+#include <cmath>
+#include <limits>
+#include <numeric>
+
+#include "type_algorithms.h"
+
+inline constexpr unsigned g_max_n = 128;
+
+template <class Real>
+std::array<Real, 8> sample_points() {
+ return {0.0, 0.1, 0.5, 1.0, 2.7, 5.6, 10.3, 13.0};
+}
+
+template <class Real>
+class CompareFloatingValues {
+private:
+ Real tol;
+
+public:
+ CompareFloatingValues() {
+ if (std::is_same_v<Real, float>)
+ tol = 1e-4f;
+ else if (std::is_same_v<Real, double>)
+ tol = 1e-9;
+ else
+ tol = 1e-12l;
+ }
+
+ 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;
+
+ return std::abs(result - expected) < (tol + std::abs(expected) * tol);
+ }
+};
+
+template <class Real, std::size_t N>
+class Polynomial {
+ static_assert(N >= 1);
+ std::array<int, N> coeffs_;
+ int scaling_;
+
+public:
+ // Polynomial P = ( a0 + a1*x + a2*x^2 + ... + a_{N-1}*x^{N-1} ) / C
+ // coeffs = [a0, a1, .., a_{N-1}]
+ // scaling_ = C
+ Polynomial(Real, const int (&coeffs)[N], int scaling) : coeffs_{std::to_array(coeffs)}, scaling_{scaling} {}
+
+ // Evaluation at value x via Horner's method.
+ Real operator()(Real x) const {
+ return std::accumulate(coeffs_.rbegin() + 1,
+ coeffs_.rend(),
+ static_cast<Real>(coeffs_.back()),
+ [x](Real acc, int coeff) -> Real { return acc * x + coeff; }) /
+ scaling_;
+ };
+};
+
+template <class Real, std::size_t N>
+Polynomial(Real, const int (&)[N], int) -> Polynomial<Real, N>;
+
+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 (unsigned n = 0; n < g_max_n; ++n)
+ for (unsigned m = 0; m < g_max_n; ++m)
+ for (Real NaN : {nl::quiet_NaN(), nl::signaling_NaN()})
+ assert(std::isnan(std::assoc_laguerre(n, m, NaN)));
+ }
+
+ { // simple sample points for n, m = 0..127 should not produce NaNs.
+ for (unsigned n = 0; n < g_max_n; ++n)
+ for (unsigned m = 0; m < g_max_n; ++m)
+ for (Real x : sample_points<Real>())
+ assert(!std::isnan(std::assoc_laguerre(n, m, x)));
+ }
+
+ { // For any $x < 0$ a domain_error exception should be thrown.
+#ifndef _LIBCPP_NO_EXCEPTIONS
+ for (unsigned n = 0; n < g_max_n; ++n)
+ for (unsigned m = 0; m < g_max_n; ++m) {
+ const auto is_domain_error = [n, m](Real x) -> bool {
+ try {
+ std::assoc_laguerre(n, m, x);
+ } catch (const std::domain_error&) {
+ return true;
+ }
+ return false;
+ };
+
+ assert(is_domain_error(-std::numeric_limits<Real>::infinity()));
+ for (Real x : sample_points<Real>())
+ if (x > 0)
+ assert(is_domain_error(-x));
+ }
+#endif // _LIBCPP_NO_EXCEPTIONS
+ }
+
+ { // compare against std::laguerre
+ for (unsigned n = 0; n < g_max_n; ++n)
+ for (Real x : sample_points<Real>())
+ assert(std::assoc_laguerre(n, 0, x) == std::laguerre(n, x));
+ }
+
+ { // check against analytic polynoms for order n = 0..3, m = 0..127
+ for (int m = 0; m < static_cast<int>(g_max_n); ++m) {
+ for (Real x : sample_points<Real>()) {
+ const Polynomial p0{Real{}, {1}, 1};
+ const Polynomial p1{Real{}, {m + 1, -1}, 1};
+ const Polynomial p2{Real{}, {(m + 1) * (m + 2), -2 * (m + 2), 1}, 2};
+ const Polynomial p3{Real{}, {(m + 1) * (m + 2) * (m + 3), -3 * (m + 2) * (m + 3), 3 * (m + 3), -1}, 6};
+
+ const CompareFloatingValues<Real> compare;
+ assert(compare(std::assoc_laguerre(0, m, x), p0(x)));
+ assert(compare(std::assoc_laguerre(1, m, x), p1(x)));
+ assert(compare(std::assoc_laguerre(2, m, x), p2(x)));
+ assert(compare(std::assoc_laguerre(3, m, x), p3(x)));
+ }
+ }
+ }
+
+ { // checks std::assoc_laguerref for bitwise equality with std::assoc_laguerre(unsigned, float)
+ if constexpr (std::is_same_v<Real, float>)
+ for (unsigned n = 0; n < g_max_n; ++n)
+ for (unsigned m = 0; m < g_max_n; ++m)
+ for (float x : sample_points<float>())
+ assert(std::assoc_laguerre(n, m, x) == std::assoc_laguerref(n, m, x));
+ }
+
+ { // checks std::assoc_laguerrel for bitwise equality with std::assoc_laguerre(unsigned, long double)
+ if constexpr (std::is_same_v<Real, long double>)
+ for (unsigned n = 0; n < g_max_n; ++n)
+ for (unsigned m = 0; m < g_max_n; ++m)
+ for (long double x : sample_points<long double>())
+ assert(std::assoc_laguerre(n, m, x) == std::assoc_laguerrel(n, m, x));
+ }
+
+#if 0
+ { // evaluation at x=0: P_n(0) = 1
+ const CompareFloatingValues<Real> compare;
+ for (unsigned n = 0; n < g_max_n; ++n)
+ assert(compare(std::assoc_laguerre(n, Real{0}), 1));
+ }
+#endif
+}
+
+struct TestFloat {
+ template <class Real>
+ static void operator()() {
+ test<Real>();
+ }
+};
+
+struct TestInt {
+ template <class Integer>
+ static void operator()() {
+ // checks that std::assoc_laguerre(unsigned, Integer) actually wraps std::assoc_laguerre(unsigned, double)
+ for (unsigned n = 0; n < g_max_n; ++n)
+ for (unsigned m = 0; m < g_max_n; ++m)
+ for (Integer x{0}; x < 20; ++x)
+ assert(std::assoc_laguerre(n, m, x) == std::assoc_laguerre(n, m, 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/test/std/numerics/c.math/laguerre.pass.cpp b/libcxx/test/std/numerics/c.math/laguerre.pass.cpp
new file mode 100644
index 0000000000000..a1e16537b53c8
--- /dev/null
+++ b/libcxx/test/std/numerics/c.math/laguerre.pass.cpp
@@ -0,0 +1,217 @@
+//===----------------------------------------------------------------------===//
+//
+// 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 laguerre(unsigned n, double x);
+// float laguerre(unsigned n, float x);
+// long double laguerre(unsigned n, long double x);
+// float laguerref(unsigned n, float x);
+// long double laguerrel(unsigned n, long double x);
+// template <class Integer>
+// double laguerre(unsigned n, Integer x);
+
+#include <array>
+#include <cassert>
+#include <cmath>
+#include <limits>
+#include <numeric>
+
+#include "type_algorithms.h"
+
+inline constexpr unsigned g_max_n = 128;
+
+template <class Real>
+std::array<Real, 8> sample_points() {
+ return {0.0, 0.1, 0.5, 1.0, 2.7, 5.6, 10.3, 13.0};
+}
+
+template <class Real>
+class CompareFloatingValues {
+private:
+ Real tol;
+
+public:
+ CompareFloatingValues() {
+ if (std::is_same_v<Real, float>)
+ tol = 1e-4f;
+ else if (std::is_same_v<Real, double>)
+ tol = 1e-9;
+ else
+ tol = 1e-12l;
+ }
+
+ 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;
+
+ return std::abs(result - expected) < (tol + std::abs(expected) * tol);
+ }
+};
+
+template <class Real, std::size_t N>
+class Polynomial {
+ static_assert(N >= 1);
+ std::array<int, N> coeffs_;
+ int scaling_;
+
+public:
+ // Polynomial P = ( a0 + a1*x + a2*x^2 + ... + a_{N-1}*x^{N-1} ) / C
+ // coeffs = [a0, a1, .., a_{N-1}]
+ // scaling_ = C
+ Polynomial(Real, const int (&coeffs)[N], int scaling) : coeffs_{std::to_array(coeffs)}, scaling_{scaling} {}
+
+ // Evaluation at value x via Horner's method.
+ Real operator()(Real x) const {
+ return std::accumulate(coeffs_.rbegin() + 1,
+ coeffs_.rend(),
+ static_cast<Real>(coeffs_.back()),
+ [x](Real acc, int coeff) -> Real { return acc * x + coeff; }) /
+ scaling_;
+ };
+};
+
+template <class Real, std::size_t N>
+Polynomial(Real, const int (&)[N], int) -> Polynomial<Real, N>;
+
+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::laguerre(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::laguerre(n, x)));
+ }
+
+ { // For any $x < 0$ a domain_error exception should be thrown.
+#ifndef _LIBCPP_NO_EXCEPTIONS
+ for (unsigned n = 0; n < g_max_n; ++n) {
+ const auto is_domain_error = [n](Real x) -> bool {
+ try {
+ std::laguerre(n, x);
+ } catch (const std::domain_error&) {
+ return true;
+ }
+ return false;
+ };
+
+ assert(is_domain_error(-std::numeric_limits<Real>::infinity()));
+ for (Real x : sample_points<Real>())
+ if (x > 0)
+ assert(is_domain_error(-x));
+ }
+#endif // _LIBCPP_NO_EXCEPTIONS
+ }
+
+ { // check against analytic polynoms for order n=0..6
+ for (Real x : sample_points<Real>()) {
+ const Polynomial p0{Real{}, {1}, 1};
+ const Polynomial p1{Real{}, {1, -1}, 1};
+ const Polynomial p2{Real{}, {2, -4, 1}, 2};
+ const Polynomial p3{Real{}, {6, -18, 9, -1}, 6};
+ const Polynomial p4{Real{}, {24, -96, 72, -16, 1}, 24};
+ const Polynomial p5{Real{}, {120, -600, 600, -200, 25, -1}, 120};
+ const Polynomial p6{Real{}, {720, -4320, 5400, -2400, 450, -36, 1}, 720};
+
+ const CompareFloatingValues<Real> compare;
+ assert(compare(std::laguerre(0, x), p0(x)));
+ assert(compare(std::laguerre(1, x), p1(x)));
+ assert(compare(std::laguerre(2, x), p2(x)));
+ assert(compare(std::laguerre(3, x), p3(x)));
+ assert(compare(std::laguerre(4, x), p4(x)));
+ assert(compare(std::laguerre(5, x), p5(x)));
+ assert(compare(std::laguerre(6, x), p6(x)));
+ }
+ }
+
+ { // checks std::laguerref for bitwise equality with std::laguerre(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::laguerre(n, x) == std::laguerref(n, x));
+ }
+
+ { // checks std::laguerrel for bitwise equality with std::laguerre(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::laguerre(n, x) == std::laguerrel(n, x));
+ }
+
+ { // evaluation at x=0: P_n(0) = 1
+ const CompareFloatingValues<Real> compare;
+ for (unsigned n = 0; n < g_max_n; ++n)
+ assert(compare(std::laguerre(n, Real{0}), 1));
+ }
+
+ { // evaluation at x=+inf
+ Real inf = std::numeric_limits<Real>::infinity();
+ for (unsigned n = 1; n < g_max_n; ++n)
+ assert(std::laguerre(n, inf) == ((n & 1) ? -inf : inf));
+ }
+
+ { // check: if overflow occurs that it is mapped to the correct infinity
+ constexpr auto check_for_overflow = [](unsigned n_threshold, Real x) {
+ for (unsigned n = 0; n < g_max_n; ++n) {
+ if (n < n_threshold)
+ assert(std::isfinite(std::laguerre(n, x)));
+ else {
+ Real inf = std::numeric_limits<Real>::infinity();
+
+ // alternating limits (+-inf) only holds for x > largest root which is ~480.5 for order n=127.
+ assert(x > 481);
+ assert(std::laguerre(n, x) == ((n & 1) ? -inf : inf));
+ }
+ }
+ };
+
+ if constexpr (std::is_same_v<Real, float>) {
+ static_assert(sizeof(float) == 4);
+ check_for_overflow(23, 500.0f);
+ } else if constexpr (std::is_same_v<Real, double>) {
+ static_assert(sizeof(double) == 8);
+ check_for_overflow(116, 20'000.0);
+ } else if constexpr (std::is_same_v<Real, long double>) {
+ static_assert(sizeof(long double) == 16);
+ check_for_overflow(50, 1e100l);
+ }
+ }
+}
+
+struct TestFloat {
+ template <class Real>
+ static void operator()() {
+ test<Real>();
+ }
+};
+
+struct TestInt {
+ template <class Integer>
+ static void operator()() {
+ // checks that std::laguerre(unsigned, Integer) actually wraps std::laguerre(unsigned, double)
+ for (unsigned n = 0; n < g_max_n; ++n)
+ for (Integer x{0}; x < 20; ++x)
+ assert(std::laguerre(n, x) == std::laguerre(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());
+}
>From 9c28e06858f7901deb86a69ff90dff59d433685f Mon Sep 17 00:00:00 2001
From: Paul <{ID}+{username}@users.noreply.github.com>
Date: Wed, 15 May 2024 20:28:13 +0200
Subject: [PATCH 2/2] add header to modulemap
---
libcxx/include/module.modulemap | 1 +
1 file changed, 1 insertion(+)
diff --git a/libcxx/include/module.modulemap b/libcxx/include/module.modulemap
index 70dac2f19846b..034796fcb6732 100644
--- a/libcxx/include/module.modulemap
+++ b/libcxx/include/module.modulemap
@@ -1488,6 +1488,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" }
More information about the libcxx-commits
mailing list