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/experimental/cmath

+// -*- C++ -*-
+//===--------------------------- cmath ------------------------------------===//
+//
+// 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_EXPERIMENTAL_CMATH
+#define _LIBCPP_EXPERIMENTAL_CMATH
+
+/*
+   experimental/cmath synopsis
+
+#include <cmath>
+
+namespace std {
+namespace experimental {
+
+double      expint(double arg);
+float       expint(float arg);
+long double expint(long double arg);
+float       expintf(float arg);
+long double expintl(long double arg);
+
+} // namespace experimental
+} // namespace std
+
+*/
+
+#include <experimental/__config>
+#include <cmath>
+
+_LIBCPP_BEGIN_NAMESPACE_EXPERIMENTAL
+
+#if _LIBCPP_STD_VER > 14
+
+template <class _Real>
+_Real __libcpp_expint(_Real __x);
+
+inline _LIBCPP_INLINE_VISIBILITY float expintf(float __lcpp_x) {
+  return _VSTD_EXPERIMENTAL::__libcpp_expint<float>(__lcpp_x);
+}
+
+inline _LIBCPP_INLINE_VISIBILITY float expint(float __lcpp_x) {
+  return _VSTD_EXPERIMENTAL::__libcpp_expint<float>(__lcpp_x);
+}
+
+inline _LIBCPP_INLINE_VISIBILITY double expint(double __lcpp_x) {
+  return _VSTD_EXPERIMENTAL::__libcpp_expint<double>(__lcpp_x);
+}
+
+inline _LIBCPP_INLINE_VISIBILITY long double expintl(long double __lcpp_x) {
+  return _VSTD_EXPERIMENTAL::__libcpp_expint<long double>(__lcpp_x);
+}
+
+inline _LIBCPP_INLINE_VISIBILITY long double expint(long double __lcpp_x) {
+  return _VSTD_EXPERIMENTAL::__libcpp_expint<long double>(__lcpp_x);
+}
+
+#endif
+
+_LIBCPP_END_NAMESPACE_EXPERIMENTAL
+
+#endif // _LIBCPP_EXPERIMENTAL_CMATH

libcxx/src/experimental/expint.cpp

