Added implementation of expint [sf.cmath]

[thebrandre]

Summary

Added implementations of std::expint to cmath as specified in "29.9.5 Mathematical special functions [sf.cmath]" of the C++17 standard.

The implementation is based on Press, William H., et al. Numerical recipes 3rd edition: The art of scientific computing. Cambridge university press, 2007.

It is tested with a table of values generated by gcc-8.3.1 and boost-1.66.0 (confirmed with msvc-192027508 and boost-1.69.0) using boost::multiprecision::cpp_bin_float_oct (256bit float).

There are also tests for analytical bounds for large arguments and arguments close to the singularity.
And we also check the corner cases when +/- infinity or 0 is passed as an argument.

The implementation is inline for the same reason as mentioned in D60869.

Related revisions:

• D60869: Added std::assoc_laguerre and std::laguerre
• D60234: Added std::assoc_legendre and std::legendre

[thebrandre]

• Moved expint from std::experimental to std -- from my understanding, that's where it belongs from C++17 onwards.
• Added _NOEXCEPT because it doesn't throw and all the other functions in math.h have it, too.
• Added the missing overloads for integral types.

[thebrandre]

• Fixed shared linking on Windows : the internal implementation libcpp_expint<T> is no longer exported via explicit instantiation of the template. Instead, there are now three functions libcpp_expintf, libcpp_expint, and libcpp_expintl annotated with _LIBCPP_FUNC_VIS. I also put them into namespace std.
• Added "UNSUPPORTED: c++98, c++03, c++11, c++14" to the unit tests because std::expint is ifdef'ed out for versions prior to c++17.

Open questions:

• The internal implementations libcpp_expintf, libcpp_expint, and __libcpp_expintl are currently exported even if _LIBCPP_STD_VER is 14 or less because the declarations with _LIBCPP_FUNC_VIS are outside the #if _LIBCPP_STD_VER > 14 ... so one might link the same binaries for all values of _LIBCPP_STD_VER. I am not quite sure if this is desired or aspired ...
• The internal implementations could be put inside extern "C". But there is no actual benefit in doing so because the rest of the library doesn't do that either. Is that correct?

Sorry for all the corrections. I am new to the project ... but I have some hope that I finally got the hang of it. ;-)

[thebrandre]

Updated the diff: no longer make changes to math.h (as specified in ISO/IEC JTC 1/SC 22/WG 21 N3060) but only to cmath (as specified in "29.9.5 Mathematical special functions [sf.cmath]") -- updated the summary accordingly.

[thebrandre]

A very small change in the synopsis in cmath: expint is now listed separately in a new category ... just as in 29.9.1 [cmath.syn].
This is also for consistency with https://reviews.llvm.org/D60234

[thebrandre]

• The implementation is inline for the same reason as mentioned in D60869.
• Moved the tests for consistency with the latest changes in cmath (lerp).
• Simplification of the tests: now using only one set of reference values generated by boost with multiprecision floats (256bit rounded to 128bit). I am pretty optimistic that every digit is correct, so that the comparison to other implementations no longer serves any purpose.
• Rebased against current master branch.

libcxx/include/cmath

 float          truncf(float x);
 long double    truncl(long double x);
 
+// C++17 Mathematical special functions [sf.cmath]
+
+floating_point expint (arithmetic x);
+float          expintf(float x);
+long double    expintl(long double x);
+
 }  // std
 
 */

 #include <__config>
 #include <math.h>
 #include <version>
+#include <cstdlib>
+#include <limits>
 
 #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
 #pragma GCC system_header

     return isfinite(__lcpp_x);
 }
 
