Index: libcxx/include/cmath =================================================================== --- libcxx/include/cmath +++ libcxx/include/cmath @@ -296,6 +296,12 @@ float truncf(float x); long double truncl(long double x); +// C++17 Mathematical special functions [sf.cmath] + +floating_point expint (arithmetic x); +float expintf(float x); +long double expintl(long double x); + } // std */ @@ -303,6 +309,8 @@ #include <__config> #include #include +#include +#include #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER) #pragma GCC system_header @@ -606,6 +614,221 @@ return isfinite(__lcpp_x); } +#if _LIBCPP_STD_VER > 14 + +// expint + +/// \brief Returns the Euler-Mascheroni constant. +/// +/// The value is taken from: +/// Sweeney, Dura W. "On the computation of Euler’s constant." +/// Mathematics of Computation 17.82 (1963): 170-178. +/// \return the Euler-Mascheroni constant \f$ \gamma = 0.5772156649\ldots \f$ +template +inline _LIBCPP_INLINE_VISIBILITY _Real __libcpp_euler_gamma() _NOEXCEPT { + return _Real(0.57721566490153286060651209008240243104215933593992l); +} + +/// \brief Evaluates std::expint(x) for x < -1. +/// +/// This function evaluates the following expression: +/// \f[ +/// \mathrm{Ei} (x) = -e^{x} \left( +/// \frac{1}{x+1-}\frac{1}{x+3-}\frac{4}{x+5-} \cdots\right) +/// \f] +/// See Press, William H., et al. Numerical recipes 3rd edition: The art of +/// scientific computing. Cambridge university press, 2007, p. 267., eq. (6.3.5) +/// The implementation uses two iterations of the form +/// \f[ +/// s^{(i)} = 1 - x + 2i - \frac{i^2}{s^{(i-1)}} +/// \f] +/// with different start values \f$ s_0^{(1)} = 1 -x + 2 \f$ and +/// \f$ s_1^{(1)} = 1 -x + 2 - \frac{1}{1-x} \f$. +/// Both of them are bound by +/// \f$ 1 - x +i < s^{(i)} \leq 1 - x + 2i \f$. +/// Therefore, the difference monotonically decreasing +/// \f[ +/// s_0^{(i)} - s_1^{(i)} = \frac{i^2}{s_0^{(i-1)} s_1^{(i-1)}} +/// \left( s_0^{(i)} - s_1^{(i)} \right). +/// \f] +/// \param [in] __x Argument \f$ x < -1\f$ of \f$\mathrm{Ei} (x) \f$. +/// \return continued fraction representation of \f$ \mathrm{Ei} (x) \f$ +/// for \f$ x < -1 \f$ +template +inline _LIBCPP_INLINE_VISIBILITY _Real +__libcpp_expint_continued_fraction(_Real __x) _NOEXCEPT { + const int __max_iter = 1000; + const _Real __mxp1 = -__x + _Real(1); + if (_VSTD::isinf(__mxp1 + _Real(2 * __max_iter))) + return _Real(0); + + // See algorithm expint.h on p. 268 of + // Press, William H., et al. Numerical recipes 3rd edition: The art of + // scientific computing. Cambridge university press, 2007. + // a -> -__t0 + // b -> __t1 + // c -> _s0 + // d -> 1/__ s1 + + _Real __h = _Real(1) / __mxp1; + _Real __s0 = __mxp1 + _Real(2); + _Real __s1 = __s0 - __h; + __h *= __s0 / __s1; + + for (int __i = 2; __i < __max_iter; ++__i) { + const _Real __t0 = _Real(__i * __i); + const _Real __t1 = __mxp1 + _Real(2 * __i); + __s0 = __t1 - __t0 / __s0; + __s1 = __t1 - __t0 / __s1; + if (__s0 == __s1) + return -__h * exp(__x); + + __h *= __s0 / __s1; + } + + // (See the estimates in the comment above.) + // The closer __x is to 0, the slower the convergence. + // We use this function only for __x < -1, which makes + // __x == -1 the worst case. + + _LIBCPP_UNREACHABLE(); +} + +/// \brief Evaluates std::expint(x) for x != 0. +/// +/// This function evaluates the following expression: +/// \f[ +/// \mathrm{Ei} (x) = \gamma + \ln \left|x\right| + \sum \limits_{n=1}^{\infty} +/// \frac{x^n}{n\cdotn!} +/// \f] +/// See Press, William H., et al. Numerical recipes 3rd edition: The art of +/// scientific computing. Cambridge university press, 2007, p. 267., eq. (6.3.6) +/// or +/// Harris, Frank E. "Tables of the exponential integral Ei (x)." +/// Mathematical Tables and Other Aids to Computation 11.57 (1957): 9-16. +/// \attention This representation converges slowly for large +/// \f$ \left|x\right| \f$! +/// \param [in] __x Argument \f$ x \neq 0 \f$ of \f$\mathrm{Ei} (x) \f$. +/// \return Evaluation of a series representation of \f$ \mathrm{Ei} (x) \f$ +/// for \f$ -1 < x < 1 \f$. +template +inline _LIBCPP_INLINE_VISIBILITY +_Real __libcpp_expint_power_series(_Real __x) _NOEXCEPT { + _Real __sum = _Real(0); + _Real __k = _Real(1); + _Real factor = _Real(1); + for (;;) { + factor *= __x / __k; + const _Real __prev = __sum; + __sum += factor / __k; + if (__prev == __sum) + return __sum + __libcpp_euler_gamma<_Real>() + + _VSTD::log(_VSTD::abs(__x)); + + __k += _Real(1); + } +} + +/// \brief Evaluates std::expint(x) for x > -std::log( +/// std::numeric_limits::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) _NOEXCEPT { + _Real __sum = _Real(0); + _Real __t = _Real(1); + _Real __k = _Real(1); + do { + __t *= __k / __x; + __sum += __t; + __k += _Real(1); + } while (__k < __x && __t > _VSTD::numeric_limits<_Real>::epsilon()); + + const _Real __overflow_limit = + (_Real(1) - _VSTD::numeric_limits<_Real>::epsilon()) * + _VSTD::log(_VSTD::numeric_limits<_Real>::max()); + + if (__x >= __overflow_limit) { + if (_VSTD::isinf(__x)) { + return _VSTD::numeric_limits<_Real>::infinity(); + } + // The first term is actually exp(x)/x in disguise. + // It is less precise but isn't subject to premature overflow. + return _VSTD::exp(__x - _VSTD::log(__x)) * (_Real(1) + __sum); + } + + return _VSTD::exp(__x) / __x * (_Real(1) + __sum); +} + +/// \brief Internal implementation of std::expint. +/// +/// The implementation is based on chapter 6.3, p. 266 et seqq. of +/// Press, William H., et al. Numerical recipes 3rd edition: The art of +/// scientific computing. Cambridge university press, 2007 +/// \tparam _Real float, double, or long double +/// \param __x Argument \f$ x \gg 1 \f$ of \f$\mathrm{Ei} (x) \f$. +/// \return std::expint(x). +template +inline _LIBCPP_INLINE_VISIBILITY +_Real __libcpp_expint(_Real __x) _NOEXCEPT { + if (__x < _Real(-1)) + return __libcpp_expint_continued_fraction(__x); + + const _Real __asymptotic = + -_VSTD::log(_VSTD::numeric_limits<_Real>::epsilon()); + if (__x < __asymptotic) + return __libcpp_expint_power_series(__x); + + return __libcpp_expint_asymptotic(__x); +} + +template +inline _LIBCPP_INLINE_VISIBILITY +typename __lazy_enable_if::value, __promote<_A1>>::type +expint(_A1 __lcpp_x) { + typedef typename __promote<_A1>::type __result_type; + return __libcpp_expint<__result_type>((__result_type)__lcpp_x); +} + +inline _LIBCPP_INLINE_VISIBILITY +float expintf(float __lcpp_x) { + return __libcpp_expint(__lcpp_x); +} + +inline _LIBCPP_INLINE_VISIBILITY +long double expintl(long double __lcpp_x) { + return __libcpp_expint(__lcpp_x); +} + +#endif // _LIBCPP_STD_VER > 14 + + #if _LIBCPP_STD_VER > 17 template constexpr Index: libcxx/test/std/numerics/c.math/c.math.expint/expint_bounds.pass.cpp =================================================================== --- /dev/null +++ libcxx/test/std/numerics/c.math/c.math.expint/expint_bounds.pass.cpp @@ -0,0 +1,209 @@ +//===----------------------------------------------------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +// UNSUPPORTED: c++98, c++03, c++11, c++14 + +// + +// floating_point expint (arithmetic x); // C++17 +// float expintf(float x); // C++17 +// long double expintl(long double x); // C++17 + +#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::expint(x))); + assert(!std::isinf(std::expint(-x))); + const T Difference = std::expint(x) - std::expint(-x); + const T LowerBound = x*(T(2) + x * x * (T(1)/ T(9) + x * x / T(300))); + const T UpperBound = T(2) * std::sinh(x); + const T AbsToleranceForDifference = + T(2) * AbsTolerance + + (std::abs(std::expint(x)) + std::abs(std::expint(-x))) * RelTolerance; + assert(isInsideInterval(Difference, LowerBound, UpperBound, T(0), + AbsToleranceForDifference)); +} + +template +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::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::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::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::expint(T(1)) + std::log(x) + std::exp(x) / x - std::exp(T(1)); + const T UpperBound = std::expint(T(1)) + std::log(x) + + (std::exp(x) - std::exp(T(1))) - (x - T(1)); + assert(isInsideInterval(std::expint(x), LowerBound, UpperBound, RelTolerance, + AbsTolerance)); +} + +template +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: + const T LogMax = std::log(std::numeric_limits::max()); + assert(!std::isinf(std::expint(LogMax))); + const T GivesInf = T(4) * LogMax * LogMax; + assert(std::isinf(std::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); +} Index: libcxx/test/std/numerics/c.math/c.math.expint/expint_comparisons.pass.cpp =================================================================== --- /dev/null +++ libcxx/test/std/numerics/c.math/c.math.expint/expint_comparisons.pass.cpp @@ -0,0 +1,195 @@ +//===----------------------------------------------------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +// UNSUPPORTED: c++98, c++03, c++11, c++14 + +// + +// floating_point expint (arithmetic x); // C++17 +// float expintf(float x); // C++17 +// long double expintl(long double x); // C++17 + +#include +#include +#include + +/// \brief Returns a map of x and std::expint(x) +/// +/// This table is generated with gcc-8.3.1 and boost-1.66.0 +/// and confirmed with msvc-192027508 and boost-1.69.0 +/// using boost::multiprecision::cpp_bin_float_oct +/// rounded to std::numeric_limits::max_digits10. +/// \note The key_type is float to make sure that it can be +/// exactly represented by all floating-point types. +/// The mapped_type is long double to avoid precision loss +/// but no overflow occurs when casting to float. +std::map generateExpintTestValues() +{ + return {{-5.000000000e+01f, -3.783264029550459018698967854021285780e-24L}, + {-4.900000000e+01f, -1.048981164236802369849688680874776641e-23L}, + {-4.800000000e+01f, -2.909664190405842339165765290916506568e-23L}, + {-4.700000000e+01f, -8.074197842725813951937477940562408284e-23L}, + {-4.600000000e+01f, -2.241531759744299745983944746971850400e-22L}, + {-4.500000000e+01f, -6.225690809462383643095308453348153894e-22L}, + {-4.400000000e+01f, -1.729959874281647757227851510531351102e-21L}, + {-4.300000000e+01f, -4.809496556950017850712335884632312421e-21L}, + {-4.200000000e+01f, -1.337790881001177514314392868306569266e-20L}, + {-4.100000000e+01f, -3.723166776459977716235516016523522486e-20L}, + {-4.000000000e+01f, -1.036773261451656972150641889145259771e-19L}, + {-3.900000000e+01f, -2.888779301522700998725095238399534096e-19L}, + {-3.800000000e+01f, -8.054106914290749868944421587855629671e-19L}, + {-3.700000000e+01f, -2.247020697588571222083353307904064249e-18L}, + {-3.600000000e+01f, -6.273339009762241588152637033384756530e-18L}, + {-3.500000000e+01f, -1.752705938994737200054832669092913667e-17L}, + {-3.400000000e+01f, -4.900676118392787715370624194027081175e-17L}, + {-3.300000000e+01f, -1.371384348448746572239655805859484591e-16L}, + {-3.200000000e+01f, -3.840961801225066831461234630025750207e-16L}, + {-3.100000000e+01f, -1.076767038616238261726105157668358207e-15L}, + {-3.000000000e+01f, -3.021552010688812544815825045153697921e-15L}, + {-2.900000000e+01f, -8.487759778353562840972762979415674942e-15L}, + {-2.800000000e+01f, -2.386941511933733162909373835426501483e-14L}, + {-2.700000000e+01f, -6.720637435262040377094883818453725297e-14L}, + {-2.600000000e+01f, -1.894685885674978242549641305135035184e-13L}, + {-2.500000000e+01f, -5.348899755340216640325471118085924169e-13L}, + {-2.400000000e+01f, -1.512305893999705766474308468724763618e-12L}, + {-2.300000000e+01f, -4.282684795665672618906770326270850700e-12L}, + {-2.200000000e+01f, -1.214937895620437265840775447119305968e-11L}, + {-2.100000000e+01f, -3.453201267146756266667883784577524543e-11L}, + {-2.000000000e+01f, -9.835525290649881690396987108894776074e-11L}, + {-1.900000000e+01f, -2.807829097060795267324152564574528010e-10L}, + {-1.800000000e+01f, -8.036090344828677657207213012471781597e-10L}, + {-1.700000000e+01f, -2.306431989821654494245004407041613170e-09L}, + {-1.600000000e+01f, -6.640487249441042785696129359642739491e-09L}, + {-1.500000000e+01f, -1.918627892147866977104309696741338816e-08L}, + {-1.400000000e+01f, -5.565631111145182114989920850651295593e-08L}, + {-1.300000000e+01f, -1.621866218801432870669074367074383031e-07L}, + {-1.200000000e+01f, -4.751081824672493932594612696661441836e-07L}, + {-1.100000000e+01f, -1.400300304247441775446564237904285784e-06L}, + {-1.000000000e+01f, -4.156968929685324277402859810278180384e-06L}, + {-9.000000000e+00f, -1.244735417800627212114856523810080257e-05L}, + {-8.000000000e+00f, -3.766562284392490177255799595075272089e-05L}, + {-7.599999905e+00f, -5.885877835522211132219909198852529146e-05L}, + {-7.199999809e+00f, -9.218813666493392496819998690716120347e-05L}, + {-7.000000000e+00f, -1.154817316103382164310114560437894989e-04L}, + {-6.800000191e+00f, -1.447577610043514442473650415536007813e-04L}, + {-6.400000095e+00f, -2.279479311702414515081458172733752388e-04L}, + {-6.000000000e+00f, -3.600824521626586592953941157717971888e-04L}, + {-5.599999905e+00f, -5.708402326224589986646571434159349354e-04L}, + {-5.199999809e+00f, -9.086218147950596268246438778225680700e-04L}, + {-5.000000000e+00f, -1.148295591275325797330561969819722076e-03L}, + {-4.800000191e+00f, -1.452993609857595605868097990432304282e-03L}, + {-4.400000095e+00f, -2.336009774551668672538799012916084760e-03L}, + {-4.000000000e+00f, -3.779352409848906478874860132466414852e-03L}, + {-3.000000000e+00f, -1.304838109419703741250074582864502295e-02L}, + {-2.000000000e+00f, -4.890051070806111956723983522804952231e-02L}, + {-1.000000000e+00f, -2.193839343955202736771637754601216490e-01L}, + {-8.000000119e-01f, -3.105965718500198010734802054347857172e-01L}, + {-6.000000238e-01f, -4.543794813815877227508100153022577629e-01L}, + {-4.000000060e-01f, -7.023801088771155253272499368438163604e-01L}, + {-2.000000030e-01f, -1.222650531983854271436328625153714048e+00L}, + {2.000000030e-01f, -8.217605697020810421808504320544601028e-01L}, + {4.000000060e-01f, 1.047652408492449863690536252242036180e-01L}, + {6.000000238e-01f, 7.698813623418547522320476636656036129e-01L}, + {8.000000119e-01f, 1.347396581375470007673701978355742968e+00L}, + {1.000000000e+00f, 1.895117816355936755466520934331634269e+00L}, + {2.000000000e+00f, 4.954234356001890163379505130227035276e+00L}, + {3.000000000e+00f, 9.933832570625416558008336019216765263e+00L}, + {4.000000000e+00f, 1.963087447005622002264572027972383885e+01L}, + {4.400000095e+00f, 2.600897503700524499824399694364832107e+01L}, + {4.800000191e+00f, 3.469789470216587397588224274169374529e+01L}, + {5.000000000e+00f, 4.018527535580317745509142179379586710e+01L}, + {5.199999809e+00f, 4.662484385677105401583505624620463816e+01L}, + {5.599999905e+00f, 6.310178136896516972405815088748916221e+01L}, + {6.000000000e+00f, 8.598976214243920480358340030799069004e+01L}, + {6.400000095e+00f, 1.179348746683721712003557710511492152e+02L}, + {6.800000191e+00f, 1.627071127553717604838712010833682320e+02L}, + {7.000000000e+00f, 1.915047433355013959530631482724569474e+02L}, + {7.199999809e+00f, 2.256877722182689852472513019665927695e+02L}, + {7.599999905e+00f, 3.145718534240306126656270737853903947e+02L}, + {8.000000000e+00f, 4.403798995348382689974245966593933924e+02L}, + {9.000000000e+00f, 1.037878290717089587657573226793622215e+03L}, + {1.000000000e+01f, 2.492228976241877759138440143998524849e+03L}, + {1.100000000e+01f, 6.071406374098611507964887284858515516e+03L}, + {1.200000000e+01f, 1.495953266639752885229246187605753281e+04L}, + {1.300000000e+01f, 3.719768849068903560439164528876347903e+04L}, + {1.400000000e+01f, 9.319251363396537129882452836392441672e+04L}, + {1.500000000e+01f, 2.349558524907683035782457458951611726e+05L}, + {1.600000000e+01f, 5.955609986708370018501610068484626103e+05L}, + {1.700000000e+01f, 1.516637894042516884432797432876246260e+06L}, + {1.800000000e+01f, 3.877904330597443502996466079950798035e+06L}, + {1.900000000e+01f, 9.950907251046844760026002538253063332e+06L}, + {2.000000000e+01f, 2.561565266405658882048112080409807183e+07L}, + {2.100000000e+01f, 6.612718635548492136250291987962721053e+07L}, + {2.200000000e+01f, 1.711446713003636684975370635366578572e+08L}, + {2.300000000e+01f, 4.439663698302712208698485234569711018e+08L}, + {2.400000000e+01f, 1.154115391849182948286759099974094906e+09L}, + {2.500000000e+01f, 3.005950906525548689841377604167345810e+09L}, + {2.600000000e+01f, 7.842940991898186370453025612588634511e+09L}, + {2.700000000e+01f, 2.049649711988081236484165294369597173e+10L}, + {2.800000000e+01f, 5.364511859231469415605083102066723533e+10L}, + {2.900000000e+01f, 1.405991957584069047339548366002198035e+11L}, + {3.000000000e+01f, 3.689732094072741970640063289108457470e+11L}, + {3.100000000e+01f, 9.694555759683939661661691071346617779e+11L}, + {3.200000000e+01f, 2.550043566357786926146742307679350278e+12L}, + {3.300000000e+01f, 6.714640184076497558707440521912031531e+12L}, + {3.400000000e+01f, 1.769803724411626854310341944075653955e+13L}, + {3.500000000e+01f, 4.669055014466159544500146290990063797e+13L}, + {3.600000000e+01f, 1.232852079912097685430891907640127096e+14L}, + {3.700000000e+01f, 3.257988998672263996790001683698134913e+14L}, + {3.800000000e+01f, 8.616388199965786544948016613356301502e+14L}, + {3.900000000e+01f, 2.280446200301902595340816714438768255e+15L}, + {4.000000000e+01f, 6.039718263611241578359231418510691294e+15L}, + {4.100000000e+01f, 1.600664914324504111069970545013283212e+16L}, + {4.200000000e+01f, 4.244796092136850759367705586573591072e+16L}, + {4.300000000e+01f, 1.126348290166966760275342380239567459e+17L}, + {4.400000000e+01f, 2.990444718632336675058132672827320204e+17L}, + {4.500000000e+01f, 7.943916035704453771510168303218332002e+17L}, + {4.600000000e+01f, 2.111342388647824195000286962846296846e+18L}, + {4.700000000e+01f, 5.614329680810343111535109717067692832e+18L}, + {4.800000000e+01f, 1.493630213112993142255380874118611276e+19L}, + {4.900000000e+01f, 3.975442747903744836006716995140844348e+19L}, + {5.000000000e+01f, 1.058563689713169096306154143322998720e+20L}}; +} + +template +void compareWithBoost(T RelTolerance, T AbsTolerance) { + const std::map Table = generateExpintTestValues(); + for (const auto& xy : Table) { + const T Value = std::expint(static_cast(xy.first)); + const T ExpectedValue = static_cast(xy.second); + const T Tolerance = std::abs(ExpectedValue) * RelTolerance + AbsTolerance; + const T Error = std::abs(Value - ExpectedValue); + assert(Error <= Tolerance); + } +} + +/// \brief Checks for a change of sign close to the value from the literature +/// for the root of expint. +/// +/// \tparam T one of float, double, or long double +/// \param [in] RelTolerance Tolerance for x; width of the interval where +/// expint changes its sign. +/// The reference value 0.372507... in its hexadecimal representation is taken +/// from Cody, W., & Thacher, H. (1969). +/// Chebyshev Approximations for the Exponential Integral Ei(x). +/// Mathematics of Computation, 23(106), 289-303. doi:10.2307/2004423 +template +void compareRootWithCody1969(T RelTolerance) { + const T Root = T(0x0.5F5CA54AD2D7F0F264C31p0L); + const T Tolerance = RelTolerance * Root; + assert(std::expint(Root - Tolerance) < T(0)); + assert(std::expint(Root + Tolerance) > T(0)); +} + +int main(int, char**) { + compareRootWithCody1969(1e-6f); + compareRootWithCody1969(1e-14); + + compareWithBoost(1e-6f, 0.0f); + compareWithBoost(1e-14, 0.0); +} Index: libcxx/test/std/numerics/c.math/c.math.expint/expint_limits.pass.cpp =================================================================== --- /dev/null +++ libcxx/test/std/numerics/c.math/c.math.expint/expint_limits.pass.cpp @@ -0,0 +1,40 @@ +//===----------------------------------------------------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +// UNSUPPORTED: c++98, c++03, c++11, c++14 + +// + +// floating_point expint (arithmetic x); // C++17 +// float expintf(float x); // C++17 +// long double expintl(long double x); // C++17 + +#include +#include +#include + +template +void testLimits() { + const T NegativeLimit = std::expint(-std::numeric_limits::infinity()); + assert(NegativeLimit <= 0); + assert(NegativeLimit >= -std::numeric_limits::min()); + + const T PositiveLimit = std::expint(std::numeric_limits::infinity()); + assert(PositiveLimit > T(0)); + assert(std::isinf(PositiveLimit)); + + const T ExpintAtZero = std::expint(T(0)); + assert(ExpintAtZero < T(0)); + assert(std::isinf(ExpintAtZero)); +} + +int main(int, char**) { + testLimits(); + testLimits(); + testLimits(); +} Index: libcxx/test/std/numerics/c.math/cmath.pass.cpp =================================================================== --- libcxx/test/std/numerics/c.math/cmath.pass.cpp +++ libcxx/test/std/numerics/c.math/cmath.pass.cpp @@ -100,6 +100,10 @@ Ambiguous tgamma(Ambiguous){ return Ambiguous(); } Ambiguous trunc(Ambiguous){ return Ambiguous(); } +#if _LIBCPP_STD_VER > 14 +Ambiguous expint(Ambiguous){ return Ambiguous(); } +#endif + template ()))> std::true_type has_abs_imp(int); template @@ -1553,6 +1557,28 @@ assert(std::trunc(1) == 1); } +#if _LIBCPP_STD_VER > 14 + +void test_expint() +{ + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); + static_assert((std::is_same::value), ""); +} + +#endif + int main(int, char**) { test_abs(); @@ -1626,5 +1652,9 @@ test_tgamma(); test_trunc(); +#if _LIBCPP_STD_VER > 14 + test_expint(); +#endif + return 0; }