1615 lines · plain
1// Copyright John Maddock 2007.2// Copyright Matt Borland 2024.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_MATH_EXPINT_HPP8#define BOOST_MATH_EXPINT_HPP9 10#ifdef _MSC_VER11#pragma once12#pragma warning(push)13#pragma warning(disable:4702) // Unreachable code (release mode only warning)14#endif15 16#include <boost/math/tools/config.hpp>17#include <boost/math/tools/cstdint.hpp>18#include <boost/math/tools/type_traits.hpp>19#include <boost/math/tools/tuple.hpp>20#include <boost/math/tools/precision.hpp>21#include <boost/math/tools/promotion.hpp>22#include <boost/math/tools/fraction.hpp>23#include <boost/math/tools/series.hpp>24#include <boost/math/policies/error_handling.hpp>25#include <boost/math/special_functions/math_fwd.hpp>26#include <boost/math/special_functions/digamma.hpp>27#include <boost/math/special_functions/log1p.hpp>28 29#if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)30//31// This is the only way we can avoid32// warning: non-standard suffix on floating constant [-Wpedantic]33// when building with -Wall -pedantic. Neither __extension__34// nor #pragma diagnostic ignored work :(35//36#pragma GCC system_header37#endif38 39namespace boost{ namespace math{40 41template <class T, class Policy>42BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type43 expint(unsigned n, T z, const Policy& /*pol*/);44 45namespace detail{46 47template <class T>48BOOST_MATH_GPU_ENABLED inline T expint_1_rational(const T& z, const boost::math::integral_constant<int, 0>&)49{50 // this function is never actually called51 BOOST_MATH_ASSERT(0);52 return z;53}54 55template <class T>56BOOST_MATH_GPU_ENABLED T expint_1_rational(const T& z, const boost::math::integral_constant<int, 53>&)57{58 BOOST_MATH_STD_USING59 T result;60 if(z <= 1)61 {62 // Maximum Deviation Found: 2.006e-1863 // Expected Error Term: 2.006e-1864 // Max error found at double precision: 2.760e-1765 // LCOV_EXCL_START66 static const T Y = 0.66373538970947265625F;67 static const T P[6] = {68 BOOST_MATH_BIG_CONSTANT(T, 53, 0.0865197248079397976498),69 BOOST_MATH_BIG_CONSTANT(T, 53, 0.0320913665303559189999),70 BOOST_MATH_BIG_CONSTANT(T, 53, -0.245088216639761496153),71 BOOST_MATH_BIG_CONSTANT(T, 53, -0.0368031736257943745142),72 BOOST_MATH_BIG_CONSTANT(T, 53, -0.00399167106081113256961),73 BOOST_MATH_BIG_CONSTANT(T, 53, -0.000111507792921197858394)74 };75 static const T Q[6] = {76 BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),77 BOOST_MATH_BIG_CONSTANT(T, 53, 0.37091387659397013215),78 BOOST_MATH_BIG_CONSTANT(T, 53, 0.056770677104207528384),79 BOOST_MATH_BIG_CONSTANT(T, 53, 0.00427347600017103698101),80 BOOST_MATH_BIG_CONSTANT(T, 53, 0.000131049900798434683324),81 BOOST_MATH_BIG_CONSTANT(T, 53, -0.528611029520217142048e-6)82 };83 // LCOV_EXCL_STOP84 result = tools::evaluate_polynomial(P, z)85 / tools::evaluate_polynomial(Q, z);86 result += z - log(z) - Y;87 }88 else if(z < -boost::math::tools::log_min_value<T>())89 {90 // Maximum Deviation Found (interpolated): 1.444e-1791 // Max error found at double precision: 3.119e-1792 // LCOV_EXCL_START93 static const T P[11] = {94 BOOST_MATH_BIG_CONSTANT(T, 53, -0.121013190657725568138e-18),95 BOOST_MATH_BIG_CONSTANT(T, 53, -0.999999999999998811143),96 BOOST_MATH_BIG_CONSTANT(T, 53, -43.3058660811817946037),97 BOOST_MATH_BIG_CONSTANT(T, 53, -724.581482791462469795),98 BOOST_MATH_BIG_CONSTANT(T, 53, -6046.8250112711035463),99 BOOST_MATH_BIG_CONSTANT(T, 53, -27182.6254466733970467),100 BOOST_MATH_BIG_CONSTANT(T, 53, -66598.2652345418633509),101 BOOST_MATH_BIG_CONSTANT(T, 53, -86273.1567711649528784),102 BOOST_MATH_BIG_CONSTANT(T, 53, -54844.4587226402067411),103 BOOST_MATH_BIG_CONSTANT(T, 53, -14751.4895786128450662),104 BOOST_MATH_BIG_CONSTANT(T, 53, -1185.45720315201027667)105 };106 static const T Q[12] = {107 BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),108 BOOST_MATH_BIG_CONSTANT(T, 53, 45.3058660811801465927),109 BOOST_MATH_BIG_CONSTANT(T, 53, 809.193214954550328455),110 BOOST_MATH_BIG_CONSTANT(T, 53, 7417.37624454689546708),111 BOOST_MATH_BIG_CONSTANT(T, 53, 38129.5594484818471461),112 BOOST_MATH_BIG_CONSTANT(T, 53, 113057.05869159631492),113 BOOST_MATH_BIG_CONSTANT(T, 53, 192104.047790227984431),114 BOOST_MATH_BIG_CONSTANT(T, 53, 180329.498380501819718),115 BOOST_MATH_BIG_CONSTANT(T, 53, 86722.3403467334749201),116 BOOST_MATH_BIG_CONSTANT(T, 53, 18455.4124737722049515),117 BOOST_MATH_BIG_CONSTANT(T, 53, 1229.20784182403048905),118 BOOST_MATH_BIG_CONSTANT(T, 53, -0.776491285282330997549)119 };120 // LCOV_EXCL_STOP121 T recip = 1 / z;122 result = 1 + tools::evaluate_polynomial(P, recip)123 / tools::evaluate_polynomial(Q, recip);124 result *= exp(-z) * recip;125 }126 else127 {128 result = 0;129 }130 return result;131}132 133template <class T>134BOOST_MATH_GPU_ENABLED T expint_1_rational(const T& z, const boost::math::integral_constant<int, 64>&)135{136 BOOST_MATH_STD_USING137 T result;138 if(z <= 1)139 {140 // Maximum Deviation Found: 3.807e-20141 // Expected Error Term: 3.807e-20142 // Max error found at long double precision: 6.249e-20143 // LCOV_EXCL_START144 static const T Y = 0.66373538970947265625F;145 static const T P[6] = {146 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0865197248079397956816),147 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0275114007037026844633),148 BOOST_MATH_BIG_CONSTANT(T, 64, -0.246594388074877139824),149 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0237624819878732642231),150 BOOST_MATH_BIG_CONSTANT(T, 64, -0.00259113319641673986276),151 BOOST_MATH_BIG_CONSTANT(T, 64, 0.30853660894346057053e-4)152 };153 static const T Q[7] = {154 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),155 BOOST_MATH_BIG_CONSTANT(T, 64, 0.317978365797784100273),156 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0393622602554758722511),157 BOOST_MATH_BIG_CONSTANT(T, 64, 0.00204062029115966323229),158 BOOST_MATH_BIG_CONSTANT(T, 64, 0.732512107100088047854e-5),159 BOOST_MATH_BIG_CONSTANT(T, 64, -0.202872781770207871975e-5),160 BOOST_MATH_BIG_CONSTANT(T, 64, 0.52779248094603709945e-7)161 };162 // LCOV_EXCL_STOP163 result = tools::evaluate_polynomial(P, z)164 / tools::evaluate_polynomial(Q, z);165 result += z - log(z) - Y;166 }167 else if(z < -boost::math::tools::log_min_value<T>())168 {169 // Maximum Deviation Found (interpolated): 2.220e-20170 // Max error found at long double precision: 1.346e-19171 // LCOV_EXCL_START172 static const T P[14] = {173 BOOST_MATH_BIG_CONSTANT(T, 64, -0.534401189080684443046e-23),174 BOOST_MATH_BIG_CONSTANT(T, 64, -0.999999999999999999905),175 BOOST_MATH_BIG_CONSTANT(T, 64, -62.1517806091379402505),176 BOOST_MATH_BIG_CONSTANT(T, 64, -1568.45688271895145277),177 BOOST_MATH_BIG_CONSTANT(T, 64, -21015.3431990874009619),178 BOOST_MATH_BIG_CONSTANT(T, 64, -164333.011755931661949),179 BOOST_MATH_BIG_CONSTANT(T, 64, -777917.270775426696103),180 BOOST_MATH_BIG_CONSTANT(T, 64, -2244188.56195255112937),181 BOOST_MATH_BIG_CONSTANT(T, 64, -3888702.98145335643429),182 BOOST_MATH_BIG_CONSTANT(T, 64, -3909822.65621952648353),183 BOOST_MATH_BIG_CONSTANT(T, 64, -2149033.9538897398457),184 BOOST_MATH_BIG_CONSTANT(T, 64, -584705.537139793925189),185 BOOST_MATH_BIG_CONSTANT(T, 64, -65815.2605361889477244),186 BOOST_MATH_BIG_CONSTANT(T, 64, -2038.82870680427258038)187 };188 static const T Q[14] = {189 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),190 BOOST_MATH_BIG_CONSTANT(T, 64, 64.1517806091379399478),191 BOOST_MATH_BIG_CONSTANT(T, 64, 1690.76044393722763785),192 BOOST_MATH_BIG_CONSTANT(T, 64, 24035.9534033068949426),193 BOOST_MATH_BIG_CONSTANT(T, 64, 203679.998633572361706),194 BOOST_MATH_BIG_CONSTANT(T, 64, 1074661.58459976978285),195 BOOST_MATH_BIG_CONSTANT(T, 64, 3586552.65020899358773),196 BOOST_MATH_BIG_CONSTANT(T, 64, 7552186.84989547621411),197 BOOST_MATH_BIG_CONSTANT(T, 64, 9853333.79353054111434),198 BOOST_MATH_BIG_CONSTANT(T, 64, 7689642.74550683631258),199 BOOST_MATH_BIG_CONSTANT(T, 64, 3385553.35146759180739),200 BOOST_MATH_BIG_CONSTANT(T, 64, 763218.072732396428725),201 BOOST_MATH_BIG_CONSTANT(T, 64, 73930.2995984054930821),202 BOOST_MATH_BIG_CONSTANT(T, 64, 2063.86994219629165937)203 };204 // LCOV_EXCL_STOP205 T recip = 1 / z;206 result = 1 + tools::evaluate_polynomial(P, recip)207 / tools::evaluate_polynomial(Q, recip);208 result *= exp(-z) * recip;209 }210 else211 {212 result = 0;213 }214 return result;215}216 217template <class T>218BOOST_MATH_GPU_ENABLED T expint_1_rational(const T& z, const boost::math::integral_constant<int, 113>&)219{220 BOOST_MATH_STD_USING221 T result;222 if(z <= 1)223 {224 // Maximum Deviation Found: 2.477e-35225 // Expected Error Term: 2.477e-35226 // Max error found at long double precision: 6.810e-35227 // LCOV_EXCL_START228 static const T Y = 0.66373538970947265625F;229 static const T P[10] = {230 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0865197248079397956434879099175975937),231 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0369066175910795772830865304506087759),232 BOOST_MATH_BIG_CONSTANT(T, 113, -0.24272036838415474665971599314725545),233 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0502166331248948515282379137550178307),234 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00768384138547489410285101483730424919),235 BOOST_MATH_BIG_CONSTANT(T, 113, -0.000612574337702109683505224915484717162),236 BOOST_MATH_BIG_CONSTANT(T, 113, -0.380207107950635046971492617061708534e-4),237 BOOST_MATH_BIG_CONSTANT(T, 113, -0.136528159460768830763009294683628406e-5),238 BOOST_MATH_BIG_CONSTANT(T, 113, -0.346839106212658259681029388908658618e-7),239 BOOST_MATH_BIG_CONSTANT(T, 113, -0.340500302777838063940402160594523429e-9)240 };241 static const T Q[10] = {242 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),243 BOOST_MATH_BIG_CONSTANT(T, 113, 0.426568827778942588160423015589537302),244 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0841384046470893490592450881447510148),245 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0100557215850668029618957359471132995),246 BOOST_MATH_BIG_CONSTANT(T, 113, 0.000799334870474627021737357294799839363),247 BOOST_MATH_BIG_CONSTANT(T, 113, 0.434452090903862735242423068552687688e-4),248 BOOST_MATH_BIG_CONSTANT(T, 113, 0.15829674748799079874182885081231252e-5),249 BOOST_MATH_BIG_CONSTANT(T, 113, 0.354406206738023762100882270033082198e-7),250 BOOST_MATH_BIG_CONSTANT(T, 113, 0.369373328141051577845488477377890236e-9),251 BOOST_MATH_BIG_CONSTANT(T, 113, -0.274149801370933606409282434677600112e-12)252 };253 // LCOV_EXCL_STOP254 result = tools::evaluate_polynomial(P, z)255 / tools::evaluate_polynomial(Q, z);256 result += z - log(z) - Y;257 }258 else if(z <= 4)259 {260 // Max error in interpolated form: 5.614e-35261 // Max error found at long double precision: 7.979e-35262 // LCOV_EXCL_START263 static const T Y = 0.70190334320068359375F;264 265 static const T P[16] = {266 BOOST_MATH_BIG_CONSTANT(T, 113, 0.298096656795020369955077350585959794),267 BOOST_MATH_BIG_CONSTANT(T, 113, 12.9314045995266142913135497455971247),268 BOOST_MATH_BIG_CONSTANT(T, 113, 226.144334921582637462526628217345501),269 BOOST_MATH_BIG_CONSTANT(T, 113, 2070.83670924261732722117682067381405),270 BOOST_MATH_BIG_CONSTANT(T, 113, 10715.1115684330959908244769731347186),271 BOOST_MATH_BIG_CONSTANT(T, 113, 30728.7876355542048019664777316053311),272 BOOST_MATH_BIG_CONSTANT(T, 113, 38520.6078609349855436936232610875297),273 BOOST_MATH_BIG_CONSTANT(T, 113, -27606.0780981527583168728339620565165),274 BOOST_MATH_BIG_CONSTANT(T, 113, -169026.485055785605958655247592604835),275 BOOST_MATH_BIG_CONSTANT(T, 113, -254361.919204983608659069868035092282),276 BOOST_MATH_BIG_CONSTANT(T, 113, -195765.706874132267953259272028679935),277 BOOST_MATH_BIG_CONSTANT(T, 113, -83352.6826013533205474990119962408675),278 BOOST_MATH_BIG_CONSTANT(T, 113, -19251.6828496869586415162597993050194),279 BOOST_MATH_BIG_CONSTANT(T, 113, -2226.64251774578542836725386936102339),280 BOOST_MATH_BIG_CONSTANT(T, 113, -109.009437301400845902228611986479816),281 BOOST_MATH_BIG_CONSTANT(T, 113, -1.51492042209561411434644938098833499)282 };283 static const T Q[16] = {284 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),285 BOOST_MATH_BIG_CONSTANT(T, 113, 46.734521442032505570517810766704587),286 BOOST_MATH_BIG_CONSTANT(T, 113, 908.694714348462269000247450058595655),287 BOOST_MATH_BIG_CONSTANT(T, 113, 9701.76053033673927362784882748513195),288 BOOST_MATH_BIG_CONSTANT(T, 113, 63254.2815292641314236625196594947774),289 BOOST_MATH_BIG_CONSTANT(T, 113, 265115.641285880437335106541757711092),290 BOOST_MATH_BIG_CONSTANT(T, 113, 732707.841188071900498536533086567735),291 BOOST_MATH_BIG_CONSTANT(T, 113, 1348514.02492635723327306628712057794),292 BOOST_MATH_BIG_CONSTANT(T, 113, 1649986.81455283047769673308781585991),293 BOOST_MATH_BIG_CONSTANT(T, 113, 1326000.828522976970116271208812099),294 BOOST_MATH_BIG_CONSTANT(T, 113, 683643.09490612171772350481773951341),295 BOOST_MATH_BIG_CONSTANT(T, 113, 217640.505137263607952365685653352229),296 BOOST_MATH_BIG_CONSTANT(T, 113, 40288.3467237411710881822569476155485),297 BOOST_MATH_BIG_CONSTANT(T, 113, 3932.89353979531632559232883283175754),298 BOOST_MATH_BIG_CONSTANT(T, 113, 169.845369689596739824177412096477219),299 BOOST_MATH_BIG_CONSTANT(T, 113, 2.17607292280092201170768401876895354)300 };301 // LCOV_EXCL_STOP302 T recip = 1 / z;303 result = Y + tools::evaluate_polynomial(P, recip)304 / tools::evaluate_polynomial(Q, recip);305 result *= exp(-z) * recip;306 }307 else if(z < -boost::math::tools::log_min_value<T>())308 {309 // Max error in interpolated form: 4.413e-35310 // Max error found at long double precision: 8.928e-35311 // LCOV_EXCL_START312 static const T P[19] = {313 BOOST_MATH_BIG_CONSTANT(T, 113, -0.559148411832951463689610809550083986e-40),314 BOOST_MATH_BIG_CONSTANT(T, 113, -0.999999999999999999999999999999999997),315 BOOST_MATH_BIG_CONSTANT(T, 113, -166.542326331163836642960118190147367),316 BOOST_MATH_BIG_CONSTANT(T, 113, -12204.639128796330005065904675153652),317 BOOST_MATH_BIG_CONSTANT(T, 113, -520807.069767086071806275022036146855),318 BOOST_MATH_BIG_CONSTANT(T, 113, -14435981.5242137970691490903863125326),319 BOOST_MATH_BIG_CONSTANT(T, 113, -274574945.737064301247496460758654196),320 BOOST_MATH_BIG_CONSTANT(T, 113, -3691611582.99810039356254671781473079),321 BOOST_MATH_BIG_CONSTANT(T, 113, -35622515944.8255047299363690814678763),322 BOOST_MATH_BIG_CONSTANT(T, 113, -248040014774.502043161750715548451142),323 BOOST_MATH_BIG_CONSTANT(T, 113, -1243190389769.53458416330946622607913),324 BOOST_MATH_BIG_CONSTANT(T, 113, -4441730126135.54739052731990368425339),325 BOOST_MATH_BIG_CONSTANT(T, 113, -11117043181899.7388524310281751971366),326 BOOST_MATH_BIG_CONSTANT(T, 113, -18976497615396.9717776601813519498961),327 BOOST_MATH_BIG_CONSTANT(T, 113, -21237496819711.1011661104761906067131),328 BOOST_MATH_BIG_CONSTANT(T, 113, -14695899122092.5161620333466757812848),329 BOOST_MATH_BIG_CONSTANT(T, 113, -5737221535080.30569711574295785864903),330 BOOST_MATH_BIG_CONSTANT(T, 113, -1077042281708.42654526404581272546244),331 BOOST_MATH_BIG_CONSTANT(T, 113, -68028222642.1941480871395695677675137)332 };333 static const T Q[20] = {334 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),335 BOOST_MATH_BIG_CONSTANT(T, 113, 168.542326331163836642960118190147311),336 BOOST_MATH_BIG_CONSTANT(T, 113, 12535.7237814586576783518249115343619),337 BOOST_MATH_BIG_CONSTANT(T, 113, 544891.263372016404143120911148640627),338 BOOST_MATH_BIG_CONSTANT(T, 113, 15454474.7241010258634446523045237762),339 BOOST_MATH_BIG_CONSTANT(T, 113, 302495899.896629522673410325891717381),340 BOOST_MATH_BIG_CONSTANT(T, 113, 4215565948.38886507646911672693270307),341 BOOST_MATH_BIG_CONSTANT(T, 113, 42552409471.7951815668506556705733344),342 BOOST_MATH_BIG_CONSTANT(T, 113, 313592377066.753173979584098301610186),343 BOOST_MATH_BIG_CONSTANT(T, 113, 1688763640223.4541980740597514904542),344 BOOST_MATH_BIG_CONSTANT(T, 113, 6610992294901.59589748057620192145704),345 BOOST_MATH_BIG_CONSTANT(T, 113, 18601637235659.6059890851321772682606),346 BOOST_MATH_BIG_CONSTANT(T, 113, 36944278231087.2571020964163402941583),347 BOOST_MATH_BIG_CONSTANT(T, 113, 50425858518481.7497071917028793820058),348 BOOST_MATH_BIG_CONSTANT(T, 113, 45508060902865.0899967797848815980644),349 BOOST_MATH_BIG_CONSTANT(T, 113, 25649955002765.3817331501988304758142),350 BOOST_MATH_BIG_CONSTANT(T, 113, 8259575619094.6518520988612711292331),351 BOOST_MATH_BIG_CONSTANT(T, 113, 1299981487496.12607474362723586264515),352 BOOST_MATH_BIG_CONSTANT(T, 113, 70242279152.8241187845178443118302693),353 BOOST_MATH_BIG_CONSTANT(T, 113, -37633302.9409263839042721539363416685)354 };355 // LCOV_EXCL_STOP356 T recip = 1 / z;357 result = 1 + tools::evaluate_polynomial(P, recip)358 / tools::evaluate_polynomial(Q, recip);359 result *= exp(-z) * recip;360 }361 else362 {363 result = 0;364 }365 return result;366}367 368 369template <class T>370struct expint_fraction371{372 typedef boost::math::pair<T,T> result_type;373 BOOST_MATH_GPU_ENABLED expint_fraction(unsigned n_, T z_) : b(n_ + z_), i(-1), n(n_){}374 BOOST_MATH_GPU_ENABLED boost::math::pair<T,T> operator()()375 {376 boost::math::pair<T,T> result = boost::math::make_pair(-static_cast<T>((i+1) * (n+i)), b);377 b += 2;378 ++i;379 return result;380 }381private:382 T b;383 int i;384 unsigned n;385};386 387template <class T, class Policy>388BOOST_MATH_GPU_ENABLED inline T expint_as_fraction(unsigned n, T z, const Policy& pol)389{390 BOOST_MATH_STD_USING391 BOOST_MATH_INSTRUMENT_VARIABLE(z)392 boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();393 expint_fraction<T> f(n, z);394 T result = tools::continued_fraction_b(395 f,396 boost::math::policies::get_epsilon<T, Policy>(),397 max_iter);398 policies::check_series_iterations<T>("boost::math::expint_continued_fraction<%1%>(unsigned,%1%)", max_iter, pol);399 BOOST_MATH_INSTRUMENT_VARIABLE(result)400 BOOST_MATH_INSTRUMENT_VARIABLE(max_iter)401 result = exp(-z) / result;402 BOOST_MATH_INSTRUMENT_VARIABLE(result)403 return result;404}405 406template <class T>407struct expint_series408{409 typedef T result_type;410 BOOST_MATH_GPU_ENABLED expint_series(unsigned k_, T z_, T x_k_, T denom_, T fact_)411 : k(k_), z(z_), x_k(x_k_), denom(denom_), fact(fact_){}412 BOOST_MATH_GPU_ENABLED T operator()()413 {414 x_k *= -z;415 denom += 1;416 fact *= ++k;417 return x_k / (denom * fact);418 }419private:420 unsigned k;421 T z;422 T x_k;423 T denom;424 T fact;425};426 427template <class T, class Policy>428BOOST_MATH_GPU_ENABLED inline T expint_as_series(unsigned n, T z, const Policy& pol)429{430 BOOST_MATH_STD_USING431 boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();432 433 BOOST_MATH_INSTRUMENT_VARIABLE(z)434 435 T result = 0;436 T x_k = -1;437 T denom = T(1) - n;438 T fact = 1;439 unsigned k = 0;440 for(; k < n - 1;)441 {442 result += x_k / (denom * fact);443 denom += 1;444 x_k *= -z;445 fact *= ++k;446 }447 BOOST_MATH_INSTRUMENT_VARIABLE(result)448 result += pow(-z, static_cast<T>(n - 1))449 * (boost::math::digamma(static_cast<T>(n), pol) - log(z)) / fact;450 BOOST_MATH_INSTRUMENT_VARIABLE(result)451 452 expint_series<T> s(k, z, x_k, denom, fact);453 result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, result);454 policies::check_series_iterations<T>("boost::math::expint_series<%1%>(unsigned,%1%)", max_iter, pol);455 BOOST_MATH_INSTRUMENT_VARIABLE(result)456 BOOST_MATH_INSTRUMENT_VARIABLE(max_iter)457 return result;458}459 460template <class T, class Policy, class Tag>461BOOST_MATH_GPU_ENABLED T expint_imp(unsigned n, T z, const Policy& pol, const Tag& tag)462{463 BOOST_MATH_STD_USING464 constexpr auto function = "boost::math::expint<%1%>(unsigned, %1%)";465 if(z < 0)466 return policies::raise_domain_error<T>(function, "Function requires z >= 0 but got %1%.", z, pol);467 if(z == 0)468 return n == 1 ? policies::raise_overflow_error<T>(function, nullptr, pol) : T(1 / (static_cast<T>(n - 1)));469 470 T result;471 472 bool f;473 if(n < 3)474 {475 f = z < T(0.5);476 }477 else478 {479 f = z < (static_cast<T>(n - 2) / static_cast<T>(n - 1));480 }481#ifdef _MSC_VER482# pragma warning(push)483# pragma warning(disable:4127) // conditional expression is constant484#endif485 if(n == 0)486 {487 result = exp(-z) / z;488 }489 else if((n == 1) && (Tag::value))490 {491 result = expint_1_rational(z, tag);492 }493 else if(f)494 {495 result = expint_as_series(n, z, pol);496 }497 else498 {499 result = expint_as_fraction(n, z, pol);500 }501#ifdef _MSC_VER502# pragma warning(pop)503#endif504 505 return result;506}507 508template <class T>509struct expint_i_series510{511 typedef T result_type;512 BOOST_MATH_GPU_ENABLED expint_i_series(T z_) : k(0), z_k(1), z(z_){}513 BOOST_MATH_GPU_ENABLED T operator()()514 {515 z_k *= z / ++k;516 return z_k / k;517 }518private:519 unsigned k;520 T z_k;521 T z;522};523 524template <class T, class Policy>525BOOST_MATH_GPU_ENABLED T expint_i_as_series(T z, const Policy& pol)526{527 BOOST_MATH_STD_USING528 T result = log(z); // (log(z) - log(1 / z)) / 2;529 result += constants::euler<T>();530 expint_i_series<T> s(z);531 boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();532 result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, result);533 policies::check_series_iterations<T>("boost::math::expint_i_series<%1%>(%1%)", max_iter, pol);534 return result;535}536 537template <class T, class Policy, class Tag>538BOOST_MATH_GPU_ENABLED T expint_i_imp(T z, const Policy& pol, const Tag& tag)539{540 constexpr auto function = "boost::math::expint<%1%>(%1%)";541 if(z < 0)542 return -expint_imp(1, T(-z), pol, tag);543 if(z == 0)544 return -policies::raise_overflow_error<T>(function, nullptr, pol); // LCOV_EXCL_LINE confirmed covered by real_concept tests545 return expint_i_as_series(z, pol);546}547 548template <class T, class Policy>549BOOST_MATH_GPU_ENABLED T expint_i_imp(T z, const Policy& pol, const boost::math::integral_constant<int, 53>& tag)550{551 BOOST_MATH_STD_USING552 constexpr auto function = "boost::math::expint<%1%>(%1%)";553 if(z < 0)554 return -expint_imp(1, T(-z), pol, tag);555 if(z == 0)556 return -policies::raise_overflow_error<T>(function, nullptr, pol);557 558 T result;559 560 if(z <= 6)561 {562 // Maximum Deviation Found: 2.852e-18563 // Expected Error Term: 2.852e-18564 // Max Error found at double precision = Poly: 2.636335e-16 Cheb: 4.187027e-16565 // LCOV_EXCL_START566 BOOST_MATH_STATIC const T P[10] = {567 BOOST_MATH_BIG_CONSTANT(T, 53, 2.98677224343598593013),568 BOOST_MATH_BIG_CONSTANT(T, 53, 0.356343618769377415068),569 BOOST_MATH_BIG_CONSTANT(T, 53, 0.780836076283730801839),570 BOOST_MATH_BIG_CONSTANT(T, 53, 0.114670926327032002811),571 BOOST_MATH_BIG_CONSTANT(T, 53, 0.0499434773576515260534),572 BOOST_MATH_BIG_CONSTANT(T, 53, 0.00726224593341228159561),573 BOOST_MATH_BIG_CONSTANT(T, 53, 0.00115478237227804306827),574 BOOST_MATH_BIG_CONSTANT(T, 53, 0.000116419523609765200999),575 BOOST_MATH_BIG_CONSTANT(T, 53, 0.798296365679269702435e-5),576 BOOST_MATH_BIG_CONSTANT(T, 53, 0.2777056254402008721e-6)577 };578 BOOST_MATH_STATIC const T Q[8] = {579 BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),580 BOOST_MATH_BIG_CONSTANT(T, 53, -1.17090412365413911947),581 BOOST_MATH_BIG_CONSTANT(T, 53, 0.62215109846016746276),582 BOOST_MATH_BIG_CONSTANT(T, 53, -0.195114782069495403315),583 BOOST_MATH_BIG_CONSTANT(T, 53, 0.0391523431392967238166),584 BOOST_MATH_BIG_CONSTANT(T, 53, -0.00504800158663705747345),585 BOOST_MATH_BIG_CONSTANT(T, 53, 0.000389034007436065401822),586 BOOST_MATH_BIG_CONSTANT(T, 53, -0.138972589601781706598e-4)587 };588 589 BOOST_MATH_STATIC_LOCAL_VARIABLE const T c1 = BOOST_MATH_BIG_CONSTANT(T, 53, 1677624236387711.0);590 BOOST_MATH_STATIC_LOCAL_VARIABLE const T c2 = BOOST_MATH_BIG_CONSTANT(T, 53, 4503599627370496.0);591 BOOST_MATH_STATIC_LOCAL_VARIABLE const T r1 = static_cast<T>(c1 / c2);592 BOOST_MATH_STATIC_LOCAL_VARIABLE const T r2 = BOOST_MATH_BIG_CONSTANT(T, 53, 0.131401834143860282009280387409357165515556574352422001206362e-16);593 BOOST_MATH_STATIC_LOCAL_VARIABLE const T r = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 53, 0.372507410781366634461991866580119133535689497771654051555657435242200120636201854384926049951548942392));594 // LCOV_EXCL_STOP595 T t = (z / 3) - 1;596 result = tools::evaluate_polynomial(P, t)597 / tools::evaluate_polynomial(Q, t);598 t = (z - r1) - r2;599 result *= t;600 if(fabs(t) < T(0.1))601 {602 result += boost::math::log1p(t / r, pol);603 }604 else605 {606 result += log(z / r);607 }608 }609 else if (z <= 10)610 {611 // Maximum Deviation Found: 6.546e-17612 // Expected Error Term: 6.546e-17613 // Max Error found at double precision = Poly: 6.890169e-17 Cheb: 6.772128e-17614 // LCOV_EXCL_START615 BOOST_MATH_STATIC_LOCAL_VARIABLE const T Y = 1.158985137939453125F;616 BOOST_MATH_STATIC const T P[8] = {617 BOOST_MATH_BIG_CONSTANT(T, 53, 0.00139324086199402804173),618 BOOST_MATH_BIG_CONSTANT(T, 53, -0.0349921221823888744966),619 BOOST_MATH_BIG_CONSTANT(T, 53, -0.0264095520754134848538),620 BOOST_MATH_BIG_CONSTANT(T, 53, -0.00761224003005476438412),621 BOOST_MATH_BIG_CONSTANT(T, 53, -0.00247496209592143627977),622 BOOST_MATH_BIG_CONSTANT(T, 53, -0.000374885917942100256775),623 BOOST_MATH_BIG_CONSTANT(T, 53, -0.554086272024881826253e-4),624 BOOST_MATH_BIG_CONSTANT(T, 53, -0.396487648924804510056e-5)625 };626 BOOST_MATH_STATIC const T Q[8] = {627 BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),628 BOOST_MATH_BIG_CONSTANT(T, 53, 0.744625566823272107711),629 BOOST_MATH_BIG_CONSTANT(T, 53, 0.329061095011767059236),630 BOOST_MATH_BIG_CONSTANT(T, 53, 0.100128624977313872323),631 BOOST_MATH_BIG_CONSTANT(T, 53, 0.0223851099128506347278),632 BOOST_MATH_BIG_CONSTANT(T, 53, 0.00365334190742316650106),633 BOOST_MATH_BIG_CONSTANT(T, 53, 0.000402453408512476836472),634 BOOST_MATH_BIG_CONSTANT(T, 53, 0.263649630720255691787e-4)635 };636 // LCOV_EXCL_STOP637 T t = z / 2 - 4;638 result = Y + tools::evaluate_polynomial(P, t)639 / tools::evaluate_polynomial(Q, t);640 result *= exp(z) / z;641 result += z;642 }643 else if(z <= 20)644 {645 // Maximum Deviation Found: 1.843e-17646 // Expected Error Term: -1.842e-17647 // Max Error found at double precision = Poly: 4.375868e-17 Cheb: 5.860967e-17648 // LCOV_EXCL_START649 BOOST_MATH_STATIC_LOCAL_VARIABLE const T Y = 1.0869731903076171875F;650 BOOST_MATH_STATIC const T P[9] = {651 BOOST_MATH_BIG_CONSTANT(T, 53, -0.00893891094356945667451),652 BOOST_MATH_BIG_CONSTANT(T, 53, -0.0484607730127134045806),653 BOOST_MATH_BIG_CONSTANT(T, 53, -0.0652810444222236895772),654 BOOST_MATH_BIG_CONSTANT(T, 53, -0.0478447572647309671455),655 BOOST_MATH_BIG_CONSTANT(T, 53, -0.0226059218923777094596),656 BOOST_MATH_BIG_CONSTANT(T, 53, -0.00720603636917482065907),657 BOOST_MATH_BIG_CONSTANT(T, 53, -0.00155941947035972031334),658 BOOST_MATH_BIG_CONSTANT(T, 53, -0.000209750022660200888349),659 BOOST_MATH_BIG_CONSTANT(T, 53, -0.138652200349182596186e-4)660 };661 BOOST_MATH_STATIC const T Q[9] = {662 BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),663 BOOST_MATH_BIG_CONSTANT(T, 53, 1.97017214039061194971),664 BOOST_MATH_BIG_CONSTANT(T, 53, 1.86232465043073157508),665 BOOST_MATH_BIG_CONSTANT(T, 53, 1.09601437090337519977),666 BOOST_MATH_BIG_CONSTANT(T, 53, 0.438873285773088870812),667 BOOST_MATH_BIG_CONSTANT(T, 53, 0.122537731979686102756),668 BOOST_MATH_BIG_CONSTANT(T, 53, 0.0233458478275769288159),669 BOOST_MATH_BIG_CONSTANT(T, 53, 0.00278170769163303669021),670 BOOST_MATH_BIG_CONSTANT(T, 53, 0.000159150281166108755531)671 };672 // LCOV_EXCL_STOP673 T t = z / 5 - 3;674 result = Y + tools::evaluate_polynomial(P, t)675 / tools::evaluate_polynomial(Q, t);676 result *= exp(z) / z;677 result += z;678 }679 else if(z <= 40)680 {681 // Maximum Deviation Found: 5.102e-18682 // Expected Error Term: 5.101e-18683 // Max Error found at double precision = Poly: 1.441088e-16 Cheb: 1.864792e-16684 // LCOV_EXCL_START685 BOOST_MATH_STATIC_LOCAL_VARIABLE const T Y = 1.03937530517578125F;686 BOOST_MATH_STATIC const T P[9] = {687 BOOST_MATH_BIG_CONSTANT(T, 53, -0.00356165148914447597995),688 BOOST_MATH_BIG_CONSTANT(T, 53, -0.0229930320357982333406),689 BOOST_MATH_BIG_CONSTANT(T, 53, -0.0449814350482277917716),690 BOOST_MATH_BIG_CONSTANT(T, 53, -0.0453759383048193402336),691 BOOST_MATH_BIG_CONSTANT(T, 53, -0.0272050837209380717069),692 BOOST_MATH_BIG_CONSTANT(T, 53, -0.00994403059883350813295),693 BOOST_MATH_BIG_CONSTANT(T, 53, -0.00207592267812291726961),694 BOOST_MATH_BIG_CONSTANT(T, 53, -0.000192178045857733706044),695 BOOST_MATH_BIG_CONSTANT(T, 53, -0.113161784705911400295e-9)696 };697 BOOST_MATH_STATIC const T Q[9] = {698 BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),699 BOOST_MATH_BIG_CONSTANT(T, 53, 2.84354408840148561131),700 BOOST_MATH_BIG_CONSTANT(T, 53, 3.6599610090072393012),701 BOOST_MATH_BIG_CONSTANT(T, 53, 2.75088464344293083595),702 BOOST_MATH_BIG_CONSTANT(T, 53, 1.2985244073998398643),703 BOOST_MATH_BIG_CONSTANT(T, 53, 0.383213198510794507409),704 BOOST_MATH_BIG_CONSTANT(T, 53, 0.0651165455496281337831),705 BOOST_MATH_BIG_CONSTANT(T, 53, 0.00488071077519227853585)706 };707 // LCOV_EXCL_STOP708 T t = z / 10 - 3;709 result = Y + tools::evaluate_polynomial(P, t)710 / tools::evaluate_polynomial(Q, t);711 result *= exp(z) / z;712 result += z;713 }714 else715 {716 // Max Error found at double precision = 3.381886e-17717 // LCOV_EXCL_START718 BOOST_MATH_STATIC_LOCAL_VARIABLE const T exp40 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 53, 2.35385266837019985407899910749034804508871617254555467236651e17));719 BOOST_MATH_STATIC_LOCAL_VARIABLE const T Y= 1.013065338134765625F;720 BOOST_MATH_STATIC const T P[6] = {721 BOOST_MATH_BIG_CONSTANT(T, 53, -0.0130653381347656243849),722 BOOST_MATH_BIG_CONSTANT(T, 53, 0.19029710559486576682),723 BOOST_MATH_BIG_CONSTANT(T, 53, 94.7365094537197236011),724 BOOST_MATH_BIG_CONSTANT(T, 53, -2516.35323679844256203),725 BOOST_MATH_BIG_CONSTANT(T, 53, 18932.0850014925993025),726 BOOST_MATH_BIG_CONSTANT(T, 53, -38703.1431362056714134)727 };728 BOOST_MATH_STATIC const T Q[7] = {729 BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),730 BOOST_MATH_BIG_CONSTANT(T, 53, 61.9733592849439884145),731 BOOST_MATH_BIG_CONSTANT(T, 53, -2354.56211323420194283),732 BOOST_MATH_BIG_CONSTANT(T, 53, 22329.1459489893079041),733 BOOST_MATH_BIG_CONSTANT(T, 53, -70126.245140396567133),734 BOOST_MATH_BIG_CONSTANT(T, 53, 54738.2833147775537106),735 BOOST_MATH_BIG_CONSTANT(T, 53, 8297.16296356518409347)736 };737 // LCOV_EXCL_STOP738 T t = 1 / z;739 result = Y + tools::evaluate_polynomial(P, t)740 / tools::evaluate_polynomial(Q, t);741 if(z < 41)742 result *= exp(z) / z;743 else744 {745 // Avoid premature overflow if we can:746 t = z - 40;747 if(t > tools::log_max_value<T>())748 {749 result = policies::raise_overflow_error<T>(function, nullptr, pol);750 }751 else752 {753 result *= exp(z - 40) / z;754 if(result > tools::max_value<T>() / exp40)755 {756 result = policies::raise_overflow_error<T>(function, nullptr, pol);757 }758 else759 {760 result *= exp40;761 }762 }763 }764 result += z;765 }766 return result;767}768 769template <class T, class Policy>770BOOST_MATH_GPU_ENABLED T expint_i_imp(T z, const Policy& pol, const boost::math::integral_constant<int, 64>& tag)771{772 BOOST_MATH_STD_USING773 constexpr auto function = "boost::math::expint<%1%>(%1%)";774 if(z < 0)775 return -expint_imp(1, T(-z), pol, tag);776 if(z == 0)777 return -policies::raise_overflow_error<T>(function, nullptr, pol);778 779 T result;780 781 if(z <= 6)782 {783 // Maximum Deviation Found: 3.883e-21784 // Expected Error Term: 3.883e-21785 // Max Error found at long double precision = Poly: 3.344801e-19 Cheb: 4.989937e-19786 787 // LCOV_EXCL_START788 static const T P[11] = {789 BOOST_MATH_BIG_CONSTANT(T, 64, 2.98677224343598593764),790 BOOST_MATH_BIG_CONSTANT(T, 64, 0.25891613550886736592),791 BOOST_MATH_BIG_CONSTANT(T, 64, 0.789323584998672832285),792 BOOST_MATH_BIG_CONSTANT(T, 64, 0.092432587824602399339),793 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0514236978728625906656),794 BOOST_MATH_BIG_CONSTANT(T, 64, 0.00658477469745132977921),795 BOOST_MATH_BIG_CONSTANT(T, 64, 0.00124914538197086254233),796 BOOST_MATH_BIG_CONSTANT(T, 64, 0.000131429679565472408551),797 BOOST_MATH_BIG_CONSTANT(T, 64, 0.11293331317982763165e-4),798 BOOST_MATH_BIG_CONSTANT(T, 64, 0.629499283139417444244e-6),799 BOOST_MATH_BIG_CONSTANT(T, 64, 0.177833045143692498221e-7)800 };801 static const T Q[9] = {802 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),803 BOOST_MATH_BIG_CONSTANT(T, 64, -1.20352377969742325748),804 BOOST_MATH_BIG_CONSTANT(T, 64, 0.66707904942606479811),805 BOOST_MATH_BIG_CONSTANT(T, 64, -0.223014531629140771914),806 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0493340022262908008636),807 BOOST_MATH_BIG_CONSTANT(T, 64, -0.00741934273050807310677),808 BOOST_MATH_BIG_CONSTANT(T, 64, 0.00074353567782087939294),809 BOOST_MATH_BIG_CONSTANT(T, 64, -0.455861727069603367656e-4),810 BOOST_MATH_BIG_CONSTANT(T, 64, 0.131515429329812837701e-5)811 };812 813 static const T c1 = BOOST_MATH_BIG_CONSTANT(T, 64, 1677624236387711.0);814 static const T c2 = BOOST_MATH_BIG_CONSTANT(T, 64, 4503599627370496.0);815 static const T r1 = c1 / c2;816 static const T r2 = BOOST_MATH_BIG_CONSTANT(T, 64, 0.131401834143860282009280387409357165515556574352422001206362e-16);817 static const T r = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.372507410781366634461991866580119133535689497771654051555657435242200120636201854384926049951548942392));818 // LCOV_EXCL_STOP819 820 T t = (z / 3) - 1;821 result = tools::evaluate_polynomial(P, t)822 / tools::evaluate_polynomial(Q, t);823 t = (z - r1) - r2;824 result *= t;825 if(fabs(t) < T(0.1))826 {827 result += boost::math::log1p(t / r, pol);828 }829 else830 {831 result += log(z / r);832 }833 }834 else if (z <= 10)835 {836 // Maximum Deviation Found: 2.622e-21837 // Expected Error Term: -2.622e-21838 // Max Error found at long double precision = Poly: 1.208328e-20 Cheb: 1.073723e-20839 // LCOV_EXCL_START840 static const T Y = 1.158985137939453125F;841 static const T P[9] = {842 BOOST_MATH_BIG_CONSTANT(T, 64, 0.00139324086199409049399),843 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0345238388952337563247),844 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0382065278072592940767),845 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0156117003070560727392),846 BOOST_MATH_BIG_CONSTANT(T, 64, -0.00383276012430495387102),847 BOOST_MATH_BIG_CONSTANT(T, 64, -0.000697070540945496497992),848 BOOST_MATH_BIG_CONSTANT(T, 64, -0.877310384591205930343e-4),849 BOOST_MATH_BIG_CONSTANT(T, 64, -0.623067256376494930067e-5),850 BOOST_MATH_BIG_CONSTANT(T, 64, -0.377246883283337141444e-6)851 };852 static const T Q[10] = {853 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),854 BOOST_MATH_BIG_CONSTANT(T, 64, 1.08073635708902053767),855 BOOST_MATH_BIG_CONSTANT(T, 64, 0.553681133533942532909),856 BOOST_MATH_BIG_CONSTANT(T, 64, 0.176763647137553797451),857 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0387891748253869928121),858 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0060603004848394727017),859 BOOST_MATH_BIG_CONSTANT(T, 64, 0.000670519492939992806051),860 BOOST_MATH_BIG_CONSTANT(T, 64, 0.4947357050100855646e-4),861 BOOST_MATH_BIG_CONSTANT(T, 64, 0.204339282037446434827e-5),862 BOOST_MATH_BIG_CONSTANT(T, 64, 0.146951181174930425744e-7)863 };864 // LCOV_EXCL_STOP865 T t = z / 2 - 4;866 result = Y + tools::evaluate_polynomial(P, t)867 / tools::evaluate_polynomial(Q, t);868 result *= exp(z) / z;869 result += z;870 }871 else if(z <= 20)872 {873 // Maximum Deviation Found: 3.220e-20874 // Expected Error Term: 3.220e-20875 // Max Error found at long double precision = Poly: 7.696841e-20 Cheb: 6.205163e-20876 877 // LCOV_EXCL_START878 static const T Y = 1.0869731903076171875F;879 static const T P[10] = {880 BOOST_MATH_BIG_CONSTANT(T, 64, -0.00893891094356946995368),881 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0487562980088748775943),882 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0670568657950041926085),883 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0509577352851442932713),884 BOOST_MATH_BIG_CONSTANT(T, 64, -0.02551800927409034206),885 BOOST_MATH_BIG_CONSTANT(T, 64, -0.00892913759760086687083),886 BOOST_MATH_BIG_CONSTANT(T, 64, -0.00224469630207344379888),887 BOOST_MATH_BIG_CONSTANT(T, 64, -0.000392477245911296982776),888 BOOST_MATH_BIG_CONSTANT(T, 64, -0.44424044184395578775e-4),889 BOOST_MATH_BIG_CONSTANT(T, 64, -0.252788029251437017959e-5)890 };891 static const T Q[10] = {892 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),893 BOOST_MATH_BIG_CONSTANT(T, 64, 2.00323265503572414261),894 BOOST_MATH_BIG_CONSTANT(T, 64, 1.94688958187256383178),895 BOOST_MATH_BIG_CONSTANT(T, 64, 1.19733638134417472296),896 BOOST_MATH_BIG_CONSTANT(T, 64, 0.513137726038353385661),897 BOOST_MATH_BIG_CONSTANT(T, 64, 0.159135395578007264547),898 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0358233587351620919881),899 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0056716655597009417875),900 BOOST_MATH_BIG_CONSTANT(T, 64, 0.000577048986213535829925),901 BOOST_MATH_BIG_CONSTANT(T, 64, 0.290976943033493216793e-4)902 };903 // LCOV_EXCL_STOP904 T t = z / 5 - 3;905 result = Y + tools::evaluate_polynomial(P, t)906 / tools::evaluate_polynomial(Q, t);907 result *= exp(z) / z;908 result += z;909 }910 else if(z <= 40)911 {912 // Maximum Deviation Found: 2.940e-21913 // Expected Error Term: -2.938e-21914 // Max Error found at long double precision = Poly: 3.419893e-19 Cheb: 3.359874e-19915 // LCOV_EXCL_START916 static const T Y = 1.03937530517578125F;917 static const T P[12] = {918 BOOST_MATH_BIG_CONSTANT(T, 64, -0.00356165148914447278177),919 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0240235006148610849678),920 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0516699967278057976119),921 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0586603078706856245674),922 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0409960120868776180825),923 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0185485073689590665153),924 BOOST_MATH_BIG_CONSTANT(T, 64, -0.00537842101034123222417),925 BOOST_MATH_BIG_CONSTANT(T, 64, -0.000920988084778273760609),926 BOOST_MATH_BIG_CONSTANT(T, 64, -0.716742618812210980263e-4),927 BOOST_MATH_BIG_CONSTANT(T, 64, -0.504623302166487346677e-9),928 BOOST_MATH_BIG_CONSTANT(T, 64, 0.712662196671896837736e-10),929 BOOST_MATH_BIG_CONSTANT(T, 64, -0.533769629702262072175e-11)930 };931 static const T Q[9] = {932 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),933 BOOST_MATH_BIG_CONSTANT(T, 64, 3.13286733695729715455),934 BOOST_MATH_BIG_CONSTANT(T, 64, 4.49281223045653491929),935 BOOST_MATH_BIG_CONSTANT(T, 64, 3.84900294427622911374),936 BOOST_MATH_BIG_CONSTANT(T, 64, 2.15205199043580378211),937 BOOST_MATH_BIG_CONSTANT(T, 64, 0.802912186540269232424),938 BOOST_MATH_BIG_CONSTANT(T, 64, 0.194793170017818925388),939 BOOST_MATH_BIG_CONSTANT(T, 64, 0.0280128013584653182994),940 BOOST_MATH_BIG_CONSTANT(T, 64, 0.00182034930799902922549)941 };942 // LCOV_EXCL_STOP943 T t = z / 10 - 3;944 result = Y + tools::evaluate_polynomial(P, t)945 / tools::evaluate_polynomial(Q, t);946 BOOST_MATH_INSTRUMENT_VARIABLE(result)947 result *= exp(z) / z;948 BOOST_MATH_INSTRUMENT_VARIABLE(result)949 result += z;950 BOOST_MATH_INSTRUMENT_VARIABLE(result)951 }952 else953 {954 // Maximum Deviation Found: 3.536e-20955 // Max Error found at long double precision = Poly: 1.310671e-19 Cheb: 8.630943e-11956 // LCOV_EXCL_START957 static const T exp40 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.35385266837019985407899910749034804508871617254555467236651e17));958 static const T Y= 1.013065338134765625F;959 static const T P[9] = {960 BOOST_MATH_BIG_CONSTANT(T, 64, -0.0130653381347656250004),961 BOOST_MATH_BIG_CONSTANT(T, 64, 0.644487780349757303739),962 BOOST_MATH_BIG_CONSTANT(T, 64, 143.995670348227433964),963 BOOST_MATH_BIG_CONSTANT(T, 64, -13918.9322758014173709),964 BOOST_MATH_BIG_CONSTANT(T, 64, 476260.975133624194484),965 BOOST_MATH_BIG_CONSTANT(T, 64, -7437102.15135982802122),966 BOOST_MATH_BIG_CONSTANT(T, 64, 53732298.8764767916542),967 BOOST_MATH_BIG_CONSTANT(T, 64, -160695051.957997452509),968 BOOST_MATH_BIG_CONSTANT(T, 64, 137839271.592778020028)969 };970 static const T Q[9] = {971 BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),972 BOOST_MATH_BIG_CONSTANT(T, 64, 27.2103343964943718802),973 BOOST_MATH_BIG_CONSTANT(T, 64, -8785.48528692879413676),974 BOOST_MATH_BIG_CONSTANT(T, 64, 397530.290000322626766),975 BOOST_MATH_BIG_CONSTANT(T, 64, -7356441.34957799368252),976 BOOST_MATH_BIG_CONSTANT(T, 64, 63050914.5343400957524),977 BOOST_MATH_BIG_CONSTANT(T, 64, -246143779.638307701369),978 BOOST_MATH_BIG_CONSTANT(T, 64, 384647824.678554961174),979 BOOST_MATH_BIG_CONSTANT(T, 64, -166288297.874583961493)980 };981 // LCOV_EXCL_STOP982 T t = 1 / z;983 result = Y + tools::evaluate_polynomial(P, t)984 / tools::evaluate_polynomial(Q, t);985 if(z < 41)986 result *= exp(z) / z;987 else988 {989 // Avoid premature overflow if we can:990 t = z - 40;991 if(t > tools::log_max_value<T>())992 {993 result = policies::raise_overflow_error<T>(function, nullptr, pol);994 }995 else996 {997 result *= exp(z - 40) / z;998 if(result > tools::max_value<T>() / exp40)999 {1000 result = policies::raise_overflow_error<T>(function, nullptr, pol);1001 }1002 else1003 {1004 result *= exp40;1005 }1006 }1007 }1008 result += z;1009 }1010 return result;1011}1012 1013template <class T, class Policy>1014BOOST_MATH_GPU_ENABLED void expint_i_imp_113a(T& result, const T& z, const Policy& pol)1015{1016 BOOST_MATH_STD_USING1017 // Maximum Deviation Found: 1.230e-361018 // Expected Error Term: -1.230e-361019 // Max Error found at long double precision = Poly: 4.355299e-34 Cheb: 7.512581e-341020 1021 // LCOV_EXCL_START1022 static const T P[15] = {1023 BOOST_MATH_BIG_CONSTANT(T, 113, 2.98677224343598593765287235997328555),1024 BOOST_MATH_BIG_CONSTANT(T, 113, -0.333256034674702967028780537349334037),1025 BOOST_MATH_BIG_CONSTANT(T, 113, 0.851831522798101228384971644036708463),1026 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0657854833494646206186773614110374948),1027 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0630065662557284456000060708977935073),1028 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00311759191425309373327784154659649232),1029 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00176213568201493949664478471656026771),1030 BOOST_MATH_BIG_CONSTANT(T, 113, -0.491548660404172089488535218163952295e-4),1031 BOOST_MATH_BIG_CONSTANT(T, 113, 0.207764227621061706075562107748176592e-4),1032 BOOST_MATH_BIG_CONSTANT(T, 113, -0.225445398156913584846374273379402765e-6),1033 BOOST_MATH_BIG_CONSTANT(T, 113, 0.996939977231410319761273881672601592e-7),1034 BOOST_MATH_BIG_CONSTANT(T, 113, 0.212546902052178643330520878928100847e-9),1035 BOOST_MATH_BIG_CONSTANT(T, 113, 0.154646053060262871360159325115980023e-9),1036 BOOST_MATH_BIG_CONSTANT(T, 113, 0.143971277122049197323415503594302307e-11),1037 BOOST_MATH_BIG_CONSTANT(T, 113, 0.306243138978114692252817805327426657e-13)1038 };1039 static const T Q[15] = {1040 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),1041 BOOST_MATH_BIG_CONSTANT(T, 113, -1.40178870313943798705491944989231793),1042 BOOST_MATH_BIG_CONSTANT(T, 113, 0.943810968269701047641218856758605284),1043 BOOST_MATH_BIG_CONSTANT(T, 113, -0.405026631534345064600850391026113165),1044 BOOST_MATH_BIG_CONSTANT(T, 113, 0.123924153524614086482627660399122762),1045 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0286364505373369439591132549624317707),1046 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00516148845910606985396596845494015963),1047 BOOST_MATH_BIG_CONSTANT(T, 113, -0.000738330799456364820380739850924783649),1048 BOOST_MATH_BIG_CONSTANT(T, 113, 0.843737760991856114061953265870882637e-4),1049 BOOST_MATH_BIG_CONSTANT(T, 113, -0.767957673431982543213661388914587589e-5),1050 BOOST_MATH_BIG_CONSTANT(T, 113, 0.549136847313854595809952100614840031e-6),1051 BOOST_MATH_BIG_CONSTANT(T, 113, -0.299801381513743676764008325949325404e-7),1052 BOOST_MATH_BIG_CONSTANT(T, 113, 0.118419479055346106118129130945423483e-8),1053 BOOST_MATH_BIG_CONSTANT(T, 113, -0.30372295663095470359211949045344607e-10),1054 BOOST_MATH_BIG_CONSTANT(T, 113, 0.382742953753485333207877784720070523e-12)1055 };1056 1057 static const T c1 = BOOST_MATH_BIG_CONSTANT(T, 113, 1677624236387711.0);1058 static const T c2 = BOOST_MATH_BIG_CONSTANT(T, 113, 4503599627370496.0);1059 static const T c3 = BOOST_MATH_BIG_CONSTANT(T, 113, 266514582277687.0);1060 static const T c4 = BOOST_MATH_BIG_CONSTANT(T, 113, 4503599627370496.0);1061 static const T c5 = BOOST_MATH_BIG_CONSTANT(T, 113, 4503599627370496.0);1062 static const T r1 = c1 / c2;1063 static const T r2 = c3 / c4 / c5;1064 static const T r3 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, 0.283806480836357377069325311780969887585024578164571984232357e-31));1065 static const T r = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, 0.372507410781366634461991866580119133535689497771654051555657435242200120636201854384926049951548942392));1066 // LCOV_EXCL_STOP1067 T t = (z / 3) - 1;1068 result = tools::evaluate_polynomial(P, t)1069 / tools::evaluate_polynomial(Q, t);1070 t = ((z - r1) - r2) - r3;1071 result *= t;1072 if(fabs(t) < 0.1)1073 {1074 result += boost::math::log1p(t / r, pol);1075 }1076 else1077 {1078 result += log(z / r);1079 }1080}1081 1082template <class T>1083BOOST_MATH_GPU_ENABLED void expint_i_113b(T& result, const T& z)1084{1085 BOOST_MATH_STD_USING1086 // Maximum Deviation Found: 7.779e-361087 // Expected Error Term: -7.779e-361088 // Max Error found at long double precision = Poly: 2.576723e-35 Cheb: 1.236001e-341089 // LCOV_EXCL_START1090 static const T Y = 1.158985137939453125F;1091 static const T P[15] = {1092 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00139324086199409049282472239613554817),1093 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0338173111691991289178779840307998955),1094 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0555972290794371306259684845277620556),1095 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0378677976003456171563136909186202177),1096 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0152221583517528358782902783914356667),1097 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00428283334203873035104248217403126905),1098 BOOST_MATH_BIG_CONSTANT(T, 113, -0.000922782631491644846511553601323435286),1099 BOOST_MATH_BIG_CONSTANT(T, 113, -0.000155513428088853161562660696055496696),1100 BOOST_MATH_BIG_CONSTANT(T, 113, -0.205756580255359882813545261519317096e-4),1101 BOOST_MATH_BIG_CONSTANT(T, 113, -0.220327406578552089820753181821115181e-5),1102 BOOST_MATH_BIG_CONSTANT(T, 113, -0.189483157545587592043421445645377439e-6),1103 BOOST_MATH_BIG_CONSTANT(T, 113, -0.122426571518570587750898968123803867e-7),1104 BOOST_MATH_BIG_CONSTANT(T, 113, -0.635187358949437991465353268374523944e-9),1105 BOOST_MATH_BIG_CONSTANT(T, 113, -0.203015132965870311935118337194860863e-10),1106 BOOST_MATH_BIG_CONSTANT(T, 113, -0.384276705503357655108096065452950822e-12)1107 };1108 static const T Q[15] = {1109 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),1110 BOOST_MATH_BIG_CONSTANT(T, 113, 1.58784732785354597996617046880946257),1111 BOOST_MATH_BIG_CONSTANT(T, 113, 1.18550755302279446339364262338114098),1112 BOOST_MATH_BIG_CONSTANT(T, 113, 0.55598993549661368604527040349702836),1113 BOOST_MATH_BIG_CONSTANT(T, 113, 0.184290888380564236919107835030984453),1114 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0459658051803613282360464632326866113),1115 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0089505064268613225167835599456014705),1116 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00139042673882987693424772855926289077),1117 BOOST_MATH_BIG_CONSTANT(T, 113, 0.000174210708041584097450805790176479012),1118 BOOST_MATH_BIG_CONSTANT(T, 113, 0.176324034009707558089086875136647376e-4),1119 BOOST_MATH_BIG_CONSTANT(T, 113, 0.142935845999505649273084545313710581e-5),1120 BOOST_MATH_BIG_CONSTANT(T, 113, 0.907502324487057260675816233312747784e-7),1121 BOOST_MATH_BIG_CONSTANT(T, 113, 0.431044337808893270797934621235918418e-8),1122 BOOST_MATH_BIG_CONSTANT(T, 113, 0.139007266881450521776529705677086902e-9),1123 BOOST_MATH_BIG_CONSTANT(T, 113, 0.234715286125516430792452741830364672e-11)1124 };1125 // LCOV_EXCL_STOP1126 T t = z / 2 - 4;1127 result = Y + tools::evaluate_polynomial(P, t)1128 / tools::evaluate_polynomial(Q, t);1129 result *= exp(z) / z;1130 result += z;1131}1132 1133template <class T>1134BOOST_MATH_GPU_ENABLED void expint_i_113c(T& result, const T& z)1135{1136 BOOST_MATH_STD_USING1137 // Maximum Deviation Found: 1.082e-341138 // Expected Error Term: 1.080e-341139 // Max Error found at long double precision = Poly: 1.958294e-34 Cheb: 2.472261e-341140 1141 // LCOV_EXCL_START1142 static const T Y = 1.091579437255859375F;1143 static const T P[17] = {1144 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00685089599550151282724924894258520532),1145 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0443313550253580053324487059748497467),1146 BOOST_MATH_BIG_CONSTANT(T, 113, -0.071538561252424027443296958795814874),1147 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0622923153354102682285444067843300583),1148 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0361631270264607478205393775461208794),1149 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0153192826839624850298106509601033261),1150 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00496967904961260031539602977748408242),1151 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00126989079663425780800919171538920589),1152 BOOST_MATH_BIG_CONSTANT(T, 113, -0.000258933143097125199914724875206326698),1153 BOOST_MATH_BIG_CONSTANT(T, 113, -0.422110326689204794443002330541441956e-4),1154 BOOST_MATH_BIG_CONSTANT(T, 113, -0.546004547590412661451073996127115221e-5),1155 BOOST_MATH_BIG_CONSTANT(T, 113, -0.546775260262202177131068692199272241e-6),1156 BOOST_MATH_BIG_CONSTANT(T, 113, -0.404157632825805803833379568956559215e-7),1157 BOOST_MATH_BIG_CONSTANT(T, 113, -0.200612596196561323832327013027419284e-8),1158 BOOST_MATH_BIG_CONSTANT(T, 113, -0.502538501472133913417609379765434153e-10),1159 BOOST_MATH_BIG_CONSTANT(T, 113, -0.326283053716799774936661568391296584e-13),1160 BOOST_MATH_BIG_CONSTANT(T, 113, 0.869226483473172853557775877908693647e-15)1161 };1162 static const T Q[15] = {1163 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),1164 BOOST_MATH_BIG_CONSTANT(T, 113, 2.23227220874479061894038229141871087),1165 BOOST_MATH_BIG_CONSTANT(T, 113, 2.40221000361027971895657505660959863),1166 BOOST_MATH_BIG_CONSTANT(T, 113, 1.65476320985936174728238416007084214),1167 BOOST_MATH_BIG_CONSTANT(T, 113, 0.816828602963895720369875535001248227),1168 BOOST_MATH_BIG_CONSTANT(T, 113, 0.306337922909446903672123418670921066),1169 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0902400121654409267774593230720600752),1170 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0212708882169429206498765100993228086),1171 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00404442626252467471957713495828165491),1172 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0006195601618842253612635241404054589),1173 BOOST_MATH_BIG_CONSTANT(T, 113, 0.755930932686543009521454653994321843e-4),1174 BOOST_MATH_BIG_CONSTANT(T, 113, 0.716004532773778954193609582677482803e-5),1175 BOOST_MATH_BIG_CONSTANT(T, 113, 0.500881663076471627699290821742924233e-6),1176 BOOST_MATH_BIG_CONSTANT(T, 113, 0.233593219218823384508105943657387644e-7),1177 BOOST_MATH_BIG_CONSTANT(T, 113, 0.554900353169148897444104962034267682e-9)1178 };1179 // LCOV_EXCL_STOP1180 T t = z / 4 - 3.5;1181 result = Y + tools::evaluate_polynomial(P, t)1182 / tools::evaluate_polynomial(Q, t);1183 result *= exp(z) / z;1184 result += z;1185}1186 1187template <class T>1188BOOST_MATH_GPU_ENABLED void expint_i_113d(T& result, const T& z)1189{1190 BOOST_MATH_STD_USING1191 // Maximum Deviation Found: 3.163e-351192 // Expected Error Term: 3.163e-351193 // Max Error found at long double precision = Poly: 4.158110e-35 Cheb: 5.385532e-351194 // LCOV_EXCL_START1195 static const T Y = 1.051731109619140625F;1196 static const T P[14] = {1197 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00144552494420652573815404828020593565),1198 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0126747451594545338365684731262912741),1199 BOOST_MATH_BIG_CONSTANT(T, 113, -0.01757394877502366717526779263438073),1200 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0126838952395506921945756139424722588),1201 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0060045057928894974954756789352443522),1202 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00205349237147226126653803455793107903),1203 BOOST_MATH_BIG_CONSTANT(T, 113, -0.000532606040579654887676082220195624207),1204 BOOST_MATH_BIG_CONSTANT(T, 113, -0.000107344687098019891474772069139014662),1205 BOOST_MATH_BIG_CONSTANT(T, 113, -0.169536802705805811859089949943435152e-4),1206 BOOST_MATH_BIG_CONSTANT(T, 113, -0.20863311729206543881826553010120078e-5),1207 BOOST_MATH_BIG_CONSTANT(T, 113, -0.195670358542116256713560296776654385e-6),1208 BOOST_MATH_BIG_CONSTANT(T, 113, -0.133291168587253145439184028259772437e-7),1209 BOOST_MATH_BIG_CONSTANT(T, 113, -0.595500337089495614285777067722823397e-9),1210 BOOST_MATH_BIG_CONSTANT(T, 113, -0.133141358866324100955927979606981328e-10)1211 };1212 static const T Q[14] = {1213 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),1214 BOOST_MATH_BIG_CONSTANT(T, 113, 1.72490783907582654629537013560044682),1215 BOOST_MATH_BIG_CONSTANT(T, 113, 1.44524329516800613088375685659759765),1216 BOOST_MATH_BIG_CONSTANT(T, 113, 0.778241785539308257585068744978050181),1217 BOOST_MATH_BIG_CONSTANT(T, 113, 0.300520486589206605184097270225725584),1218 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0879346899691339661394537806057953957),1219 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0200802415843802892793583043470125006),1220 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00362842049172586254520256100538273214),1221 BOOST_MATH_BIG_CONSTANT(T, 113, 0.000519731362862955132062751246769469957),1222 BOOST_MATH_BIG_CONSTANT(T, 113, 0.584092147914050999895178697392282665e-4),1223 BOOST_MATH_BIG_CONSTANT(T, 113, 0.501851497707855358002773398333542337e-5),1224 BOOST_MATH_BIG_CONSTANT(T, 113, 0.313085677467921096644895738538865537e-6),1225 BOOST_MATH_BIG_CONSTANT(T, 113, 0.127552010539733113371132321521204458e-7),1226 BOOST_MATH_BIG_CONSTANT(T, 113, 0.25737310826983451144405899970774587e-9)1227 };1228 // LCOV_EXCL_STOP1229 T t = z / 4 - 5.5;1230 result = Y + tools::evaluate_polynomial(P, t)1231 / tools::evaluate_polynomial(Q, t);1232 BOOST_MATH_INSTRUMENT_VARIABLE(result)1233 result *= exp(z) / z;1234 BOOST_MATH_INSTRUMENT_VARIABLE(result)1235 result += z;1236 BOOST_MATH_INSTRUMENT_VARIABLE(result)1237}1238 1239template <class T>1240BOOST_MATH_GPU_ENABLED void expint_i_113e(T& result, const T& z)1241{1242 BOOST_MATH_STD_USING1243 // Maximum Deviation Found: 7.972e-361244 // Expected Error Term: 7.962e-361245 // Max Error found at long double precision = Poly: 1.711721e-34 Cheb: 3.100018e-341246 // LCOV_EXCL_START1247 static const T Y = 1.032726287841796875F;1248 static const T P[15] = {1249 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00141056919297307534690895009969373233),1250 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0123384175302540291339020257071411437),1251 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0298127270706864057791526083667396115),1252 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0390686759471630584626293670260768098),1253 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0338226792912607409822059922949035589),1254 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0211659736179834946452561197559654582),1255 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0100428887460879377373158821400070313),1256 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00370717396015165148484022792801682932),1257 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0010768667551001624764329000496561659),1258 BOOST_MATH_BIG_CONSTANT(T, 113, -0.000246127328761027039347584096573123531),1259 BOOST_MATH_BIG_CONSTANT(T, 113, -0.437318110527818613580613051861991198e-4),1260 BOOST_MATH_BIG_CONSTANT(T, 113, -0.587532682329299591501065482317771497e-5),1261 BOOST_MATH_BIG_CONSTANT(T, 113, -0.565697065670893984610852937110819467e-6),1262 BOOST_MATH_BIG_CONSTANT(T, 113, -0.350233957364028523971768887437839573e-7),1263 BOOST_MATH_BIG_CONSTANT(T, 113, -0.105428907085424234504608142258423505e-8)1264 };1265 static const T Q[16] = {1266 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),1267 BOOST_MATH_BIG_CONSTANT(T, 113, 3.17261315255467581204685605414005525),1268 BOOST_MATH_BIG_CONSTANT(T, 113, 4.85267952971640525245338392887217426),1269 BOOST_MATH_BIG_CONSTANT(T, 113, 4.74341914912439861451492872946725151),1270 BOOST_MATH_BIG_CONSTANT(T, 113, 3.31108463283559911602405970817931801),1271 BOOST_MATH_BIG_CONSTANT(T, 113, 1.74657006336994649386607925179848899),1272 BOOST_MATH_BIG_CONSTANT(T, 113, 0.718255607416072737965933040353653244),1273 BOOST_MATH_BIG_CONSTANT(T, 113, 0.234037553177354542791975767960643864),1274 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0607470145906491602476833515412605389),1275 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0125048143774226921434854172947548724),1276 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00201034366420433762935768458656609163),1277 BOOST_MATH_BIG_CONSTANT(T, 113, 0.000244823338417452367656368849303165721),1278 BOOST_MATH_BIG_CONSTANT(T, 113, 0.213511655166983177960471085462540807e-4),1279 BOOST_MATH_BIG_CONSTANT(T, 113, 0.119323998465870686327170541547982932e-5),1280 BOOST_MATH_BIG_CONSTANT(T, 113, 0.322153582559488797803027773591727565e-7),1281 BOOST_MATH_BIG_CONSTANT(T, 113, -0.161635525318683508633792845159942312e-16)1282 };1283 // LCOV_EXCL_STOP1284 T t = z / 8 - 4.25;1285 result = Y + tools::evaluate_polynomial(P, t)1286 / tools::evaluate_polynomial(Q, t);1287 BOOST_MATH_INSTRUMENT_VARIABLE(result)1288 result *= exp(z) / z;1289 BOOST_MATH_INSTRUMENT_VARIABLE(result)1290 result += z;1291 BOOST_MATH_INSTRUMENT_VARIABLE(result)1292}1293 1294template <class T>1295BOOST_MATH_GPU_ENABLED void expint_i_113f(T& result, const T& z)1296{1297 BOOST_MATH_STD_USING1298 // Maximum Deviation Found: 4.469e-361299 // Expected Error Term: 4.468e-361300 // Max Error found at long double precision = Poly: 1.288958e-35 Cheb: 2.304586e-351301 // LCOV_EXCL_START1302 static const T Y = 1.0216197967529296875F;1303 static const T P[12] = {1304 BOOST_MATH_BIG_CONSTANT(T, 113, -0.000322999116096627043476023926572650045),1305 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00385606067447365187909164609294113346),1306 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00686514524727568176735949971985244415),1307 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00606260649593050194602676772589601799),1308 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00334382362017147544335054575436194357),1309 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00126108534260253075708625583630318043),1310 BOOST_MATH_BIG_CONSTANT(T, 113, -0.000337881489347846058951220431209276776),1311 BOOST_MATH_BIG_CONSTANT(T, 113, -0.648480902304640018785370650254018022e-4),1312 BOOST_MATH_BIG_CONSTANT(T, 113, -0.87652644082970492211455290209092766e-5),1313 BOOST_MATH_BIG_CONSTANT(T, 113, -0.794712243338068631557849449519994144e-6),1314 BOOST_MATH_BIG_CONSTANT(T, 113, -0.434084023639508143975983454830954835e-7),1315 BOOST_MATH_BIG_CONSTANT(T, 113, -0.107839681938752337160494412638656696e-8)1316 };1317 static const T Q[12] = {1318 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),1319 BOOST_MATH_BIG_CONSTANT(T, 113, 2.09913805456661084097134805151524958),1320 BOOST_MATH_BIG_CONSTANT(T, 113, 2.07041755535439919593503171320431849),1321 BOOST_MATH_BIG_CONSTANT(T, 113, 1.26406517226052371320416108604874734),1322 BOOST_MATH_BIG_CONSTANT(T, 113, 0.529689923703770353961553223973435569),1323 BOOST_MATH_BIG_CONSTANT(T, 113, 0.159578150879536711042269658656115746),1324 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0351720877642000691155202082629857131),1325 BOOST_MATH_BIG_CONSTANT(T, 113, 0.00565313621289648752407123620997063122),1326 BOOST_MATH_BIG_CONSTANT(T, 113, 0.000646920278540515480093843570291218295),1327 BOOST_MATH_BIG_CONSTANT(T, 113, 0.499904084850091676776993523323213591e-4),1328 BOOST_MATH_BIG_CONSTANT(T, 113, 0.233740058688179614344680531486267142e-5),1329 BOOST_MATH_BIG_CONSTANT(T, 113, 0.498800627828842754845418576305379469e-7)1330 };1331 // LCOV_EXCL_STOP1332 T t = z / 7 - 7;1333 result = Y + tools::evaluate_polynomial(P, t)1334 / tools::evaluate_polynomial(Q, t);1335 BOOST_MATH_INSTRUMENT_VARIABLE(result)1336 result *= exp(z) / z;1337 BOOST_MATH_INSTRUMENT_VARIABLE(result)1338 result += z;1339 BOOST_MATH_INSTRUMENT_VARIABLE(result)1340}1341 1342template <class T>1343BOOST_MATH_GPU_ENABLED void expint_i_113g(T& result, const T& z)1344{1345 BOOST_MATH_STD_USING1346 // Maximum Deviation Found: 5.588e-351347 // Expected Error Term: -5.566e-351348 // Max Error found at long double precision = Poly: 9.976345e-35 Cheb: 8.358865e-351349 // LCOV_EXCL_START1350 static const T Y = 1.015148162841796875F;1351 static const T P[11] = {1352 BOOST_MATH_BIG_CONSTANT(T, 113, -0.000435714784725086961464589957142615216),1353 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00432114324353830636009453048419094314),1354 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0100740363285526177522819204820582424),1355 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0116744115827059174392383504427640362),1356 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00816145387784261141360062395898644652),1357 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00371380272673500791322744465394211508),1358 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00112958263488611536502153195005736563),1359 BOOST_MATH_BIG_CONSTANT(T, 113, -0.000228316462389404645183269923754256664),1360 BOOST_MATH_BIG_CONSTANT(T, 113, -0.29462181955852860250359064291292577e-4),1361 BOOST_MATH_BIG_CONSTANT(T, 113, -0.21972450610957417963227028788460299e-5),1362 BOOST_MATH_BIG_CONSTANT(T, 113, -0.720558173805289167524715527536874694e-7)1363 };1364 static const T Q[11] = {1365 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),1366 BOOST_MATH_BIG_CONSTANT(T, 113, 2.95918362458402597039366979529287095),1367 BOOST_MATH_BIG_CONSTANT(T, 113, 3.96472247520659077944638411856748924),1368 BOOST_MATH_BIG_CONSTANT(T, 113, 3.15563251550528513747923714884142131),1369 BOOST_MATH_BIG_CONSTANT(T, 113, 1.64674612007093983894215359287448334),1370 BOOST_MATH_BIG_CONSTANT(T, 113, 0.58695020129846594405856226787156424),1371 BOOST_MATH_BIG_CONSTANT(T, 113, 0.144358385319329396231755457772362793),1372 BOOST_MATH_BIG_CONSTANT(T, 113, 0.024146911506411684815134916238348063),1373 BOOST_MATH_BIG_CONSTANT(T, 113, 0.0026257132337460784266874572001650153),1374 BOOST_MATH_BIG_CONSTANT(T, 113, 0.000167479843750859222348869769094711093),1375 BOOST_MATH_BIG_CONSTANT(T, 113, 0.475673638665358075556452220192497036e-5)1376 };1377 // LCOV_EXCL_STOP1378 T t = z / 14 - 5;1379 result = Y + tools::evaluate_polynomial(P, t)1380 / tools::evaluate_polynomial(Q, t);1381 BOOST_MATH_INSTRUMENT_VARIABLE(result)1382 result *= exp(z) / z;1383 BOOST_MATH_INSTRUMENT_VARIABLE(result)1384 result += z;1385 BOOST_MATH_INSTRUMENT_VARIABLE(result)1386}1387 1388template <class T>1389BOOST_MATH_GPU_ENABLED void expint_i_113h(T& result, const T& z)1390{1391 BOOST_MATH_STD_USING1392 // Maximum Deviation Found: 4.448e-361393 // Expected Error Term: 4.445e-361394 // Max Error found at long double precision = Poly: 2.058532e-35 Cheb: 2.165465e-271395 // LCOV_EXCL_START1396 static const T Y= 1.00849151611328125F;1397 static const T P[9] = {1398 BOOST_MATH_BIG_CONSTANT(T, 113, -0.0084915161132812500000001440233607358),1399 BOOST_MATH_BIG_CONSTANT(T, 113, 1.84479378737716028341394223076147872),1400 BOOST_MATH_BIG_CONSTANT(T, 113, -130.431146923726715674081563022115568),1401 BOOST_MATH_BIG_CONSTANT(T, 113, 4336.26945491571504885214176203512015),1402 BOOST_MATH_BIG_CONSTANT(T, 113, -76279.0031974974730095170437591004177),1403 BOOST_MATH_BIG_CONSTANT(T, 113, 729577.956271997673695191455111727774),1404 BOOST_MATH_BIG_CONSTANT(T, 113, -3661928.69330208734947103004900349266),1405 BOOST_MATH_BIG_CONSTANT(T, 113, 8570600.041606912735872059184527855),1406 BOOST_MATH_BIG_CONSTANT(T, 113, -6758379.93672362080947905580906028645)1407 };1408 static const T Q[10] = {1409 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),1410 BOOST_MATH_BIG_CONSTANT(T, 113, -99.4868026047611434569541483506091713),1411 BOOST_MATH_BIG_CONSTANT(T, 113, 3879.67753690517114249705089803055473),1412 BOOST_MATH_BIG_CONSTANT(T, 113, -76495.82413252517165830203774900806),1413 BOOST_MATH_BIG_CONSTANT(T, 113, 820773.726408311894342553758526282667),1414 BOOST_MATH_BIG_CONSTANT(T, 113, -4803087.64956923577571031564909646579),1415 BOOST_MATH_BIG_CONSTANT(T, 113, 14521246.227703545012713173740895477),1416 BOOST_MATH_BIG_CONSTANT(T, 113, -19762752.0196769712258527849159393044),1417 BOOST_MATH_BIG_CONSTANT(T, 113, 8354144.67882768405803322344185185517),1418 BOOST_MATH_BIG_CONSTANT(T, 113, 355076.853106511136734454134915432571)1419 };1420 // LCOV_EXCL_STOP1421 T t = 1 / z;1422 result = Y + tools::evaluate_polynomial(P, t)1423 / tools::evaluate_polynomial(Q, t);1424 result *= exp(z) / z;1425 result += z;1426}1427 1428template <class T, class Policy>1429BOOST_MATH_GPU_ENABLED T expint_i_imp(T z, const Policy& pol, const boost::math::integral_constant<int, 113>& tag)1430{1431 BOOST_MATH_STD_USING1432 constexpr auto function = "boost::math::expint<%1%>(%1%)";1433 if(z < 0)1434 return -expint_imp(1, T(-z), pol, tag);1435 if(z == 0)1436 return -policies::raise_overflow_error<T>(function, nullptr, pol);1437 1438 T result;1439 1440 if(z <= 6)1441 {1442 expint_i_imp_113a(result, z, pol);1443 }1444 else if (z <= 10)1445 {1446 expint_i_113b(result, z);1447 }1448 else if(z <= 18)1449 {1450 expint_i_113c(result, z);1451 }1452 else if(z <= 26)1453 {1454 expint_i_113d(result, z);1455 }1456 else if(z <= 42)1457 {1458 expint_i_113e(result, z);1459 }1460 else if(z <= 56)1461 {1462 expint_i_113f(result, z);1463 }1464 else if(z <= 84)1465 {1466 expint_i_113g(result, z);1467 }1468 else if(z <= 210)1469 {1470 expint_i_113h(result, z);1471 }1472 else // z > 2101473 {1474 // Maximum Deviation Found: 3.963e-371475 // Expected Error Term: 3.963e-371476 // Max Error found at long double precision = Poly: 1.248049e-36 Cheb: 2.843486e-291477 1478 static const T exp40 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, 2.35385266837019985407899910749034804508871617254555467236651e17));1479 static const T Y= 1.00252532958984375F;1480 static const T P[8] = {1481 BOOST_MATH_BIG_CONSTANT(T, 113, -0.00252532958984375000000000000000000085),1482 BOOST_MATH_BIG_CONSTANT(T, 113, 1.16591386866059087390621952073890359),1483 BOOST_MATH_BIG_CONSTANT(T, 113, -67.8483431314018462417456828499277579),1484 BOOST_MATH_BIG_CONSTANT(T, 113, 1567.68688154683822956359536287575892),1485 BOOST_MATH_BIG_CONSTANT(T, 113, -17335.4683325819116482498725687644986),1486 BOOST_MATH_BIG_CONSTANT(T, 113, 93632.6567462673524739954389166550069),1487 BOOST_MATH_BIG_CONSTANT(T, 113, -225025.189335919133214440347510936787),1488 BOOST_MATH_BIG_CONSTANT(T, 113, 175864.614717440010942804684741336853)1489 };1490 static const T Q[9] = {1491 BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),1492 BOOST_MATH_BIG_CONSTANT(T, 113, -65.6998869881600212224652719706425129),1493 BOOST_MATH_BIG_CONSTANT(T, 113, 1642.73850032324014781607859416890077),1494 BOOST_MATH_BIG_CONSTANT(T, 113, -19937.2610222467322481947237312818575),1495 BOOST_MATH_BIG_CONSTANT(T, 113, 124136.267326632742667972126625064538),1496 BOOST_MATH_BIG_CONSTANT(T, 113, -384614.251466704550678760562965502293),1497 BOOST_MATH_BIG_CONSTANT(T, 113, 523355.035910385688578278384032026998),1498 BOOST_MATH_BIG_CONSTANT(T, 113, -217809.552260834025885677791936351294),1499 BOOST_MATH_BIG_CONSTANT(T, 113, -8555.81719551123640677261226549550872)1500 };1501 T t = 1 / z;1502 result = Y + tools::evaluate_polynomial(P, t)1503 / tools::evaluate_polynomial(Q, t);1504 if(z < 41)1505 result *= exp(z) / z;1506 else1507 {1508 // Avoid premature overflow if we can:1509 t = z - 40;1510 if(t > tools::log_max_value<T>())1511 {1512 result = policies::raise_overflow_error<T>(function, nullptr, pol);1513 }1514 else1515 {1516 result *= exp(z - 40) / z;1517 if(result > tools::max_value<T>() / exp40)1518 {1519 result = policies::raise_overflow_error<T>(function, nullptr, pol);1520 }1521 else1522 {1523 result *= exp40;1524 }1525 }1526 }1527 result += z;1528 }1529 return result;1530}1531 1532template <class T, class Policy>1533BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type1534 expint_forwarder(T z, const Policy& /*pol*/, boost::math::true_type const&)1535{1536 typedef typename tools::promote_args<T>::type result_type;1537 typedef typename policies::evaluation<result_type, Policy>::type value_type;1538 typedef typename policies::precision<result_type, Policy>::type precision_type;1539 typedef typename policies::normalise<1540 Policy,1541 policies::promote_float<false>,1542 policies::promote_double<false>,1543 policies::discrete_quantile<>,1544 policies::assert_undefined<> >::type forwarding_policy;1545 typedef boost::math::integral_constant<int,1546 precision_type::value <= 0 ? 0 :1547 precision_type::value <= 53 ? 53 :1548 precision_type::value <= 64 ? 64 :1549 precision_type::value <= 113 ? 113 : 01550 > tag_type;1551 1552 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expint_i_imp(static_cast<value_type>(z), forwarding_policy(), tag_type()), "boost::math::expint<%1%>(%1%)");1553}1554 1555template <class T>1556BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type1557expint_forwarder(unsigned n, T z, const boost::math::false_type&)1558{1559 return boost::math::expint(n, z, policies::policy<>());1560}1561 1562} // namespace detail1563 1564template <class T, class Policy>1565BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type1566 expint(unsigned n, T z, const Policy& /*pol*/)1567{1568 typedef typename tools::promote_args<T>::type result_type;1569 typedef typename policies::evaluation<result_type, Policy>::type value_type;1570 typedef typename policies::precision<result_type, Policy>::type precision_type;1571 typedef typename policies::normalise<1572 Policy,1573 policies::promote_float<false>,1574 policies::promote_double<false>,1575 policies::discrete_quantile<>,1576 policies::assert_undefined<> >::type forwarding_policy;1577 typedef boost::math::integral_constant<int,1578 precision_type::value <= 0 ? 0 :1579 precision_type::value <= 53 ? 53 :1580 precision_type::value <= 64 ? 64 :1581 precision_type::value <= 113 ? 113 : 01582 > tag_type;1583 1584 return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expint_imp(1585 n,1586 static_cast<value_type>(z),1587 forwarding_policy(),1588 tag_type()), "boost::math::expint<%1%>(unsigned, %1%)");1589}1590 1591template <class T, class U>1592BOOST_MATH_GPU_ENABLED inline typename detail::expint_result<T, U>::type1593 expint(T const z, U const u)1594{1595 typedef typename policies::is_policy<U>::type tag_type;1596 return detail::expint_forwarder(z, u, tag_type());1597}1598 1599template <class T>1600BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type1601 expint(T z)1602{1603 return expint(z, policies::policy<>());1604}1605 1606}} // namespaces1607 1608#ifdef _MSC_VER1609#pragma warning(pop)1610#endif1611 1612#endif // BOOST_MATH_EXPINT_HPP1613 1614 1615