+#if _LIBCPP_STD_VER > 14
+
+// expint
+
+/// \brief Returns the Euler-Mascheroni constant.
+///
+/// The value is taken from:
+/// Sweeney, Dura W. "On the computation of Euler’s constant."
+/// Mathematics of Computation 17.82 (1963): 170-178.
+/// \return the Euler-Mascheroni constant \f$ \gamma = 0.5772156649\ldots \f$
+template <class _Real>
+inline _LIBCPP_INLINE_VISIBILITY _Real __libcpp_euler_gamma() _NOEXCEPT {
+  return _Real(0.57721566490153286060651209008240243104215933593992l);
+}
+
+/// \brief Evaluates std::expint(x) for x < -1.
+///
+/// This function evaluates the following expression:
+/// \f[
+/// \mathrm{Ei} (x) = -e^{x} \left(
+/// \frac{1}{x+1-}\frac{1}{x+3-}\frac{4}{x+5-} \cdots\right)
+/// \f]
+/// See Press, William H., et al. Numerical recipes 3rd edition: The art of
+/// scientific computing. Cambridge university press, 2007, p. 267., eq. (6.3.5)
+/// The implementation uses two iterations of the form
+/// \f[
+/// s^{(i)} = 1 - x + 2i - \frac{i^2}{s^{(i-1)}}
+/// \f]
+/// with different start values \f$ s_0^{(1)} = 1 -x + 2 \f$ and
+/// \f$ s_1^{(1)} = 1 -x + 2 - \frac{1}{1-x} \f$.
+/// Both of them are bound by
+/// \f$ 1 - x +i < s^{(i)} \leq 1 - x + 2i \f$.
+/// Therefore, the difference monotonically decreasing
+/// \f[
+/// s_0^{(i)} - s_1^{(i)} = \frac{i^2}{s_0^{(i-1)} s_1^{(i-1)}}
+/// \left( s_0^{(i)} - s_1^{(i)} \right).
+/// \f]
+/// \param [in] __x Argument \f$ x < -1\f$ of \f$\mathrm{Ei} (x) \f$.
+/// \return continued fraction representation of \f$ \mathrm{Ei} (x) \f$
+/// for \f$ x < -1 \f$
+template <class _Real>
+inline _LIBCPP_INLINE_VISIBILITY _Real
+__libcpp_expint_continued_fraction(_Real __x) _NOEXCEPT {
+  const int __max_iter = 1000;
+  const _Real __mxp1 = -__x + _Real(1);
+  if (_VSTD::isinf(__mxp1 + _Real(2 * __max_iter)))
+    return _Real(0);
+
+  // See algorithm expint.h on p. 268 of
+  // Press, William H., et al. Numerical recipes 3rd edition: The art of
+  // scientific computing. Cambridge university press, 2007.
+  // a -> -__t0
+  // b -> __t1
+  // c -> _s0
+  // d -> 1/__ s1
+
+  _Real __h = _Real(1) / __mxp1;
+  _Real __s0 = __mxp1 + _Real(2);
+  _Real __s1 = __s0 - __h;
+  __h *= __s0 / __s1;
+
+  for (int __i = 2; __i < __max_iter; ++__i) {
+    const _Real __t0 = _Real(__i * __i);
+    const _Real __t1 = __mxp1 + _Real(2 * __i);
+    __s0 = __t1 - __t0 / __s0;
+    __s1 = __t1 - __t0 / __s1;
+    if (__s0 == __s1)
+      return -__h * exp(__x);
+
+    __h *= __s0 / __s1;
+  }
+
+  // (See the estimates in the comment above.)
+  // The closer __x is to 0, the slower the convergence.
+  // We use this function only for __x < -1, which makes
+  // __x == -1 the worst case.
+
+  _LIBCPP_UNREACHABLE();
+}
+
+/// \brief Evaluates std::expint(x) for x != 0.
+///
+/// This function evaluates the following expression:
+/// \f[
+/// \mathrm{Ei} (x) = \gamma + \ln \left|x\right| + \sum \limits_{n=1}^{\infty}
+/// \frac{x^n}{n\cdotn!}
+/// \f]
+/// See Press, William H., et al. Numerical recipes 3rd edition: The art of
+/// scientific computing. Cambridge university press, 2007, p. 267., eq. (6.3.6)
+/// or
+/// Harris, Frank E. "Tables of the exponential integral Ei (x)."
+/// Mathematical Tables and Other Aids to Computation 11.57 (1957): 9-16.
+/// \attention This representation converges slowly for large
+/// \f$ \left|x\right| \f$!
+/// \param [in] __x Argument \f$ x \neq 0 \f$ of \f$\mathrm{Ei} (x) \f$.
+/// \return Evaluation of a series representation of \f$ \mathrm{Ei} (x) \f$
+/// for \f$ -1 < x < 1 \f$.
+template <class _Real>
+inline _LIBCPP_INLINE_VISIBILITY
+_Real __libcpp_expint_power_series(_Real __x) _NOEXCEPT {
+  _Real __sum = _Real(0);
+  _Real __k = _Real(1);
+  _Real factor = _Real(1);
+  for (;;) {
+    factor *= __x / __k;
+    const _Real __prev = __sum;
+    __sum += factor / __k;
+    if (__prev == __sum)
+      return __sum + __libcpp_euler_gamma<_Real>() +
+             _VSTD::log(_VSTD::abs(__x));
+
+    __k += _Real(1);
+  }
+}
+
+/// \brief Evaluates std::expint(x) for x > -std::log(
+/// std::numeric_limits<T>::epsilon()).
+///
+/// This function evaluates the following expression:
+/// \f[
+/// \mathrm{Ei} (x) \approx \frac{\exp (x)}{x} \left(
+/// \sum \limits_{k=0}^{n-1} \frac{k!}{x^k} + R_{n} (x) \right).
+/// \f]
+/// where \f$ n = \lceil x \rceil - 1\f$.
+/// The remainder \f$ R_{n} (x) \f$ is bound by
+/// \f$ 0 \leq R_{n} (x) \leq \frac{n!}{n^n} \f$.
+/// Since all the summands are positive and the sum is greater than 1,
+/// one can introduce a cutoff condition at \f$ R_{n} (x) < \epsilon \f$.
+/// Using Stirling's approximation, one can see that this condition fulfilled
+/// for \f$ n > - \log \epsilon \f$
+/// \f[
+/// \frac{ R_{n} (x)}{\mathrm{Ei} (x)} \geq .
+/// \f]
+/// See eqs (2.3) and (2.4) from:
+/// Cody, W., & Thacher, H. (1969).
+/// Chebyshev Approximations for the Exponential Integral Ei(x).
+/// Mathematics of Computation, 23(106), 289-303. doi:10.2307/2004423
+/// or eq. (6.3.11) from:
+/// Press, William H., et al. Numerical recipes 3rd edition: The art of
+/// scientific computing. Cambridge university press, 2007, p. 269.
+/// \tparam _Real float, double, or long double
+/// \param __x Argument \f$ x \gg 1 \f$ of \f$\mathrm{Ei} (x) \f$.
+/// \return Evaluation of an asymptotic series representation of
+/// \f$ \mathrm{Ei} (x) \f$ for \f$ x \gg 1 \f$.
+template <class _Real>
+inline _LIBCPP_INLINE_VISIBILITY
+_Real __libcpp_expint_asymptotic(_Real __x) _NOEXCEPT {
+  _Real __sum = _Real(0);
+  _Real __t = _Real(1);
+  _Real __k = _Real(1);
+  do {
+    __t *= __k / __x;
+    __sum += __t;
+    __k += _Real(1);
+  } while (__k < __x && __t > _VSTD::numeric_limits<_Real>::epsilon());
+
+  const _Real __overflow_limit =
+      (_Real(1) - _VSTD::numeric_limits<_Real>::epsilon()) *
+      _VSTD::log(_VSTD::numeric_limits<_Real>::max());
+
+  if (__x >= __overflow_limit) {
+    if (_VSTD::isinf(__x)) {
+      return _VSTD::numeric_limits<_Real>::infinity();
+    }
+    // The first term is actually exp(x)/x in disguise.
+    // It is less precise but isn't subject to premature overflow.
+    return _VSTD::exp(__x - _VSTD::log(__x)) * (_Real(1) + __sum);
+  }
+
+  return _VSTD::exp(__x) / __x * (_Real(1) + __sum);
+}
+
+/// \brief Internal implementation of std::expint.
+///
+/// The implementation is based on chapter 6.3, p. 266 et seqq. of
+/// Press, William H., et al. Numerical recipes 3rd edition: The art of
+/// scientific computing. Cambridge university press, 2007
+/// \tparam  _Real float, double, or long double
+/// \param __x Argument \f$ x \gg 1 \f$ of \f$\mathrm{Ei} (x) \f$.
+/// \return std::expint(x).
+template <class _Real>
+inline _LIBCPP_INLINE_VISIBILITY
+_Real __libcpp_expint(_Real __x) _NOEXCEPT {
+  if (__x < _Real(-1))
+    return __libcpp_expint_continued_fraction(__x);
+
+  const _Real __asymptotic =
+      -_VSTD::log(_VSTD::numeric_limits<_Real>::epsilon());
+  if (__x < __asymptotic)
+    return __libcpp_expint_power_series(__x);
+
+  return __libcpp_expint_asymptotic(__x);
+}
+
+template <class _A1>
+inline _LIBCPP_INLINE_VISIBILITY
+typename __lazy_enable_if<is_arithmetic<_A1>::value, __promote<_A1>>::type
+expint(_A1 __lcpp_x) {
+  typedef typename __promote<_A1>::type __result_type;
+  return __libcpp_expint<__result_type>((__result_type)__lcpp_x);
+}
+
+inline _LIBCPP_INLINE_VISIBILITY
+float expintf(float __lcpp_x) {
+  return __libcpp_expint<float>(__lcpp_x);
+}
+
+inline _LIBCPP_INLINE_VISIBILITY
+long double expintl(long double __lcpp_x) {
+  return __libcpp_expint<long double>(__lcpp_x);
+}
+
+#endif // _LIBCPP_STD_VER > 14
+
+
 #if _LIBCPP_STD_VER > 17
 template <typename _Fp>
 constexpr

