diff --git a/libcxx/include/experimental/cmath b/libcxx/include/experimental/cmath new file mode 100644 --- /dev/null +++ b/libcxx/include/experimental/cmath @@ -0,0 +1,66 @@ +// -*- 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 + +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 +#include + +_LIBCPP_BEGIN_NAMESPACE_EXPERIMENTAL + +#if _LIBCPP_STD_VER > 14 + +template +_Real __libcpp_expint(_Real __x); + +inline _LIBCPP_INLINE_VISIBILITY float expintf(float __lcpp_x) { + return _VSTD_EXPERIMENTAL::__libcpp_expint(__lcpp_x); +} + +inline _LIBCPP_INLINE_VISIBILITY float expint(float __lcpp_x) { + return _VSTD_EXPERIMENTAL::__libcpp_expint(__lcpp_x); +} + +inline _LIBCPP_INLINE_VISIBILITY double expint(double __lcpp_x) { + return _VSTD_EXPERIMENTAL::__libcpp_expint(__lcpp_x); +} + +inline _LIBCPP_INLINE_VISIBILITY long double expintl(long double __lcpp_x) { + return _VSTD_EXPERIMENTAL::__libcpp_expint(__lcpp_x); +} + +inline _LIBCPP_INLINE_VISIBILITY long double expint(long double __lcpp_x) { + return _VSTD_EXPERIMENTAL::__libcpp_expint(__lcpp_x); +} + +#endif + +_LIBCPP_END_NAMESPACE_EXPERIMENTAL + +#endif // _LIBCPP_EXPERIMENTAL_CMATH diff --git a/libcxx/src/experimental/expint.cpp b/libcxx/src/experimental/expint.cpp new file mode 100644 --- /dev/null +++ b/libcxx/src/experimental/expint.cpp @@ -0,0 +1,212 @@ +//===------------------------ 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, and its explicit instantiations. +/// +//===----------------------------------------------------------------------===// + +#include +#include +#include + +_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 +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 +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 +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::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 +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 +_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); +template double __libcpp_expint(double); +template long double __libcpp_expint(long double); + +_LIBCPP_END_NAMESPACE_EXPERIMENTAL diff --git a/libcxx/test/std/experimental/numerics/c.math/expint_bounds.pass.cpp b/libcxx/test/std/experimental/numerics/c.math/expint_bounds.pass.cpp new file mode 100644 --- /dev/null +++ b/libcxx/test/std/experimental/numerics/c.math/expint_bounds.pass.cpp @@ -0,0 +1,213 @@ +//===----------------------------------------------------------------------===// +// +// 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 +// +//===----------------------------------------------------------------------===// + +// + +// 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 +#include +#include + +template +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 +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 +void testSingularity(T RelTolerance, T AbsTolerance) { + for (T x : {T(1e-3), T(1e-2), T(0.1), T(1), T(10)}) + testSingularity(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 +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 +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::min()); + assert(std::experimental::expint(LogMin) < std::numeric_limits::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 +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 +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 +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 +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::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::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 +void testAllBounds(T RelTolerance, T AbsTolerance) { + testSingularity(RelTolerance, AbsTolerance); + testBoundsNegativeX(RelTolerance, AbsTolerance); + testBoundsPositiveSmallX(RelTolerance, AbsTolerance); + testBoundsPositiveLargeX(RelTolerance, AbsTolerance); +} + +int main(int, char**) { + testAllBounds(1e-6f, 0.0f); + testAllBounds(1e-14, 0.0); +} diff --git a/libcxx/test/std/experimental/numerics/c.math/expint_comparisons.pass.cpp b/libcxx/test/std/experimental/numerics/c.math/expint_comparisons.pass.cpp new file mode 100644 --- /dev/null +++ b/libcxx/test/std/experimental/numerics/c.math/expint_comparisons.pass.cpp @@ -0,0 +1,585 @@ +//===----------------------------------------------------------------------===// +// +// 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 +// +//===----------------------------------------------------------------------===// + +// + +// 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 +#include + +/// \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 +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 +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 +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 +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 +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 +void testReferenceValues(T RelTolerance) { + const auto CompareValues = [RelTolerance](long double x, long double y) { + const T ExpintAtX = std::experimental::expint(static_cast(x)); + const T Tolerance = RelTolerance * std::abs(static_cast(y)); + assert(std::abs(ExpintAtX - static_cast(y)) < Tolerance); + }; + + compareWithBoost(CompareValues); + compareWithGCC(CompareValues); + compareWithMSVC(CompareValues); + compareWithHarris1957(CompareValues); + + compareRootWithCody1969(RelTolerance); +} + +int main(int, char**) { + testReferenceValues(1e-6f); + testReferenceValues(1e-14); +} diff --git a/libcxx/test/std/experimental/numerics/c.math/expint_limits.pass.cpp b/libcxx/test/std/experimental/numerics/c.math/expint_limits.pass.cpp new file mode 100644 --- /dev/null +++ b/libcxx/test/std/experimental/numerics/c.math/expint_limits.pass.cpp @@ -0,0 +1,41 @@ +//===----------------------------------------------------------------------===// +// +// 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 +// +//===----------------------------------------------------------------------===// + +// + +// 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 +#include + +template +void testLimits() { + const T NegativeLimit = + std::experimental::expint(-std::numeric_limits::infinity()); + assert(NegativeLimit <= 0); + assert(NegativeLimit >= -std::numeric_limits::min()); + + const T PositiveLimit = + std::experimental::expint(std::numeric_limits::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(); + testLimits(); + testLimits(); +}