+//===------------------------ expint.cpp ------------------------*- 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
+//
+//===----------------------------------------------------------------------===//
+///
+/// \file
+/// This file contains the internal implementation of std::expint,
+/// __libcpp_expint<T>, and its explicit instantiations.
+///
+//===----------------------------------------------------------------------===//
+
+#include <experimental/cmath>
+#include <limits>
+#include <cstdlib>
+
+_LIBCPP_BEGIN_NAMESPACE_EXPERIMENTAL
+
+/// \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() {
+  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) {
+  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, we have to test for.
+
+  _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) {
+  _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) {
+  _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>
+_Real __libcpp_expint(_Real __x) {
+  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 float __libcpp_expint<float>(float);
+template double __libcpp_expint<double>(double);
+template long double __libcpp_expint<long double>(long double);
+
+_LIBCPP_END_NAMESPACE_EXPERIMENTAL

libcxx/test/std/experimental/numerics/c.math/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
+//
+//===----------------------------------------------------------------------===//
+
+// <cmath>
+
+// double      expint(double arg);
+// float       expint(float arg);
+// long double expint(long double arg);
+// float       expintf(float arg);
+// long double expintl(long double arg);
+
+#include <experimental/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::experimental::expint(x)));
+  assert(!std::isinf(std::experimental::expint(-x)));
+  const T Difference =
+      std::experimental::expint(x) - std::experimental::expint(-x);
+  const T LowerBound = T(2) * x + x * x * x / T(9) + x * x * x * x * x / T(300);
+  const T UpperBound = T(2) * std::sinh(x);
+  const T AbsToleranceForDifference =
+      T(2) * AbsTolerance + (std::abs(std::experimental::expint(x)) +
+                             std::abs(std::experimental::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::experimental::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::experimental::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::experimental::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::experimental::expint(T(1)) + std::log(x) +
+                       std::exp(x) / x - std::exp(T(1));
+  const T UpperBound = std::experimental::expint(T(1)) + std::log(x) +
+                       (std::exp(x) - std::exp(T(1))) - (x - T(1));
+  assert(isInsideInterval(std::experimental::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:
+  // boost's implementation raises overflow errors for these cases.
+  const T LogMax = std::log(std::numeric_limits<T>::max());
+  // fails with libstdc++
+  assert(!std::isinf(std::experimental::expint(LogMax)));
+  const T GivesInf = T(4) * LogMax * LogMax;
+  assert(std::isinf(std::experimental::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/experimental/numerics/c.math/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
+//
+//===----------------------------------------------------------------------===//
+
+// <cmath>
+
+// double      expint(double arg);
+// float       expint(float arg);
+// long double expint(long double arg);
+// float       expintf(float arg);
+// long double expintl(long double arg);
+
+#include <experimental/cmath>
+#include <cassert>
+
+/// \brief Calls the provided functor for pairs of x and std::expint(x)
+/// generated with boost 1.66.0 (long double version) and gcc 8.3.1 (linux)
+///
+/// \tparam BinaryFunction binary function with two long double arguments.
+/// \param ValueSink function that will be invoked for all pairs of
+/// x, std::expint(x).
+template <class BinaryFunction>
+void compareWithBoost(BinaryFunction&& ValueSink) {
+  ValueSink(-50.0L, -3.783264029550459018674e-24L);
+  ValueSink(-49.0L, -1.048981164236802369801e-23L);
+  ValueSink(-48.0L, -2.909664190405842339156e-23L);
+  ValueSink(-47.0L, -8.074197842725813951364e-23L);
+  ValueSink(-46.0L, -2.241531759744299745718e-22L);
+  ValueSink(-45.0L, -6.225690809462383642327e-22L);
+  ValueSink(-44.0L, -1.729959874281647757243e-21L);
+  ValueSink(-43.0L, -4.809496556950017851307e-21L);
+  ValueSink(-42.0L, -1.337790881001177514266e-20L);
+  ValueSink(-41.0L, -3.723166776459977716474e-20L);
+  ValueSink(-40.0L, -1.036773261451656972253e-19L);
+  ValueSink(-39.0L, -2.888779301522700998649e-19L);
+  ValueSink(-38.0L, -8.054106914290749868900e-19L);
+  ValueSink(-37.0L, -2.247020697588571222102e-18L);
+  ValueSink(-36.0L, -6.273339009762241587931e-18L);
+  ValueSink(-35.0L, -1.752705938994737200182e-17L);
+  ValueSink(-34.0L, -4.900676118392787715785e-17L);
+  ValueSink(-33.0L, -1.371384348448746572367e-16L);
+  ValueSink(-32.0L, -3.840961801225066831459e-16L);
+  ValueSink(-31.0L, -1.076767038616238261652e-15L);
+  ValueSink(-30.0L, -3.021552010688812544709e-15L);
+  ValueSink(-29.0L, -8.487759778353562840572e-15L);
+  ValueSink(-28.0L, -2.386941511933733162938e-14L);
+  ValueSink(-27.0L, -6.720637435262040377100e-14L);
+  ValueSink(-26.0L, -1.894685885674978242531e-13L);
+  ValueSink(-25.0L, -5.348899755340216640372e-13L);
+  ValueSink(-24.0L, -1.512305893999705766465e-12L);
+  ValueSink(-23.0L, -4.282684795665672618532e-12L);
+  ValueSink(-22.0L, -1.214937895620437265937e-11L);
+  ValueSink(-21.0L, -3.453201267146756266397e-11L);
+  ValueSink(-20.0L, -9.835525290649881690041e-11L);
+  ValueSink(-19.0L, -2.807829097060795267271e-10L);
+  ValueSink(-18.0L, -8.036090344828677657003e-10L);
+  ValueSink(-17.0L, -2.306431989821654494292e-09L);
+  ValueSink(-16.0L, -6.640487249441042786072e-09L);
+  ValueSink(-15.0L, -1.918627892147866977352e-08L);
+  ValueSink(-14.0L, -5.565631111145182114846e-08L);
+  ValueSink(-13.0L, -1.621866218801432870747e-07L);
+  ValueSink(-12.0L, -4.751081824672493932258e-07L);
+  ValueSink(-11.0L, -1.400300304247441775445e-06L);
+  ValueSink(-10.0L, -4.156968929685324277509e-06L);
+  ValueSink(-9.0L, -1.244735417800627212206e-05L);
+  ValueSink(-8.0L, -3.766562284392490177484e-05L);
+  ValueSink(-7.6L, -5.885877207538391849259e-05L);
+  ValueSink(-7.2L, -9.218811688716204236178e-05L);
+  ValueSink(-7.0L, -1.154817316103382164251e-04L);
+  ValueSink(-6.8L, -1.447577922448995782596e-04L);
+  ValueSink(-6.4L, -2.279479559293728467645e-04L);
+  ValueSink(-6.0L, -3.600824521626586592747e-04L);
+  ValueSink(-5.6L, -5.708401696482116743476e-04L);
+  ValueSink(-5.2L, -9.086216124486595848182e-04L);
+  ValueSink(-5.0L, -1.148295591275325797266e-03L);
+  ValueSink(-4.8L, -1.452993936878324219980e-03L);
+  ValueSink(-4.4L, -2.336010040655828993748e-03L);
+  ValueSink(-4.0L, -3.779352409848906478814e-03L);
+  ValueSink(-3.0L, -1.304838109419703741302e-02L);
+  ValueSink(-2.0L, -4.890051070806111956887e-02L);
+  ValueSink(-1.0L, -2.193839343955202737243e-01L);
+  ValueSink(-0.8L, -3.105965785455430345835e-01L);
+  ValueSink(-0.6L, -4.543795031894021081936e-01L);
+  ValueSink(-0.4L, -7.023801188656624785609e-01L);
+  ValueSink(-0.2L, -1.222650544183893088403e+00L);
+  ValueSink(0.2L, -8.217605879024003157421e-01L);
+  ValueSink(0.4L, 1.047652186193247932463e-01L);
+  ValueSink(0.6L, 7.698812899373594371543e-01L);
+  ValueSink(0.8L, 1.347396548212325938080e+00L);
+  ValueSink(1.0L, 1.895117816355936755501e+00L);
+  ValueSink(2.0L, 4.954234356001890163517e+00L);
+  ValueSink(3.0L, 9.933832570625416557783e+00L);
+  ValueSink(4.0L, 1.963087447005622002312e+01L);
+  ValueSink(4.4L, 2.600897327160514598335e+01L);
+  ValueSink(4.8L, 3.469788987377532763628e+01L);
+  ValueSink(5.0L, 4.018527535580317745431e+01L);
+  ValueSink(5.2L, 4.662485050579674776550e+01L);
+  ValueSink(5.6L, 6.310178597429926152310e+01L);
+  ValueSink(6.0L, 8.598976214243920485752e+01L);
+  ValueSink(6.4L, 1.179348657001818873394e+02L);
+  ValueSink(6.8L, 1.627070875714314148602e+02L);
+  ValueSink(7.0L, 1.915047433355013959538e+02L);
+  ValueSink(7.2L, 2.256878077010638103811e+02L);
+  ValueSink(7.6L, 3.145718784980834850373e+02L);
+  ValueSink(8.0L, 4.403798995348382689929e+02L);
+  ValueSink(9.0L, 1.037878290717089587614e+03L);
+  ValueSink(10.0L, 2.492228976241877759001e+03L);
+  ValueSink(11.0L, 6.071406374098611507861e+03L);
+  ValueSink(12.0L, 1.495953266639752885325e+04L);
+  ValueSink(13.0L, 3.719768849068903559996e+04L);
+  ValueSink(14.0L, 9.319251363396537130512e+04L);
+  ValueSink(15.0L, 2.349558524907683035963e+05L);
+  ValueSink(16.0L, 5.955609986708370018391e+05L);
+  ValueSink(17.0L, 1.516637894042516884497e+06L);
+  ValueSink(18.0L, 3.877904330597443502938e+06L);
+  ValueSink(19.0L, 9.950907251046844759912e+06L);
+  ValueSink(20.0L, 2.561565266405658882286e+07L);
+  ValueSink(21.0L, 6.612718635548492135422e+07L);
+  ValueSink(22.0L, 1.711446713003636685171e+08L);
+  ValueSink(23.0L, 4.439663698302712208824e+08L);
+  ValueSink(24.0L, 1.154115391849182948354e+09L);
+  ValueSink(25.0L, 3.005950906525548689999e+09L);
+  ValueSink(26.0L, 7.842940991898186370265e+09L);
+  ValueSink(27.0L, 2.049649711988081236556e+10L);
+  ValueSink(28.0L, 5.364511859231469415501e+10L);
+  ValueSink(29.0L, 1.405991957584069047421e+11L);
+  ValueSink(30.0L, 3.689732094072741970718e+11L);
+  ValueSink(31.0L, 9.694555759683939662576e+11L);
+  ValueSink(32.0L, 2.550043566357786926031e+12L);
+  ValueSink(33.0L, 6.714640184076497558594e+12L);
+  ValueSink(34.0L, 1.769803724411626854324e+13L);
+  ValueSink(35.0L, 4.669055014466159544373e+13L);
+  ValueSink(36.0L, 1.232852079912097685471e+14L);
+  ValueSink(37.0L, 3.257988998672263996887e+14L);
+  ValueSink(38.0L, 8.616388199965786545410e+14L);
+  ValueSink(39.0L, 2.280446200301902595459e+15L);
+  ValueSink(40.0L, 6.039718263611241578125e+15L);
+  ValueSink(41.0L, 1.600664914324504111523e+16L);
+  ValueSink(42.0L, 4.244796092136850758594e+16L);
+  ValueSink(43.0L, 1.126348290166966760156e+17L);
+  ValueSink(44.0L, 2.990444718632336675938e+17L);
+  ValueSink(45.0L, 7.943916035704453771250e+17L);
+  ValueSink(46.0L, 2.111342388647824194875e+18L);
+  ValueSink(47.0L, 5.614329680810343111000e+18L);
+  ValueSink(48.0L, 1.493630213112993142300e+19L);
+  ValueSink(49.0L, 3.975442747903744836000e+19L);
+  ValueSink(50.0L, 1.058563689713169096240e+20L);
+}
+
+/// \brief Calls the provided functor for pairs of x and std::expint(x)
+/// generated with gcc 8.3.1 (linux, long double version).
+///
+/// \tparam BinaryFunction binary function with two long double arguments.
+/// \param ValueSink function that will be invoked for all pairs of
+/// x, std::expint(x).
+template <class BinaryFunction>
+void compareWithGCC(BinaryFunction&& ValueSink) {
+  ValueSink(-50.0L, -3.783264029550459017956e-24L);
+  ValueSink(-49.0L, -1.048981164236802369658e-23L);
+  ValueSink(-48.0L, -2.909664190405842339156e-23L);
+  ValueSink(-47.0L, -8.074197842725813953659e-23L);
+  ValueSink(-46.0L, -2.241531759744299745948e-22L);
+  ValueSink(-45.0L, -6.225690809462383642786e-22L);
+  ValueSink(-44.0L, -1.729959874281647757059e-21L);
+  ValueSink(-43.0L, -4.809496556950017850940e-21L);
+  ValueSink(-42.0L, -1.337790881001177514486e-20L);
+  ValueSink(-41.0L, -3.723166776459977716180e-20L);
+  ValueSink(-40.0L, -1.036773261451656972253e-19L);
+  ValueSink(-39.0L, -2.888779301522700998649e-19L);
+  ValueSink(-38.0L, -8.054106914290749867959e-19L);
+  ValueSink(-37.0L, -2.247020697588571222102e-18L);
+  ValueSink(-36.0L, -6.273339009762241589059e-18L);
+  ValueSink(-35.0L, -1.752705938994737199881e-17L);
+  ValueSink(-34.0L, -4.900676118392787714882e-17L);
+  ValueSink(-33.0L, -1.371384348448746572247e-16L);
+  ValueSink(-32.0L, -3.840961801225066830496e-16L);
+  ValueSink(-31.0L, -1.076767038616238261652e-15L);
+  ValueSink(-30.0L, -3.021552010688812544901e-15L);
+  ValueSink(-29.0L, -8.487759778353562839801e-15L);
+  ValueSink(-28.0L, -2.386941511933733163246e-14L);
+  ValueSink(-27.0L, -6.720637435262040376484e-14L);
+  ValueSink(-26.0L, -1.894685885674978242162e-13L);
+  ValueSink(-25.0L, -5.348899755340216640865e-13L);
+  ValueSink(-24.0L, -1.512305893999705766268e-12L);
+  ValueSink(-23.0L, -4.282684795665672618927e-12L);
+  ValueSink(-22.0L, -1.214937895620437265858e-11L);
+  ValueSink(-21.0L, -3.453201267146756266712e-11L);
+  ValueSink(-20.0L, -9.835525290649881691934e-11L);
+  ValueSink(-19.0L, -2.807829097060795268028e-10L);
+  ValueSink(-18.0L, -8.036090344828677657508e-10L);
+  ValueSink(-17.0L, -2.306431989821654494494e-09L);
+  ValueSink(-16.0L, -6.640487249441042785668e-09L);
+  ValueSink(-15.0L, -1.918627892147866976706e-08L);
+  ValueSink(-14.0L, -5.565631111145182119369e-08L);
+  ValueSink(-13.0L, -1.621866218801432870747e-07L);
+  ValueSink(-12.0L, -4.751081824672493933550e-07L);
+  ValueSink(-11.0L, -1.400300304247441775031e-06L);
+  ValueSink(-10.0L, -4.156968929685324276268e-06L);
+  ValueSink(-9.0L, -1.244735417800627212619e-05L);
+  ValueSink(-8.0L, -3.766562284392490177815e-05L);
+  ValueSink(-7.6L, -5.885877207538391845951e-05L);
+  ValueSink(-7.2L, -9.218811688716204240148e-05L);
+  ValueSink(-7.0L, -1.154817316103382164251e-04L);
+  ValueSink(-6.8L, -1.447577922448995783522e-04L);
+  ValueSink(-6.4L, -2.279479559293728466322e-04L);
+  ValueSink(-6.0L, -3.600824521626586591688e-04L);
+  ValueSink(-5.6L, -5.708401696482116741887e-04L);
+  ValueSink(-5.2L, -9.086216124486595849240e-04L);
+  ValueSink(-5.0L, -1.148295591275325797160e-03L);
+  ValueSink(-4.8L, -1.452993936878324219239e-03L);
+  ValueSink(-4.4L, -2.336010040655828995865e-03L);
+  ValueSink(-4.0L, -3.779352409848906479661e-03L);
+  ValueSink(-3.0L, -1.304838109419703740879e-02L);
+  ValueSink(-2.0L, -4.890051070806111957565e-02L);
+  ValueSink(-1.0L, -2.193839343955202732906e-01L);
+  ValueSink(-0.8L, -3.105965785455430346377e-01L);
+  ValueSink(-0.6L, -4.543795031894021081123e-01L);
+  ValueSink(-0.4L, -7.023801188656624785609e-01L);
+  ValueSink(-0.2L, -1.222650544183893088403e+00L);
+  ValueSink(0.2L, -8.217605879024003156337e-01L);
+  ValueSink(0.4L, 1.047652186193247931785e-01L);
+  ValueSink(0.6L, 7.698812899373594372085e-01L);
+  ValueSink(0.8L, 1.347396548212325938189e+00L);
+  ValueSink(1.0L, 1.895117816355936755501e+00L);
+  ValueSink(2.0L, 4.954234356001890163083e+00L);
+  ValueSink(3.0L, 9.933832570625416555181e+00L);
+  ValueSink(4.0L, 1.963087447005622002659e+01L);
+  ValueSink(4.4L, 2.600897327160514597814e+01L);
+  ValueSink(4.8L, 3.469788987377532763282e+01L);
+  ValueSink(5.0L, 4.018527535580317745084e+01L);
+  ValueSink(5.2L, 4.662485050579674775162e+01L);
+  ValueSink(5.6L, 6.310178597429926148840e+01L);
+  ValueSink(6.0L, 8.598976214243920479507e+01L);
+  ValueSink(6.4L, 1.179348657001818873463e+02L);
+  ValueSink(6.8L, 1.627070875714314148464e+02L);
+  ValueSink(7.0L, 1.915047433355013959122e+02L);
+  ValueSink(7.2L, 2.256878077010638103811e+02L);
+  ValueSink(7.6L, 3.145718784980834849818e+02L);
+  ValueSink(8.0L, 4.403798995348382690485e+02L);
+  ValueSink(9.0L, 1.037878290717089587836e+03L);
+  ValueSink(10.0L, 2.492228976241877758557e+03L);
+  ValueSink(11.0L, 6.071406374098611509194e+03L);
+  ValueSink(12.0L, 1.495953266639752885414e+04L);
+  ValueSink(13.0L, 3.719768849068903560351e+04L);
+  ValueSink(14.0L, 9.319251363396537127670e+04L);
+  ValueSink(15.0L, 2.349558524907683035678e+05L);
+  ValueSink(16.0L, 5.955609986708370020096e+05L);
+  ValueSink(17.0L, 1.516637894042516884497e+06L);
+  ValueSink(18.0L, 3.877904330597443502484e+06L);
+  ValueSink(19.0L, 9.950907251046844759912e+06L);
+  ValueSink(20.0L, 2.561565266405658881740e+07L);
+  ValueSink(21.0L, 6.612718635548492133603e+07L);
+  ValueSink(22.0L, 1.711446713003636685462e+08L);
+  ValueSink(23.0L, 4.439663698302712207660e+08L);
+  ValueSink(24.0L, 1.154115391849182948354e+09L);
+  ValueSink(25.0L, 3.005950906525548688835e+09L);
+  ValueSink(26.0L, 7.842940991898186370730e+09L);
+  ValueSink(27.0L, 2.049649711988081235811e+10L);
+  ValueSink(28.0L, 5.364511859231469415501e+10L);
+  ValueSink(29.0L, 1.405991957584069047570e+11L);
+  ValueSink(30.0L, 3.689732094072741971314e+11L);
+  ValueSink(31.0L, 9.694555759683939658403e+11L);
+  ValueSink(32.0L, 2.550043566357786926746e+12L);
+  ValueSink(33.0L, 6.714640184076497560024e+12L);
+  ValueSink(34.0L, 1.769803724411626854134e+13L);
+  ValueSink(35.0L, 4.669055014466159545135e+13L);
+  ValueSink(36.0L, 1.232852079912097685471e+14L);
+  ValueSink(37.0L, 3.257988998672263995972e+14L);
+  ValueSink(38.0L, 8.616388199965786544800e+14L);
+  ValueSink(39.0L, 2.280446200301902594727e+15L);
+  ValueSink(40.0L, 6.039718263611241579102e+15L);
+  ValueSink(41.0L, 1.600664914324504110840e+16L);
+  ValueSink(42.0L, 4.244796092136850758594e+16L);
+  ValueSink(43.0L, 1.126348290166966760156e+17L);
+  ValueSink(44.0L, 2.990444718632336677188e+17L);
+  ValueSink(45.0L, 7.943916035704453773125e+17L);
+  ValueSink(46.0L, 2.111342388647824196000e+18L);
+  ValueSink(47.0L, 5.614329680810343109000e+18L);
+  ValueSink(48.0L, 1.493630213112993141900e+19L);
+  ValueSink(49.0L, 3.975442747903744833200e+19L);
+  ValueSink(50.0L, 1.058563689713169095760e+20L);
+}
+
+/// \brief Calls the provided functor for pairs of x and std::expint(x)
+/// generated with msvc 191627027 (double version).
+///
+/// \tparam BinaryFunction binary function with two long double arguments.
+/// \param ValueSink function that will be invoked for all pairs of
+/// x, std::expint(x).
+template <class BinaryFunction>
+void compareWithMSVC(BinaryFunction&& ValueSink) {
+  //Generated with msvc 191627027
+  ValueSink(-50.0L, -3.78326402955045910e-24L);
+  ValueSink(-49.0L, -1.04898116423680238e-23L);
+  ValueSink(-48.0L, -2.90966419040584234e-23L);
+  ValueSink(-47.0L, -8.07419784272581273e-23L);
+  ValueSink(-46.0L, -2.24153175974429984e-22L);
+  ValueSink(-45.0L, -6.22569080946238381e-22L);
+  ValueSink(-44.0L, -1.72995987428164762e-21L);
+  ValueSink(-43.0L, -4.80949655695001805e-21L);
+  ValueSink(-42.0L, -1.33779088100117746e-20L);
+  ValueSink(-41.0L, -3.72316677645997796e-20L);
+  ValueSink(-40.0L, -1.03677326145165702e-19L);
+  ValueSink(-39.0L, -2.88877930152270071e-19L);
+  ValueSink(-38.0L, -8.05410691429074992e-19L);
+  ValueSink(-37.0L, -2.24702069758857136e-18L);
+  ValueSink(-36.0L, -6.27333900976224210e-18L);
+  ValueSink(-35.0L, -1.75270593899473712e-17L);
+  ValueSink(-34.0L, -4.90067611839278744e-17L);
+  ValueSink(-33.0L, -1.37138434844874676e-16L);
+  ValueSink(-32.0L, -3.84096180122506713e-16L);
+  ValueSink(-31.0L, -1.07676703861623826e-15L);
+  ValueSink(-30.0L, -3.02155201068881243e-15L);
+  ValueSink(-29.0L, -8.48775977835356183e-15L);
+  ValueSink(-28.0L, -2.38694151193373298e-14L);
+  ValueSink(-27.0L, -6.72063743526203904e-14L);
+  ValueSink(-26.0L, -1.89468588567497827e-13L);
+  ValueSink(-25.0L, -5.34889975534021666e-13L);
+  ValueSink(-24.0L, -1.51230589399970585e-12L);
+  ValueSink(-23.0L, -4.28268479566567217e-12L);
+  ValueSink(-22.0L, -1.21493789562043715e-11L);
+  ValueSink(-21.0L, -3.45320126714675594e-11L);
+  ValueSink(-20.0L, -9.83552529064988154e-11L);
+  ValueSink(-19.0L, -2.80782909706079540e-10L);
+  ValueSink(-18.0L, -8.03609034482867693e-10L);
+  ValueSink(-17.0L, -2.30643198982165473e-09L);
+  ValueSink(-16.0L, -6.64048724944104267e-09L);
+  ValueSink(-15.0L, -1.91862789214786695e-08L);
+  ValueSink(-14.0L, -5.56563111114518230e-08L);
+  ValueSink(-13.0L, -1.62186621880143277e-07L);
+  ValueSink(-12.0L, -4.75108182467249308e-07L);
+  ValueSink(-11.0L, -1.40030030424744183e-06L);
+  ValueSink(-10.0L, -4.15696892968532464e-06L);
+  ValueSink(-9.0L, -1.24473541780062723e-05L);
+  ValueSink(-8.0L, -3.76656228439249064e-05L);
+  ValueSink(-7.6L, -5.88587720753839093e-05L);
+  ValueSink(-7.2L, -9.21881168871620238e-05L);
+  ValueSink(-7.0L, -1.15481731610338203e-04L);
+  ValueSink(-6.8L, -1.44757792244899590e-04L);
+  ValueSink(-6.4L, -2.27947955929372760e-04L);
+  ValueSink(-6.0L, -3.60082452162658673e-04L);
+  ValueSink(-5.6L, -5.70840169648211796e-04L);
+  ValueSink(-5.2L, -9.08621612448659468e-04L);
+  ValueSink(-5.0L, -1.14829559127532571e-03L);
+  ValueSink(-4.8L, -1.45299393687832457e-03L);
+  ValueSink(-4.4L, -2.33601004065582757e-03L);
+  ValueSink(-4.0L, -3.77935240984890626e-03L);
+  ValueSink(-3.0L, -1.30483810941970351e-02L);
+  ValueSink(-2.0L, -4.89005107080611248e-02L);
+  ValueSink(-1.0L, -2.19383934395520258e-01L);
+  ValueSink(-0.8L, -3.10596578545543012e-01L);
+  ValueSink(-0.6L, -4.54379503189402012e-01L);
+  ValueSink(-0.4L, -7.02380118865662428e-01L);
+  ValueSink(-0.2L, -1.22265054418389285e+00L);
+  ValueSink(0.2L, -8.21760587902400141e-01L);
+  ValueSink(0.4L, 1.04765218619324890e-01L);
+  ValueSink(0.6L, 7.69881289937359381e-01L);
+  ValueSink(0.8L, 1.34739654821232602e+00L);
+  ValueSink(1.0L, 1.89511781635593657e+00L);
+  ValueSink(2.0L, 4.95423435600188977e+00L);
+  ValueSink(3.0L, 9.93383257062541603e+00L);
+  ValueSink(4.0L, 1.96308744700562166e+01L);
+  ValueSink(4.4L, 2.60089732716051500e+01L);
+  ValueSink(4.8L, 3.46978898737753170e+01L);
+  ValueSink(5.0L, 4.01852753558031850e+01L);
+  ValueSink(5.2L, 4.66248505057967506e+01L);
+  ValueSink(5.6L, 6.31017859742992542e+01L);
+  ValueSink(6.0L, 8.59897621424392042e+01L);
+  ValueSink(6.4L, 1.17934865700181916e+02L);
+  ValueSink(6.8L, 1.62707087571431401e+02L);
+  ValueSink(7.0L, 1.91504743335501388e+02L);
+  ValueSink(7.2L, 2.25687807701063832e+02L);
+  ValueSink(7.6L, 3.14571878498083493e+02L);
+  ValueSink(8.0L, 4.40379899534838273e+02L);
+  ValueSink(9.0L, 1.03787829071708961e+03L);
+  ValueSink(10.0L, 2.49222897624187817e+03L);
+  ValueSink(11.0L, 6.07140637409861210e+03L);
+  ValueSink(12.0L, 1.49595326663975284e+04L);
+  ValueSink(13.0L, 3.71976884906890409e+04L);
+  ValueSink(14.0L, 9.31925136339653691e+04L);
+  ValueSink(15.0L, 2.34955852490768331e+05L);
+  ValueSink(16.0L, 5.95560998670836911e+05L);
+  ValueSink(17.0L, 1.51663789404251729e+06L);
+  ValueSink(18.0L, 3.87790433059744304e+06L);
+  ValueSink(19.0L, 9.95090725104684383e+06L);
+  ValueSink(20.0L, 2.56156526640565880e+07L);
+  ValueSink(21.0L, 6.61271863554849029e+07L);
+  ValueSink(22.0L, 1.71144671300363690e+08L);
+  ValueSink(23.0L, 4.43966369830271244e+08L);
+  ValueSink(24.0L, 1.15411539184918308e+09L);
+  ValueSink(25.0L, 3.00595090652554893e+09L);
+  ValueSink(26.0L, 7.84294099189818668e+09L);
+  ValueSink(27.0L, 2.04964971198808098e+10L);
+  ValueSink(28.0L, 5.36451185923146973e+10L);
+  ValueSink(29.0L, 1.40599195758406891e+11L);
+  ValueSink(30.0L, 3.68973209407274170e+11L);
+  ValueSink(31.0L, 9.69455575968393921e+11L);
+  ValueSink(32.0L, 2.55004356635778662e+12L);
+  ValueSink(33.0L, 6.71464018407649609e+12L);
+  ValueSink(34.0L, 1.76980372441162656e+13L);
+  ValueSink(35.0L, 4.66905501446615938e+13L);
+  ValueSink(36.0L, 1.23285207991209750e+14L);
+  ValueSink(37.0L, 3.25798899867226438e+14L);
+  ValueSink(38.0L, 8.61638819996578750e+14L);
+  ValueSink(39.0L, 2.28044620030190250e+15L);
+  ValueSink(40.0L, 6.03971826361124100e+15L);
+  ValueSink(41.0L, 1.60066491432450440e+16L);
+  ValueSink(42.0L, 4.24479609213685120e+16L);
+  ValueSink(43.0L, 1.12634829016696704e+17L);
+  ValueSink(44.0L, 2.99044471863233664e+17L);
+  ValueSink(45.0L, 7.94391603570445312e+17L);
+  ValueSink(46.0L, 2.11134238864782413e+18L);
+  ValueSink(47.0L, 5.61432968081034240e+18L);
+  ValueSink(48.0L, 1.49363021311299318e+19L);
+  ValueSink(49.0L, 3.97544274790374523e+19L);
+  ValueSink(50.0L, 1.05856368971316904e+20L);
+}
+
+/// \brief Calls the provided functor for the argument-value pairs from
+/// Harris, Frank E. "Tables of the exponential integral Ei (x)."
+/// Mathematical Tables and Other Aids to Computation 11.57 (1957): 9-16.
+///
+/// The values for +7.6 and -7.6 are omitted because they do not coincide
+/// with the values generated by gcc.
+/// \tparam BinaryFunction binary function with two long double arguments.
+/// \param ValueSink function that will be invoked for all pairs of
+/// x, expint(x).
+template <class BinaryFunction>
+void compareWithHarris1957(BinaryFunction&& ValueSink) {
+  ValueSink(-50.0L, -3.78326402955045902e-24L);
+  ValueSink(-49.0L, -1.04898116423680237e-23L);
+  ValueSink(-48.0L, -2.90966419040584234e-23L);
+  ValueSink(-47.0L, -8.07419784272581395e-23L);
+  ValueSink(-46.0L, -2.24153175974429975e-22L);
+  ValueSink(-45.0L, -6.22569080946238364e-22L);
+  ValueSink(-44.0L, -1.72995987428164776e-21L);
+  ValueSink(-43.0L, -4.80949655695001785e-21L);
+  ValueSink(-42.0L, -1.33779088100117751e-20L);
+  ValueSink(-41.0L, -3.72316677645997772e-20L);
+  ValueSink(-40.0L, -1.03677326145165697e-19L);
+  ValueSink(-39.0L, -2.88877930152270100e-19L);
+  ValueSink(-38.0L, -8.05410691429074987e-19L);
+  ValueSink(-37.0L, -2.24702069758857122e-18L);
+  ValueSink(-36.0L, -6.27333900976224159e-18L);
+  ValueSink(-35.0L, -1.75270593899473720e-17L);
+  ValueSink(-34.0L, -4.90067611839278771e-17L);
+  ValueSink(-33.0L, -1.37138434844874657e-16L);
+  ValueSink(-32.0L, -3.84096180122506683e-16L);
+  ValueSink(-31.0L, -1.07676703861623826e-15L);
+  ValueSink(-30.0L, -3.02155201068881254e-15L);
+  ValueSink(-29.0L, -8.48775977835356284e-15L);
+  ValueSink(-28.0L, -2.38694151193373316e-14L);
+  ValueSink(-27.0L, -6.72063743526204038e-14L);
+  ValueSink(-26.0L, -1.89468588567497824e-13L);
+  ValueSink(-25.0L, -5.34889975534021664e-13L);
+  ValueSink(-24.0L, -1.51230589399970577e-12L);
+  ValueSink(-23.0L, -4.28268479566567262e-12L);
+  ValueSink(-22.0L, -1.21493789562043727e-11L);
+  ValueSink(-21.0L, -3.45320126714675627e-11L);
+  ValueSink(-20.0L, -9.83552529064988169e-11L);
+  ValueSink(-19.0L, -2.80782909706079527e-10L);
+  ValueSink(-18.0L, -8.03609034482867766e-10L);
+  ValueSink(-17.0L, -2.30643198982165449e-09L);
+  ValueSink(-16.0L, -6.64048724944104278e-09L);
+  ValueSink(-15.0L, -1.91862789214786698e-08L);
+  ValueSink(-14.0L, -5.56563111114518212e-08L);
+  ValueSink(-13.0L, -1.62186621880143287e-07L);
+  ValueSink(-12.0L, -4.75108182467249393e-07L);
+  ValueSink(-11.0L, -1.40030030424744178e-06L);
+  ValueSink(-10.0L, -4.15696892968532428e-06L);
+  ValueSink(-9.0L, -1.24473541780062721e-05L);
+  ValueSink(-8.0L, -3.76656228439249018e-05L);
+  //ValueSink( -7.6L, -5.88587720753838951e-05L);
+  ValueSink(-7.2L, -9.21881168871620423e-05L);
+  ValueSink(-6.8L, -1.44757792244899578e-04L);
+  ValueSink(-6.4L, -2.27947955929372847e-04L);
+  ValueSink(-6.0L, -3.60082452162658659e-04L);
+  ValueSink(-5.6L, -5.70840169648211674e-04L);
+  ValueSink(-5.2L, -9.08621612448659585e-04L);
+  ValueSink(-4.8L, -1.45299393687832422e-03L);
+  ValueSink(-4.4L, -2.33601004065582899e-03L);
+  ValueSink(-4.0L, -3.77935240984890648e-03L);
+  ValueSink(-3.0L, -1.30483810941970374e-02L);
+  ValueSink(-2.0L, -4.89005107080611196e-02L);
+  ValueSink(-1.0L, -2.19383934395520274e-01L);
+  ValueSink(1.0L, 1.89511781635593676e+00L);
+  ValueSink(2.0L, 4.95423435600189016e+00L);
+  ValueSink(3.0L, 9.93383257062541656e+00L);
+  ValueSink(4.0L, 1.96308744700562200e+01L);
+  ValueSink(4.4L, 2.60089732716051460e+01L);
+  ValueSink(4.8L, 3.46978898737753276e+01L);
+  ValueSink(5.2L, 4.66248505057967478e+01L);
+  ValueSink(5.6L, 6.31017859742992615e+01L);
+  ValueSink(6.0L, 8.59897621424392048e+01L);
+  ValueSink(6.4L, 1.17934865700181887e+02L);
+  ValueSink(6.8L, 1.62707087571431415e+02L);
+  ValueSink(7.2L, 2.25687807701063810e+02L);
+  //ValueSink(  7.6L,  3.14571878498083578e+02L);
+  ValueSink(8.0L, 4.40379899534838269e+02L);
+  ValueSink(9.0L, 1.03787829071708959e+03L);
+  ValueSink(10.0L, 2.49222897624187776e+03L);
+  ValueSink(11.0L, 6.07140637409861151e+03L);
+  ValueSink(12.0L, 1.49595326663975289e+04L);
+  ValueSink(13.0L, 3.71976884906890356e+04L);
+  ValueSink(14.0L, 9.31925136339653713e+04L);
+  ValueSink(15.0L, 2.34955852490768304e+05L);
+  ValueSink(16.0L, 5.95560998670837002e+05L);
+  ValueSink(17.0L, 1.51663789404251688e+06L);
+  ValueSink(18.0L, 3.87790433059744350e+06L);
+  ValueSink(19.0L, 9.95090725104684476e+06L);
+  ValueSink(20.0L, 2.56156526640565888e+07L);
+  ValueSink(21.0L, 6.61271863554849213e+07L);
+  ValueSink(22.0L, 1.71144671300363668e+08L);
+  ValueSink(23.0L, 4.43966369830271221e+08L);
+  ValueSink(24.0L, 1.15411539184918295e+09L);
+  ValueSink(25.0L, 3.00595090652554869e+09L);
+  ValueSink(26.0L, 7.84294099189818637e+09L);
+  ValueSink(27.0L, 2.04964971198808124e+10L);
+  ValueSink(28.0L, 5.36451185923146942e+10L);
+  ValueSink(29.0L, 1.40599195758406905e+11L);
+  ValueSink(30.0L, 3.68973209407274197e+11L);
+  ValueSink(31.0L, 9.69455575968393966e+11L);
+  ValueSink(32.0L, 2.55004356635778693e+12L);
+  ValueSink(33.0L, 6.71464018407649756e+12L);
+  ValueSink(34.0L, 1.76980372441162685e+13L);
+  ValueSink(35.0L, 4.66905501446615954e+13L);
+  ValueSink(36.0L, 1.23285207991209769e+14L);
+  ValueSink(37.0L, 3.25798899867226400e+14L);
+  ValueSink(38.0L, 8.61638819996578654e+14L);
+  ValueSink(39.0L, 2.28044620030190260e+15L);
+  ValueSink(40.0L, 6.03971826361124158e+15L);
+  ValueSink(41.0L, 1.60066491432450411e+16L);
+  ValueSink(42.0L, 4.24479609213685076e+16L);
+  ValueSink(43.0L, 1.12634829016696676e+17L);
+  ValueSink(44.0L, 2.99044471863233668e+17L);
+  ValueSink(45.0L, 7.94391603570445377e+17L);
+  ValueSink(46.0L, 2.11134238864782419e+18L);
+  ValueSink(47.0L, 5.61432968081034311e+18L);
+  ValueSink(48.0L, 1.49363021311299314e+19L);
+  ValueSink(49.0L, 3.97544274790374484e+19L);
+  ValueSink(50.0L, 1.05856368971316910e+20L);
+}
+
+/// \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::experimental::expint(Root - Tolerance) < T(0));
+  assert(std::experimental::expint(Root + Tolerance) > T(0));
+}
+
+template <class T>
+void testReferenceValues(T RelTolerance) {
+  const auto CompareValues = [RelTolerance](long double x, long double y) {
+    const T ExpintAtX = std::experimental::expint(static_cast<T>(x));
+    const T Tolerance = RelTolerance * std::abs(static_cast<T>(y));
+    assert(std::abs(ExpintAtX - static_cast<T>(y)) < Tolerance);
+  };
+
+  compareWithBoost(CompareValues);
+  compareWithGCC(CompareValues);
+  compareWithMSVC(CompareValues);
+  compareWithHarris1957(CompareValues);
+
+  compareRootWithCody1969(RelTolerance);
+}
+
+int main(int, char**) {
+  testReferenceValues<float>(1e-6f);
+  testReferenceValues<double>(1e-14);
+}

libcxx/test/std/experimental/numerics/c.math/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
+//
+//===----------------------------------------------------------------------===//
+
+// <cmath>
+
+// double      expint(double arg);
+// float       expint(float arg);
+// long double expint(long double arg);
+// float       expintf(float arg);
+// long double expintl(long double arg);
+
+#include <experimental/cmath>
+#include <cassert>
+
+template <class T>
+void testLimits() {
+  const T NegativeLimit =
+      std::experimental::expint(-std::numeric_limits<T>::infinity());
+  assert(NegativeLimit <= 0);
+  assert(NegativeLimit >= -std::numeric_limits<T>::min());
+
+  const T PositiveLimit =
+      std::experimental::expint(std::numeric_limits<T>::infinity());
+  assert(PositiveLimit > T(0));
+  assert(std::isinf(PositiveLimit));
+
+  const T ExpintAtZero = std::experimental::expint(T(0));
+  assert(ExpintAtZero < T(0));
+  assert(std::isinf(ExpintAtZero));
+}
+
+int main(int, char**) {
+  testLimits<float>();
+  testLimits<double>();
+  testLimits<long double>();
+}