libcxx/test/std/numerics/c.math/c.math.expint/expint_bounds.pass.cpp

+//===----------------------------------------------------------------------===//
+//
+// 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++98, c++03, c++11, c++14
+
+// <cmath>
+
+// floating_point expint (arithmetic x);  // C++17
+// float          expintf(float x);       // C++17
+// long double    expintl(long double x); // C++17
+
+#include <cmath>
+#include <cassert>
+#include <initializer_list>
+
+template <class T>
+bool isInsideInterval(T Value, T LowerBound, T UpperBound, T RelTolerance,
+                      T AbsTolerance) {
+  // avoid inf times 0 because RelTolerance might be 0.
+  const T LowerBoundWithTol =
+      std::isinf(LowerBound)
+          ? LowerBound
+          : (LowerBound - AbsTolerance - std::abs(LowerBound) * RelTolerance);
+  const T UpperBoundWithTol =
+      std::isinf(UpperBound)
+          ? UpperBound
+          : (UpperBound + AbsTolerance + std::abs(UpperBound) * RelTolerance);
+  return Value >= LowerBoundWithTol && Value <= UpperBoundWithTol;
+}
+
+/// \brief Check the values of expint(x) - expint(-x) for x>0 against
+/// analytical upper and lower bounds.
+///
+/// \note This test is based on the formula
+/// \f[
+/// \mathrm{Ei}(x) - \mathrm{Ei}(-x) =
+/// 2 \int \limits_{0}^{x} \frac{\sinh(t)}{t} \,\mathrm{dt}.
+/// \f]
+/// One immediately gets the following relations
+/// \f[
+/// \Rightarrow \mathrm{Ei} (x) - \mathrm{Ei} (-x) \leq 2 \sinh (x)
+/// \f]
+/// and
+/// \f[
+/// \Rightarrow \mathrm{Ei} (x) - \mathrm{Ei} (-x) \geq
+/// 2 \sum \limits_{n=0}^{N} \frac{x^{2n+1}}{(2n+1)!(2n+1)}
+/// \quad \forall N \geq 0.
+/// \f]
+/// The test is confined to \f$ N = 5 \f$.
+/// \param [in] x x and -x are used as arguments of expint, x must be > 0.
+/// \param [in] RelTolerance relative tolerance for the value of
+/// expint(x) - expint(-x)
+/// \param [in] AbsTolerance absolute tolerance for the value of
+/// expint(x) - expint(-x)
+template <class T>
+void testSingularity(T x, T RelTolerance, T AbsTolerance) {
+  // Precondition of the test.
+  assert(!std::isinf(std::expint(x)));
+  assert(!std::isinf(std::expint(-x)));
+  const T Difference = std::expint(x) - std::expint(-x);
+  const T LowerBound = x*(T(2) + x * x * (T(1)/ T(9) +  x * x / T(300)));
+  const T UpperBound = T(2) * std::sinh(x);
+  const T AbsToleranceForDifference =
+      T(2) * AbsTolerance +
+      (std::abs(std::expint(x)) + std::abs(std::expint(-x))) * RelTolerance;
+  assert(isInsideInterval(Difference, LowerBound, UpperBound, T(0),
+                          AbsToleranceForDifference));
+}
+
+template <class T>
+void testSingularity(T RelTolerance, T AbsTolerance) {
+  for (T x : {T(1e-3), T(1e-2), T(0.1), T(1), T(10)})
+    testSingularity<T>(x, RelTolerance, AbsTolerance);
+}
+
+/// \brief Check the values of expint for x<0 against
+/// analytical upper and lower bounds.
+///
+/// This test is based on 5.1.20, p. 229 of
+/// Abramowitz, M., & Stegun, I. A. (1965).
+/// Handbook of mathematical functions
+/// with formulas, graphs, and mathematical table. In US Department of Commerce.
+/// National Bureau of Standards Applied Mathematics series 55.
+/// \param [in] x argument of expint, must be < 0.
+/// \param [in] RelTolerance relative tolerance for the value of expint(x)
+/// \param [in] AbsTolerance absolute tolerance for the value of expint(x)
+template <class T>
+void testBoundsNegativeX(T x, T RelTolerance, T AbsTolerance) {
+  assert(x < T(0));
+  const T LowerBound = -std::exp(x) * std::log(T(1) - T(1) / x);
+  const T UpperBound = -std::exp(x) * std::log(T(1) - T(2) / x) / T(2);
+  assert(isInsideInterval(std::expint(x), LowerBound, UpperBound, RelTolerance,
+                          AbsTolerance));
+}
+
+template <class T>
+void testBoundsNegativeX(T RelTolerance, T AbsTolerance) {
+  for (T x : {T(0.01), T(0.1), T(1), T(10), T(20), T(100), T(500)})
+    testBoundsNegativeX(-x, RelTolerance, AbsTolerance);
+
+  const T LogMin = std::log(std::numeric_limits<T>::min());
+  assert(std::expint(LogMin) < std::numeric_limits<T>::min());
+}
+
+/// \brief Checks the values of expint for x > 0 against
+/// analytical upper and lower bounds. The bounds are narrow for x < 1.
+///
+/// Derived from 5.1.40, p. 230 of
+/// Abramowitz, M., & Stegun, I. A. (1965).
+/// Handbook of mathematical functions
+/// with formulas, graphs, and mathematical table. In US Department of Commerce.
+/// for \f$ x > 0 \f$
+/// \f[ \mathrm{Ei} (x) - \Gamma - \ln x
+///         = \int \limits_{0}^{x} \frac{e^t-1}{t} \mathrm{d}t
+/// > \int \limits_{0}^{x} \left(1 + \tfrac{1}{2}t \right)\mathrm{d}t
+/// = x + \tfrac{1}{4}x^2 \f]
+/// \tparam T one of float, double, or long double
+/// \param [in] x argument of expint, must be > 0.
+/// \param [in] RelTolerance relative tolerance for the value of expint(x)
+/// \param [in] AbsTolerance absolute tolerance for the value of expint(x)
+template <class T>
+void testBoundsPositiveSmallX(T x, T RelTolerance, T AbsTolerance) {
+  assert(x > T(0));
+  const T Gamma = T(0.57721566490153286060651209008240243104215933593992l);
+  // works for all x > 0 but it is only a very rough estimate for x > 1
+  const T LowerBound = Gamma + std::log(x) + x + x * x / T(4);
+  const T UpperBound = Gamma - T(1) + std::log(x) + std::exp(x);
+  assert(isInsideInterval(std::expint(x), LowerBound, UpperBound, RelTolerance,
+                          AbsTolerance));
+}
+
+template <class T>
+void testBoundsPositiveSmallX(T RelTolerance, T AbsTolerance) {
+  for (T x : {T(0.01), T(0.1), T(1), T(10), T(20), T(50)})
+    testBoundsPositiveSmallX(x, RelTolerance, AbsTolerance);
+}
+
+/// \brief Checks analytic bounds for x >= 1.
+///
+/// For \f$ x \geq 1 \f$, one can write \f$ \mathrm{Ei} (x) \f$ as
+/// \f[
+/// \mathrm{Ei} (x) = \gamma + \ln x +
+/// \int \limits_{0}^{x} \frac{e^t}{t} \left( 1- e^{-t} \right)\mathrm{d}t
+/// = \mathrm{Ei} \left(1\right) + \ln x +
+/// \int \limits_{1}^{x} \frac{e^t}{t} \left( 1- e^{-t} \right)\mathrm{d}t.
+/// \f]
+/// From that equation, one derives the lower bound
+/// \f[
+/// \mathrm{Ei} (x) \geq \mathrm{Ei} \left(1\right) + \ln x +
+/// \int \limits_{1}^{x} \frac{e^t}{t} \left(1- \tfrac{1}{t}\right)\mathrm{d}t
+/// = \mathrm{Ei} \left(1\right) + \ln x + \frac{e^x}{x} - e,
+/// \f]
+/// and the upper bound
+/// \f[
+/// \mathrm{Ei} (x) \leq \mathrm{Ei} \left(1\right) + \ln x +
+/// \int \limits_{1}^{x} e^t \left( 1- e^{-t} \right) \mathrm{d}t
+/// = \mathrm{Ei} \left(1\right) + \ln x + e^x - e + x - 1
+/// \f]
+/// \note Lower and upper bound are identical for x == 1.
+/// \tparam T one of float, double, or long double
+/// \param [in] x argument of expint, must be >= 1.
+/// \param [in] RelTolerance relative tolerance for the value of expint(x)
+/// \param [in] AbsTolerance absolute tolerance for the value of expint(x)
+template <class T>
+void testBoundsPositiveLargeX(T x, T RelTolerance, T AbsTolerance) {
+  assert(x > T(1));
+  // avoid inf times 0 because RelTolerance might be 0.
+  const T LowerBound =
+      std::expint(T(1)) + std::log(x) + std::exp(x) / x - std::exp(T(1));
+  const T UpperBound = std::expint(T(1)) + std::log(x) +
+                       (std::exp(x) - std::exp(T(1))) - (x - T(1));
+  assert(isInsideInterval(std::expint(x), LowerBound, UpperBound, RelTolerance,
+                          AbsTolerance));
+}
+
+template <class T>
+void testBoundsPositiveLargeX(T RelTolerance, T AbsTolerance) {
+  for (T x : {T(1.1), T(1.5), T(2), T(10), T(20), T(50)})
+    testBoundsPositiveLargeX(x, RelTolerance, AbsTolerance);
+
+  // that's where the implementation switches to the asymptotic representation
+  const T Asymptotic = -std::log(std::numeric_limits<T>::epsilon());
+  for (int i = 1; i < 5; ++i)
+    testBoundsPositiveLargeX(T(i) * Asymptotic, RelTolerance, AbsTolerance);
+
+  // Test these separately because the bounds overflow:
+  const T LogMax = std::log(std::numeric_limits<T>::max());
+  assert(!std::isinf(std::expint(LogMax)));
+  const T GivesInf = T(4) * LogMax * LogMax;
+  assert(std::isinf(std::expint(GivesInf)));
+}
+
+template <class T>
+void testAllBounds(T RelTolerance, T AbsTolerance) {
+  testSingularity(RelTolerance, AbsTolerance);
+  testBoundsNegativeX(RelTolerance, AbsTolerance);
+  testBoundsPositiveSmallX(RelTolerance, AbsTolerance);
+  testBoundsPositiveLargeX(RelTolerance, AbsTolerance);
+}
+
+int main(int, char**) {
+  testAllBounds<float>(1e-6f, 0.0f);
+  testAllBounds<double>(1e-14, 0.0);
+}

libcxx/test/std/numerics/c.math/c.math.expint/expint_comparisons.pass.cpp

+//===----------------------------------------------------------------------===//
+//
+// 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++98, c++03, c++11, c++14
+
+// <cmath>
+
+// floating_point expint (arithmetic x);  // C++17
+// float          expintf(float x);       // C++17
+// long double    expintl(long double x); // C++17
+
+#include <cmath>
+#include <cassert>
+#include <map>
+
+/// \brief Returns a map of x and std::expint(x)
+///
+/// This table is generated with gcc-8.3.1 and boost-1.66.0
+/// and confirmed with msvc-192027508 and boost-1.69.0
+/// using boost::multiprecision::cpp_bin_float_oct
+/// rounded to std::numeric_limits<float128>::max_digits10.
+/// \note The key_type is float to make sure that it can be
+/// exactly represented by all floating-point types.
+/// The mapped_type is long double to avoid precision loss
+/// but no overflow occurs when casting to float.
+std::map<float, long double> generateExpintTestValues()
+{
+  return {{-5.000000000e+01f, -3.783264029550459018698967854021285780e-24L},
+      {-4.900000000e+01f, -1.048981164236802369849688680874776641e-23L},
+      {-4.800000000e+01f, -2.909664190405842339165765290916506568e-23L},
+      {-4.700000000e+01f, -8.074197842725813951937477940562408284e-23L},
+      {-4.600000000e+01f, -2.241531759744299745983944746971850400e-22L},
+      {-4.500000000e+01f, -6.225690809462383643095308453348153894e-22L},
+      {-4.400000000e+01f, -1.729959874281647757227851510531351102e-21L},
+      {-4.300000000e+01f, -4.809496556950017850712335884632312421e-21L},
+      {-4.200000000e+01f, -1.337790881001177514314392868306569266e-20L},
+      {-4.100000000e+01f, -3.723166776459977716235516016523522486e-20L},
+      {-4.000000000e+01f, -1.036773261451656972150641889145259771e-19L},
+      {-3.900000000e+01f, -2.888779301522700998725095238399534096e-19L},
+      {-3.800000000e+01f, -8.054106914290749868944421587855629671e-19L},
+      {-3.700000000e+01f, -2.247020697588571222083353307904064249e-18L},
+      {-3.600000000e+01f, -6.273339009762241588152637033384756530e-18L},
+      {-3.500000000e+01f, -1.752705938994737200054832669092913667e-17L},
+      {-3.400000000e+01f, -4.900676118392787715370624194027081175e-17L},
+      {-3.300000000e+01f, -1.371384348448746572239655805859484591e-16L},
+      {-3.200000000e+01f, -3.840961801225066831461234630025750207e-16L},
+      {-3.100000000e+01f, -1.076767038616238261726105157668358207e-15L},
+      {-3.000000000e+01f, -3.021552010688812544815825045153697921e-15L},
+      {-2.900000000e+01f, -8.487759778353562840972762979415674942e-15L},
+      {-2.800000000e+01f, -2.386941511933733162909373835426501483e-14L},
+      {-2.700000000e+01f, -6.720637435262040377094883818453725297e-14L},
+      {-2.600000000e+01f, -1.894685885674978242549641305135035184e-13L},
+      {-2.500000000e+01f, -5.348899755340216640325471118085924169e-13L},
+      {-2.400000000e+01f, -1.512305893999705766474308468724763618e-12L},
+      {-2.300000000e+01f, -4.282684795665672618906770326270850700e-12L},
+      {-2.200000000e+01f, -1.214937895620437265840775447119305968e-11L},
+      {-2.100000000e+01f, -3.453201267146756266667883784577524543e-11L},
+      {-2.000000000e+01f, -9.835525290649881690396987108894776074e-11L},
+      {-1.900000000e+01f, -2.807829097060795267324152564574528010e-10L},
+      {-1.800000000e+01f, -8.036090344828677657207213012471781597e-10L},
+      {-1.700000000e+01f, -2.306431989821654494245004407041613170e-09L},
+      {-1.600000000e+01f, -6.640487249441042785696129359642739491e-09L},
+      {-1.500000000e+01f, -1.918627892147866977104309696741338816e-08L},
+      {-1.400000000e+01f, -5.565631111145182114989920850651295593e-08L},
+      {-1.300000000e+01f, -1.621866218801432870669074367074383031e-07L},
+      {-1.200000000e+01f, -4.751081824672493932594612696661441836e-07L},
+      {-1.100000000e+01f, -1.400300304247441775446564237904285784e-06L},
+      {-1.000000000e+01f, -4.156968929685324277402859810278180384e-06L},
+      {-9.000000000e+00f, -1.244735417800627212114856523810080257e-05L},
+      {-8.000000000e+00f, -3.766562284392490177255799595075272089e-05L},
+      {-7.599999905e+00f, -5.885877835522211132219909198852529146e-05L},
+      {-7.199999809e+00f, -9.218813666493392496819998690716120347e-05L},
+      {-7.000000000e+00f, -1.154817316103382164310114560437894989e-04L},
+      {-6.800000191e+00f, -1.447577610043514442473650415536007813e-04L},
+      {-6.400000095e+00f, -2.279479311702414515081458172733752388e-04L},
+      {-6.000000000e+00f, -3.600824521626586592953941157717971888e-04L},
+      {-5.599999905e+00f, -5.708402326224589986646571434159349354e-04L},
+      {-5.199999809e+00f, -9.086218147950596268246438778225680700e-04L},
+      {-5.000000000e+00f, -1.148295591275325797330561969819722076e-03L},
+      {-4.800000191e+00f, -1.452993609857595605868097990432304282e-03L},
+      {-4.400000095e+00f, -2.336009774551668672538799012916084760e-03L},
+      {-4.000000000e+00f, -3.779352409848906478874860132466414852e-03L},
+      {-3.000000000e+00f, -1.304838109419703741250074582864502295e-02L},
+      {-2.000000000e+00f, -4.890051070806111956723983522804952231e-02L},
+      {-1.000000000e+00f, -2.193839343955202736771637754601216490e-01L},
+      {-8.000000119e-01f, -3.105965718500198010734802054347857172e-01L},
+      {-6.000000238e-01f, -4.543794813815877227508100153022577629e-01L},
+      {-4.000000060e-01f, -7.023801088771155253272499368438163604e-01L},
+      {-2.000000030e-01f, -1.222650531983854271436328625153714048e+00L},
+      {2.000000030e-01f, -8.217605697020810421808504320544601028e-01L},
+      {4.000000060e-01f, 1.047652408492449863690536252242036180e-01L},
+      {6.000000238e-01f, 7.698813623418547522320476636656036129e-01L},
+      {8.000000119e-01f, 1.347396581375470007673701978355742968e+00L},
+      {1.000000000e+00f, 1.895117816355936755466520934331634269e+00L},
+      {2.000000000e+00f, 4.954234356001890163379505130227035276e+00L},
+      {3.000000000e+00f, 9.933832570625416558008336019216765263e+00L},
+      {4.000000000e+00f, 1.963087447005622002264572027972383885e+01L},
+      {4.400000095e+00f, 2.600897503700524499824399694364832107e+01L},
+      {4.800000191e+00f, 3.469789470216587397588224274169374529e+01L},
+      {5.000000000e+00f, 4.018527535580317745509142179379586710e+01L},
+      {5.199999809e+00f, 4.662484385677105401583505624620463816e+01L},
+      {5.599999905e+00f, 6.310178136896516972405815088748916221e+01L},
+      {6.000000000e+00f, 8.598976214243920480358340030799069004e+01L},
+      {6.400000095e+00f, 1.179348746683721712003557710511492152e+02L},
+      {6.800000191e+00f, 1.627071127553717604838712010833682320e+02L},
+      {7.000000000e+00f, 1.915047433355013959530631482724569474e+02L},
+      {7.199999809e+00f, 2.256877722182689852472513019665927695e+02L},
+      {7.599999905e+00f, 3.145718534240306126656270737853903947e+02L},
+      {8.000000000e+00f, 4.403798995348382689974245966593933924e+02L},
+      {9.000000000e+00f, 1.037878290717089587657573226793622215e+03L},
+      {1.000000000e+01f, 2.492228976241877759138440143998524849e+03L},
+      {1.100000000e+01f, 6.071406374098611507964887284858515516e+03L},
+      {1.200000000e+01f, 1.495953266639752885229246187605753281e+04L},
+      {1.300000000e+01f, 3.719768849068903560439164528876347903e+04L},
+      {1.400000000e+01f, 9.319251363396537129882452836392441672e+04L},
+      {1.500000000e+01f, 2.349558524907683035782457458951611726e+05L},
+      {1.600000000e+01f, 5.955609986708370018501610068484626103e+05L},
+      {1.700000000e+01f, 1.516637894042516884432797432876246260e+06L},
+      {1.800000000e+01f, 3.877904330597443502996466079950798035e+06L},
+      {1.900000000e+01f, 9.950907251046844760026002538253063332e+06L},
+      {2.000000000e+01f, 2.561565266405658882048112080409807183e+07L},
+      {2.100000000e+01f, 6.612718635548492136250291987962721053e+07L},
+      {2.200000000e+01f, 1.711446713003636684975370635366578572e+08L},
+      {2.300000000e+01f, 4.439663698302712208698485234569711018e+08L},
+      {2.400000000e+01f, 1.154115391849182948286759099974094906e+09L},
+      {2.500000000e+01f, 3.005950906525548689841377604167345810e+09L},
+      {2.600000000e+01f, 7.842940991898186370453025612588634511e+09L},
+      {2.700000000e+01f, 2.049649711988081236484165294369597173e+10L},
+      {2.800000000e+01f, 5.364511859231469415605083102066723533e+10L},
+      {2.900000000e+01f, 1.405991957584069047339548366002198035e+11L},
+      {3.000000000e+01f, 3.689732094072741970640063289108457470e+11L},
+      {3.100000000e+01f, 9.694555759683939661661691071346617779e+11L},
+      {3.200000000e+01f, 2.550043566357786926146742307679350278e+12L},
+      {3.300000000e+01f, 6.714640184076497558707440521912031531e+12L},
+      {3.400000000e+01f, 1.769803724411626854310341944075653955e+13L},
+      {3.500000000e+01f, 4.669055014466159544500146290990063797e+13L},
+      {3.600000000e+01f, 1.232852079912097685430891907640127096e+14L},
+      {3.700000000e+01f, 3.257988998672263996790001683698134913e+14L},
+      {3.800000000e+01f, 8.616388199965786544948016613356301502e+14L},
+      {3.900000000e+01f, 2.280446200301902595340816714438768255e+15L},
+      {4.000000000e+01f, 6.039718263611241578359231418510691294e+15L},
+      {4.100000000e+01f, 1.600664914324504111069970545013283212e+16L},
+      {4.200000000e+01f, 4.244796092136850759367705586573591072e+16L},
+      {4.300000000e+01f, 1.126348290166966760275342380239567459e+17L},
+      {4.400000000e+01f, 2.990444718632336675058132672827320204e+17L},
+      {4.500000000e+01f, 7.943916035704453771510168303218332002e+17L},
+      {4.600000000e+01f, 2.111342388647824195000286962846296846e+18L},
+      {4.700000000e+01f, 5.614329680810343111535109717067692832e+18L},
+      {4.800000000e+01f, 1.493630213112993142255380874118611276e+19L},
+      {4.900000000e+01f, 3.975442747903744836006716995140844348e+19L},
+      {5.000000000e+01f, 1.058563689713169096306154143322998720e+20L}};
+}
+
+template <class T>
+void compareWithBoost(T RelTolerance, T AbsTolerance) {
+  const std::map<float, long double> Table = generateExpintTestValues();
+  for (const auto& xy : Table) {
+    const T Value = std::expint(static_cast<T>(xy.first));
+    const T ExpectedValue = static_cast<T>(xy.second);
+    const T Tolerance = std::abs(ExpectedValue) * RelTolerance + AbsTolerance;
+    const T Error = std::abs(Value - ExpectedValue);
+    assert(Error <= Tolerance);
+  }
+}
+
+/// \brief Checks for a change of sign close to the value from the literature
+/// for the root of expint.
+///
+/// \tparam T one of float, double, or long double
+/// \param [in] RelTolerance Tolerance for x; width of the interval where
+/// expint changes its sign.
+/// The reference value 0.372507... in its hexadecimal representation is taken
+/// from Cody, W., & Thacher, H. (1969).
+/// Chebyshev Approximations for the Exponential Integral Ei(x).
+/// Mathematics of Computation, 23(106), 289-303. doi:10.2307/2004423
+template <class T>
+void compareRootWithCody1969(T RelTolerance) {
+  const T Root = T(0x0.5F5CA54AD2D7F0F264C31p0L);
+  const T Tolerance = RelTolerance * Root;
+  assert(std::expint(Root - Tolerance) < T(0));
+  assert(std::expint(Root + Tolerance) > T(0));
+}
+
+int main(int, char**) {
+  compareRootWithCody1969<float>(1e-6f);
+  compareRootWithCody1969<double>(1e-14);
+
+  compareWithBoost<float>(1e-6f, 0.0f);
+  compareWithBoost<double>(1e-14, 0.0);
+}

libcxx/test/std/numerics/c.math/c.math.expint/expint_limits.pass.cpp

+//===----------------------------------------------------------------------===//
+//
+// 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++98, c++03, c++11, c++14
+
+// <cmath>
+
+// floating_point expint (arithmetic x);  // C++17
+// float          expintf(float x);       // C++17
+// long double    expintl(long double x); // C++17
+
+#include <cmath>
+#include <cassert>
+#include <limits>
+
+template <class T>
+void testLimits() {
+  const T NegativeLimit = std::expint(-std::numeric_limits<T>::infinity());
+  assert(NegativeLimit <= 0);
+  assert(NegativeLimit >= -std::numeric_limits<T>::min());
+
+  const T PositiveLimit = std::expint(std::numeric_limits<T>::infinity());
+  assert(PositiveLimit > T(0));
+  assert(std::isinf(PositiveLimit));
+
+  const T ExpintAtZero = std::expint(T(0));
+  assert(ExpintAtZero < T(0));
+  assert(std::isinf(ExpintAtZero));
+}
+
+int main(int, char**) {
+  testLimits<float>();
+  testLimits<double>();
+  testLimits<long double>();
+}

libcxx/test/std/numerics/c.math/cmath.pass.cpp

 Ambiguous tgamma(Ambiguous){ return Ambiguous(); }
 Ambiguous trunc(Ambiguous){ return Ambiguous(); }
 
+#if _LIBCPP_STD_VER > 14
+Ambiguous expint(Ambiguous){ return Ambiguous(); }
+#endif
+
 template <class T, class = decltype(std::abs(std::declval<T>()))>
 std::true_type has_abs_imp(int);
 template <class T>

     assert(std::trunc(1) == 1);
 }
 
+#if _LIBCPP_STD_VER > 14
+
+void test_expint()
+{
+  static_assert((std::is_same<decltype(std::expint((float)0)), float>::value), "");
+  static_assert((std::is_same<decltype(std::expint((bool)0)), double>::value), "");
+  static_assert((std::is_same<decltype(std::expint((unsigned short)0)), double>::value), "");
+  static_assert((std::is_same<decltype(std::expint((int)0)), double>::value), "");
+  static_assert((std::is_same<decltype(std::expint((unsigned int)0)), double>::value), "");
+  static_assert((std::is_same<decltype(std::expint((long)0)), double>::value), "");
+  static_assert((std::is_same<decltype(std::expint((unsigned long)0)), double>::value), "");
+  static_assert((std::is_same<decltype(std::expint((long long)0)), double>::value), "");
+  static_assert((std::is_same<decltype(std::expint((unsigned long long)0)), double>::value), "");
+  static_assert((std::is_same<decltype(std::expint((double)0)), double>::value), "");
+  static_assert((std::is_same<decltype(std::expint((long double)0)), long double>::value), "");
+  static_assert((std::is_same<decltype(std::expintf(0)), float>::value), "");
+  static_assert((std::is_same<decltype(std::expintl(0)), long double>::value), "");
+  static_assert((std::is_same<decltype(expint(Ambiguous())), Ambiguous>::value), "");
+}
+
+#endif
+
 int main(int, char**)
 {
     test_abs();

     test_tgamma();
     test_trunc();
 
+#if _LIBCPP_STD_VER > 14
+    test_expint();
+#endif
+
   return 0;
 }