math/special_functions/expint.hpp

0001 //  Copyright John Maddock 2007.
0002 //  Copyright Matt Borland 2024.
0003 //  Use, modification and distribution are subject to the
0004 //  Boost Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_MATH_EXPINT_HPP
0008 #define BOOST_MATH_EXPINT_HPP
0009
0010 #ifdef _MSC_VER
0011 #pragma once
0012 #pragma warning(push)
0013 #pragma warning(disable:4702) // Unreachable code (release mode only warning)
0014 #endif
0015
0016 #include <boost/math/tools/config.hpp>
0017 #include <boost/math/tools/cstdint.hpp>
0018 #include <boost/math/tools/type_traits.hpp>
0019 #include <boost/math/tools/tuple.hpp>
0020 #include <boost/math/tools/precision.hpp>
0021 #include <boost/math/tools/promotion.hpp>
0022 #include <boost/math/tools/fraction.hpp>
0023 #include <boost/math/tools/series.hpp>
0024 #include <boost/math/policies/error_handling.hpp>
0025 #include <boost/math/special_functions/math_fwd.hpp>
0026 #include <boost/math/special_functions/digamma.hpp>
0027 #include <boost/math/special_functions/log1p.hpp>
0028
0029 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
0030 //
0031 // This is the only way we can avoid
0032 // warning: non-standard suffix on floating constant [-Wpedantic]
0033 // when building with -Wall -pedantic.  Neither __extension__
0034 // nor #pragma diagnostic ignored work :(
0035 //
0036 #pragma GCC system_header
0037 #endif
0038
0039 namespace boost{ namespace math{
0040
0041 template <class T, class Policy>
0042 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type
0043    expint(unsigned n, T z, const Policy& /*pol*/);
0044
0045 namespace detail{
0046
0047 template <class T>
0048 BOOST_MATH_GPU_ENABLED inline T expint_1_rational(const T& z, const boost::math::integral_constant<int, 0>&)
0049 {
0050    // this function is never actually called
0051    BOOST_MATH_ASSERT(0);
0052    return z;
0053 }
0054
0055 template <class T>
0056 BOOST_MATH_GPU_ENABLED T expint_1_rational(const T& z, const boost::math::integral_constant<int, 53>&)
0057 {
0058    BOOST_MATH_STD_USING
0059    T result;
0060    if(z <= 1)
0061    {
0062       // Maximum Deviation Found:                     2.006e-18
0063       // Expected Error Term:                         2.006e-18
0064       // Max error found at double precision:         2.760e-17
0065       static const T Y = 0.66373538970947265625F;
0066       static const T P[6] = {
0067          BOOST_MATH_BIG_CONSTANT(T, 53, 0.0865197248079397976498),
0068          BOOST_MATH_BIG_CONSTANT(T, 53, 0.0320913665303559189999),
0069          BOOST_MATH_BIG_CONSTANT(T, 53, -0.245088216639761496153),
0070          BOOST_MATH_BIG_CONSTANT(T, 53, -0.0368031736257943745142),
0071          BOOST_MATH_BIG_CONSTANT(T, 53, -0.00399167106081113256961),
0072          BOOST_MATH_BIG_CONSTANT(T, 53, -0.000111507792921197858394)
0073       };
0074       static const T Q[6] = {
0075          BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),
0076          BOOST_MATH_BIG_CONSTANT(T, 53, 0.37091387659397013215),
0077          BOOST_MATH_BIG_CONSTANT(T, 53, 0.056770677104207528384),
0078          BOOST_MATH_BIG_CONSTANT(T, 53, 0.00427347600017103698101),
0079          BOOST_MATH_BIG_CONSTANT(T, 53, 0.000131049900798434683324),
0080          BOOST_MATH_BIG_CONSTANT(T, 53, -0.528611029520217142048e-6)
0081       };
0082       result = tools::evaluate_polynomial(P, z)
0083          / tools::evaluate_polynomial(Q, z);
0084       result += z - log(z) - Y;
0085    }
0086    else if(z < -boost::math::tools::log_min_value<T>())
0087    {
0088       // Maximum Deviation Found (interpolated):      1.444e-17
0089       // Max error found at double precision:         3.119e-17
0090       static const T P[11] = {
0091          BOOST_MATH_BIG_CONSTANT(T, 53, -0.121013190657725568138e-18),
0092          BOOST_MATH_BIG_CONSTANT(T, 53, -0.999999999999998811143),
0093          BOOST_MATH_BIG_CONSTANT(T, 53, -43.3058660811817946037),
0094          BOOST_MATH_BIG_CONSTANT(T, 53, -724.581482791462469795),
0095          BOOST_MATH_BIG_CONSTANT(T, 53, -6046.8250112711035463),
0096          BOOST_MATH_BIG_CONSTANT(T, 53, -27182.6254466733970467),
0097          BOOST_MATH_BIG_CONSTANT(T, 53, -66598.2652345418633509),
0098          BOOST_MATH_BIG_CONSTANT(T, 53, -86273.1567711649528784),
0099          BOOST_MATH_BIG_CONSTANT(T, 53, -54844.4587226402067411),
0100          BOOST_MATH_BIG_CONSTANT(T, 53, -14751.4895786128450662),
0101          BOOST_MATH_BIG_CONSTANT(T, 53, -1185.45720315201027667)
0102       };
0103       static const T Q[12] = {
0104          BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),
0105          BOOST_MATH_BIG_CONSTANT(T, 53, 45.3058660811801465927),
0106          BOOST_MATH_BIG_CONSTANT(T, 53, 809.193214954550328455),
0107          BOOST_MATH_BIG_CONSTANT(T, 53, 7417.37624454689546708),
0108          BOOST_MATH_BIG_CONSTANT(T, 53, 38129.5594484818471461),
0109          BOOST_MATH_BIG_CONSTANT(T, 53, 113057.05869159631492),
0110          BOOST_MATH_BIG_CONSTANT(T, 53, 192104.047790227984431),
0111          BOOST_MATH_BIG_CONSTANT(T, 53, 180329.498380501819718),
0112          BOOST_MATH_BIG_CONSTANT(T, 53, 86722.3403467334749201),
0113          BOOST_MATH_BIG_CONSTANT(T, 53, 18455.4124737722049515),
0114          BOOST_MATH_BIG_CONSTANT(T, 53, 1229.20784182403048905),
0115          BOOST_MATH_BIG_CONSTANT(T, 53, -0.776491285282330997549)
0116       };
0117       T recip = 1 / z;
0118       result = 1 + tools::evaluate_polynomial(P, recip)
0119          / tools::evaluate_polynomial(Q, recip);
0120       result *= exp(-z) * recip;
0121    }
0122    else
0123    {
0124       result = 0;
0125    }
0126    return result;
0127 }
0128
0129 template <class T>
0130 BOOST_MATH_GPU_ENABLED T expint_1_rational(const T& z, const boost::math::integral_constant<int, 64>&)
0131 {
0132    BOOST_MATH_STD_USING
0133    T result;
0134    if(z <= 1)
0135    {
0136       // Maximum Deviation Found:                     3.807e-20
0137       // Expected Error Term:                         3.807e-20
0138       // Max error found at long double precision:    6.249e-20
0139
0140       static const T Y = 0.66373538970947265625F;
0141       static const T P[6] = {
0142          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0865197248079397956816),
0143          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0275114007037026844633),
0144          BOOST_MATH_BIG_CONSTANT(T, 64, -0.246594388074877139824),
0145          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0237624819878732642231),
0146          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00259113319641673986276),
0147          BOOST_MATH_BIG_CONSTANT(T, 64, 0.30853660894346057053e-4)
0148       };
0149       static const T Q[7] = {
0150          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0151          BOOST_MATH_BIG_CONSTANT(T, 64, 0.317978365797784100273),
0152          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0393622602554758722511),
0153          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00204062029115966323229),
0154          BOOST_MATH_BIG_CONSTANT(T, 64, 0.732512107100088047854e-5),
0155          BOOST_MATH_BIG_CONSTANT(T, 64, -0.202872781770207871975e-5),
0156          BOOST_MATH_BIG_CONSTANT(T, 64, 0.52779248094603709945e-7)
0157       };
0158       result = tools::evaluate_polynomial(P, z)
0159          / tools::evaluate_polynomial(Q, z);
0160       result += z - log(z) - Y;
0161    }
0162    else if(z < -boost::math::tools::log_min_value<T>())
0163    {
0164       // Maximum Deviation Found (interpolated):     2.220e-20
0165       // Max error found at long double precision:   1.346e-19
0166       static const T P[14] = {
0167          BOOST_MATH_BIG_CONSTANT(T, 64, -0.534401189080684443046e-23),
0168          BOOST_MATH_BIG_CONSTANT(T, 64, -0.999999999999999999905),
0169          BOOST_MATH_BIG_CONSTANT(T, 64, -62.1517806091379402505),
0170          BOOST_MATH_BIG_CONSTANT(T, 64, -1568.45688271895145277),
0171          BOOST_MATH_BIG_CONSTANT(T, 64, -21015.3431990874009619),
0172          BOOST_MATH_BIG_CONSTANT(T, 64, -164333.011755931661949),
0173          BOOST_MATH_BIG_CONSTANT(T, 64, -777917.270775426696103),
0174          BOOST_MATH_BIG_CONSTANT(T, 64, -2244188.56195255112937),
0175          BOOST_MATH_BIG_CONSTANT(T, 64, -3888702.98145335643429),
0176          BOOST_MATH_BIG_CONSTANT(T, 64, -3909822.65621952648353),
0177          BOOST_MATH_BIG_CONSTANT(T, 64, -2149033.9538897398457),
0178          BOOST_MATH_BIG_CONSTANT(T, 64, -584705.537139793925189),
0179          BOOST_MATH_BIG_CONSTANT(T, 64, -65815.2605361889477244),
0180          BOOST_MATH_BIG_CONSTANT(T, 64, -2038.82870680427258038)
0181       };
0182       static const T Q[14] = {
0183          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0184          BOOST_MATH_BIG_CONSTANT(T, 64, 64.1517806091379399478),
0185          BOOST_MATH_BIG_CONSTANT(T, 64, 1690.76044393722763785),
0186          BOOST_MATH_BIG_CONSTANT(T, 64, 24035.9534033068949426),
0187          BOOST_MATH_BIG_CONSTANT(T, 64, 203679.998633572361706),
0188          BOOST_MATH_BIG_CONSTANT(T, 64, 1074661.58459976978285),
0189          BOOST_MATH_BIG_CONSTANT(T, 64, 3586552.65020899358773),
0190          BOOST_MATH_BIG_CONSTANT(T, 64, 7552186.84989547621411),
0191          BOOST_MATH_BIG_CONSTANT(T, 64, 9853333.79353054111434),
0192          BOOST_MATH_BIG_CONSTANT(T, 64, 7689642.74550683631258),
0193          BOOST_MATH_BIG_CONSTANT(T, 64, 3385553.35146759180739),
0194          BOOST_MATH_BIG_CONSTANT(T, 64, 763218.072732396428725),
0195          BOOST_MATH_BIG_CONSTANT(T, 64, 73930.2995984054930821),
0196          BOOST_MATH_BIG_CONSTANT(T, 64, 2063.86994219629165937)
0197       };
0198       T recip = 1 / z;
0199       result = 1 + tools::evaluate_polynomial(P, recip)
0200          / tools::evaluate_polynomial(Q, recip);
0201       result *= exp(-z) * recip;
0202    }
0203    else
0204    {
0205       result = 0;
0206    }
0207    return result;
0208 }
0209
0210 template <class T>
0211 BOOST_MATH_GPU_ENABLED T expint_1_rational(const T& z, const boost::math::integral_constant<int, 113>&)
0212 {
0213    BOOST_MATH_STD_USING
0214    T result;
0215    if(z <= 1)
0216    {
0217       // Maximum Deviation Found:                     2.477e-35
0218       // Expected Error Term:                         2.477e-35
0219       // Max error found at long double precision:    6.810e-35
0220
0221       static const T Y = 0.66373538970947265625F;
0222       static const T P[10] = {
0223          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0865197248079397956434879099175975937),
0224          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0369066175910795772830865304506087759),
0225          BOOST_MATH_BIG_CONSTANT(T, 113, -0.24272036838415474665971599314725545),
0226          BOOST_MATH_BIG_CONSTANT(T, 113, -0.0502166331248948515282379137550178307),
0227          BOOST_MATH_BIG_CONSTANT(T, 113, -0.00768384138547489410285101483730424919),
0228          BOOST_MATH_BIG_CONSTANT(T, 113, -0.000612574337702109683505224915484717162),
0229          BOOST_MATH_BIG_CONSTANT(T, 113, -0.380207107950635046971492617061708534e-4),
0230          BOOST_MATH_BIG_CONSTANT(T, 113, -0.136528159460768830763009294683628406e-5),
0231          BOOST_MATH_BIG_CONSTANT(T, 113, -0.346839106212658259681029388908658618e-7),
0232          BOOST_MATH_BIG_CONSTANT(T, 113, -0.340500302777838063940402160594523429e-9)
0233       };
0234       static const T Q[10] = {
0235          BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0236          BOOST_MATH_BIG_CONSTANT(T, 113, 0.426568827778942588160423015589537302),
0237          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0841384046470893490592450881447510148),
0238          BOOST_MATH_BIG_CONSTANT(T, 113, 0.0100557215850668029618957359471132995),
0239          BOOST_MATH_BIG_CONSTANT(T, 113, 0.000799334870474627021737357294799839363),
0240          BOOST_MATH_BIG_CONSTANT(T, 113, 0.434452090903862735242423068552687688e-4),
0241          BOOST_MATH_BIG_CONSTANT(T, 113, 0.15829674748799079874182885081231252e-5),
0242          BOOST_MATH_BIG_CONSTANT(T, 113, 0.354406206738023762100882270033082198e-7),
0243          BOOST_MATH_BIG_CONSTANT(T, 113, 0.369373328141051577845488477377890236e-9),
0244          BOOST_MATH_BIG_CONSTANT(T, 113, -0.274149801370933606409282434677600112e-12)
0245       };
0246       result = tools::evaluate_polynomial(P, z)
0247          / tools::evaluate_polynomial(Q, z);
0248       result += z - log(z) - Y;
0249    }
0250    else if(z <= 4)
0251    {
0252       // Max error in interpolated form:             5.614e-35
0253       // Max error found at long double precision:   7.979e-35
0254
0255       static const T Y = 0.70190334320068359375F;
0256
0257       static const T P[16] = {
0258          BOOST_MATH_BIG_CONSTANT(T, 113, 0.298096656795020369955077350585959794),
0259          BOOST_MATH_BIG_CONSTANT(T, 113, 12.9314045995266142913135497455971247),
0260          BOOST_MATH_BIG_CONSTANT(T, 113, 226.144334921582637462526628217345501),
0261          BOOST_MATH_BIG_CONSTANT(T, 113, 2070.83670924261732722117682067381405),
0262          BOOST_MATH_BIG_CONSTANT(T, 113, 10715.1115684330959908244769731347186),
0263          BOOST_MATH_BIG_CONSTANT(T, 113, 30728.7876355542048019664777316053311),
0264          BOOST_MATH_BIG_CONSTANT(T, 113, 38520.6078609349855436936232610875297),
0265          BOOST_MATH_BIG_CONSTANT(T, 113, -27606.0780981527583168728339620565165),
0266          BOOST_MATH_BIG_CONSTANT(T, 113, -169026.485055785605958655247592604835),
0267          BOOST_MATH_BIG_CONSTANT(T, 113, -254361.919204983608659069868035092282),
0268          BOOST_MATH_BIG_CONSTANT(T, 113, -195765.706874132267953259272028679935),
0269          BOOST_MATH_BIG_CONSTANT(T, 113, -83352.6826013533205474990119962408675),
0270          BOOST_MATH_BIG_CONSTANT(T, 113, -19251.6828496869586415162597993050194),
0271          BOOST_MATH_BIG_CONSTANT(T, 113, -2226.64251774578542836725386936102339),
0272          BOOST_MATH_BIG_CONSTANT(T, 113, -109.009437301400845902228611986479816),
0273          BOOST_MATH_BIG_CONSTANT(T, 113, -1.51492042209561411434644938098833499)
0274       };
0275       static const T Q[16] = {
0276          BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0277          BOOST_MATH_BIG_CONSTANT(T, 113, 46.734521442032505570517810766704587),
0278          BOOST_MATH_BIG_CONSTANT(T, 113, 908.694714348462269000247450058595655),
0279          BOOST_MATH_BIG_CONSTANT(T, 113, 9701.76053033673927362784882748513195),
0280          BOOST_MATH_BIG_CONSTANT(T, 113, 63254.2815292641314236625196594947774),
0281          BOOST_MATH_BIG_CONSTANT(T, 113, 265115.641285880437335106541757711092),
0282          BOOST_MATH_BIG_CONSTANT(T, 113, 732707.841188071900498536533086567735),
0283          BOOST_MATH_BIG_CONSTANT(T, 113, 1348514.02492635723327306628712057794),
0284          BOOST_MATH_BIG_CONSTANT(T, 113, 1649986.81455283047769673308781585991),
0285          BOOST_MATH_BIG_CONSTANT(T, 113, 1326000.828522976970116271208812099),
0286          BOOST_MATH_BIG_CONSTANT(T, 113, 683643.09490612171772350481773951341),
0287          BOOST_MATH_BIG_CONSTANT(T, 113, 217640.505137263607952365685653352229),
0288          BOOST_MATH_BIG_CONSTANT(T, 113, 40288.3467237411710881822569476155485),
0289          BOOST_MATH_BIG_CONSTANT(T, 113, 3932.89353979531632559232883283175754),
0290          BOOST_MATH_BIG_CONSTANT(T, 113, 169.845369689596739824177412096477219),
0291          BOOST_MATH_BIG_CONSTANT(T, 113, 2.17607292280092201170768401876895354)
0292       };
0293       T recip = 1 / z;
0294       result = Y + tools::evaluate_polynomial(P, recip)
0295          / tools::evaluate_polynomial(Q, recip);
0296       result *= exp(-z) * recip;
0297    }
0298    else if(z < -boost::math::tools::log_min_value<T>())
0299    {
0300       // Max error in interpolated form:             4.413e-35
0301       // Max error found at long double precision:   8.928e-35
0302
0303       static const T P[19] = {
0304          BOOST_MATH_BIG_CONSTANT(T, 113, -0.559148411832951463689610809550083986e-40),
0305          BOOST_MATH_BIG_CONSTANT(T, 113, -0.999999999999999999999999999999999997),
0306          BOOST_MATH_BIG_CONSTANT(T, 113, -166.542326331163836642960118190147367),
0307          BOOST_MATH_BIG_CONSTANT(T, 113, -12204.639128796330005065904675153652),
0308          BOOST_MATH_BIG_CONSTANT(T, 113, -520807.069767086071806275022036146855),
0309          BOOST_MATH_BIG_CONSTANT(T, 113, -14435981.5242137970691490903863125326),
0310          BOOST_MATH_BIG_CONSTANT(T, 113, -274574945.737064301247496460758654196),
0311          BOOST_MATH_BIG_CONSTANT(T, 113, -3691611582.99810039356254671781473079),
0312          BOOST_MATH_BIG_CONSTANT(T, 113, -35622515944.8255047299363690814678763),
0313          BOOST_MATH_BIG_CONSTANT(T, 113, -248040014774.502043161750715548451142),
0314          BOOST_MATH_BIG_CONSTANT(T, 113, -1243190389769.53458416330946622607913),
0315          BOOST_MATH_BIG_CONSTANT(T, 113, -4441730126135.54739052731990368425339),
0316          BOOST_MATH_BIG_CONSTANT(T, 113, -11117043181899.7388524310281751971366),
0317          BOOST_MATH_BIG_CONSTANT(T, 113, -18976497615396.9717776601813519498961),
0318          BOOST_MATH_BIG_CONSTANT(T, 113, -21237496819711.1011661104761906067131),
0319          BOOST_MATH_BIG_CONSTANT(T, 113, -14695899122092.5161620333466757812848),
0320          BOOST_MATH_BIG_CONSTANT(T, 113, -5737221535080.30569711574295785864903),
0321          BOOST_MATH_BIG_CONSTANT(T, 113, -1077042281708.42654526404581272546244),
0322          BOOST_MATH_BIG_CONSTANT(T, 113, -68028222642.1941480871395695677675137)
0323       };
0324       static const T Q[20] = {
0325          BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
0326          BOOST_MATH_BIG_CONSTANT(T, 113, 168.542326331163836642960118190147311),
0327          BOOST_MATH_BIG_CONSTANT(T, 113, 12535.7237814586576783518249115343619),
0328          BOOST_MATH_BIG_CONSTANT(T, 113, 544891.263372016404143120911148640627),
0329          BOOST_MATH_BIG_CONSTANT(T, 113, 15454474.7241010258634446523045237762),
0330          BOOST_MATH_BIG_CONSTANT(T, 113, 302495899.896629522673410325891717381),
0331          BOOST_MATH_BIG_CONSTANT(T, 113, 4215565948.38886507646911672693270307),
0332          BOOST_MATH_BIG_CONSTANT(T, 113, 42552409471.7951815668506556705733344),
0333          BOOST_MATH_BIG_CONSTANT(T, 113, 313592377066.753173979584098301610186),
0334          BOOST_MATH_BIG_CONSTANT(T, 113, 1688763640223.4541980740597514904542),
0335          BOOST_MATH_BIG_CONSTANT(T, 113, 6610992294901.59589748057620192145704),
0336          BOOST_MATH_BIG_CONSTANT(T, 113, 18601637235659.6059890851321772682606),
0337          BOOST_MATH_BIG_CONSTANT(T, 113, 36944278231087.2571020964163402941583),
0338          BOOST_MATH_BIG_CONSTANT(T, 113, 50425858518481.7497071917028793820058),
0339          BOOST_MATH_BIG_CONSTANT(T, 113, 45508060902865.0899967797848815980644),
0340          BOOST_MATH_BIG_CONSTANT(T, 113, 25649955002765.3817331501988304758142),
0341          BOOST_MATH_BIG_CONSTANT(T, 113, 8259575619094.6518520988612711292331),
0342          BOOST_MATH_BIG_CONSTANT(T, 113, 1299981487496.12607474362723586264515),
0343          BOOST_MATH_BIG_CONSTANT(T, 113, 70242279152.8241187845178443118302693),
0344          BOOST_MATH_BIG_CONSTANT(T, 113, -37633302.9409263839042721539363416685)
0345       };
0346       T recip = 1 / z;
0347       result = 1 + tools::evaluate_polynomial(P, recip)
0348          / tools::evaluate_polynomial(Q, recip);
0349       result *= exp(-z) * recip;
0350    }
0351    else
0352    {
0353       result = 0;
0354    }
0355    return result;
0356 }
0357
0358
0359 template <class T>
0360 struct expint_fraction
0361 {
0362    typedef boost::math::pair<T,T> result_type;
0363    BOOST_MATH_GPU_ENABLED expint_fraction(unsigned n_, T z_) : b(n_ + z_), i(-1), n(n_){}
0364    BOOST_MATH_GPU_ENABLED boost::math::pair<T,T> operator()()
0365    {
0366       boost::math::pair<T,T> result = boost::math::make_pair(-static_cast<T>((i+1) * (n+i)), b);
0367       b += 2;
0368       ++i;
0369       return result;
0370    }
0371 private:
0372    T b;
0373    int i;
0374    unsigned n;
0375 };
0376
0377 template <class T, class Policy>
0378 BOOST_MATH_GPU_ENABLED inline T expint_as_fraction(unsigned n, T z, const Policy& pol)
0379 {
0380    BOOST_MATH_STD_USING
0381    BOOST_MATH_INSTRUMENT_VARIABLE(z)
0382    boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0383    expint_fraction<T> f(n, z);
0384    T result = tools::continued_fraction_b(
0385       f,
0386       boost::math::policies::get_epsilon<T, Policy>(),
0387       max_iter);
0388    policies::check_series_iterations<T>("boost::math::expint_continued_fraction<%1%>(unsigned,%1%)", max_iter, pol);
0389    BOOST_MATH_INSTRUMENT_VARIABLE(result)
0390    BOOST_MATH_INSTRUMENT_VARIABLE(max_iter)
0391    result = exp(-z) / result;
0392    BOOST_MATH_INSTRUMENT_VARIABLE(result)
0393    return result;
0394 }
0395
0396 template <class T>
0397 struct expint_series
0398 {
0399    typedef T result_type;
0400    BOOST_MATH_GPU_ENABLED expint_series(unsigned k_, T z_, T x_k_, T denom_, T fact_)
0401       : k(k_), z(z_), x_k(x_k_), denom(denom_), fact(fact_){}
0402    BOOST_MATH_GPU_ENABLED T operator()()
0403    {
0404       x_k *= -z;
0405       denom += 1;
0406       fact *= ++k;
0407       return x_k / (denom * fact);
0408    }
0409 private:
0410    unsigned k;
0411    T z;
0412    T x_k;
0413    T denom;
0414    T fact;
0415 };
0416
0417 template <class T, class Policy>
0418 BOOST_MATH_GPU_ENABLED inline T expint_as_series(unsigned n, T z, const Policy& pol)
0419 {
0420    BOOST_MATH_STD_USING
0421    boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0422
0423    BOOST_MATH_INSTRUMENT_VARIABLE(z)
0424
0425    T result = 0;
0426    T x_k = -1;
0427    T denom = T(1) - n;
0428    T fact = 1;
0429    unsigned k = 0;
0430    for(; k < n - 1;)
0431    {
0432       result += x_k / (denom * fact);
0433       denom += 1;
0434       x_k *= -z;
0435       fact *= ++k;
0436    }
0437    BOOST_MATH_INSTRUMENT_VARIABLE(result)
0438    result += pow(-z, static_cast<T>(n - 1))
0439       * (boost::math::digamma(static_cast<T>(n), pol) - log(z)) / fact;
0440    BOOST_MATH_INSTRUMENT_VARIABLE(result)
0441
0442    expint_series<T> s(k, z, x_k, denom, fact);
0443    result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, result);
0444    policies::check_series_iterations<T>("boost::math::expint_series<%1%>(unsigned,%1%)", max_iter, pol);
0445    BOOST_MATH_INSTRUMENT_VARIABLE(result)
0446    BOOST_MATH_INSTRUMENT_VARIABLE(max_iter)
0447    return result;
0448 }
0449
0450 template <class T, class Policy, class Tag>
0451 BOOST_MATH_GPU_ENABLED T expint_imp(unsigned n, T z, const Policy& pol, const Tag& tag)
0452 {
0453    BOOST_MATH_STD_USING
0454    constexpr auto function = "boost::math::expint<%1%>(unsigned, %1%)";
0455    if(z < 0)
0456       return policies::raise_domain_error<T>(function, "Function requires z >= 0 but got %1%.", z, pol);
0457    if(z == 0)
0458       return n == 1 ? policies::raise_overflow_error<T>(function, nullptr, pol) : T(1 / (static_cast<T>(n - 1)));
0459
0460    T result;
0461
0462    bool f;
0463    if(n < 3)
0464    {
0465       f = z < T(0.5);
0466    }
0467    else
0468    {
0469       f = z < (static_cast<T>(n - 2) / static_cast<T>(n - 1));
0470    }
0471 #ifdef _MSC_VER
0472 #  pragma warning(push)
0473 #  pragma warning(disable:4127) // conditional expression is constant
0474 #endif
0475    if(n == 0)
0476    {
0477       result = exp(-z) / z;
0478    }
0479    else if((n == 1) && (Tag::value))
0480    {
0481       result = expint_1_rational(z, tag);
0482    }
0483    else if(f)
0484    {
0485       result = expint_as_series(n, z, pol);
0486    }
0487    else
0488    {
0489       result = expint_as_fraction(n, z, pol);
0490    }
0491 #ifdef _MSC_VER
0492 #  pragma warning(pop)
0493 #endif
0494
0495    return result;
0496 }
0497
0498 template <class T>
0499 struct expint_i_series
0500 {
0501    typedef T result_type;
0502    BOOST_MATH_GPU_ENABLED expint_i_series(T z_) : k(0), z_k(1), z(z_){}
0503    BOOST_MATH_GPU_ENABLED T operator()()
0504    {
0505       z_k *= z / ++k;
0506       return z_k / k;
0507    }
0508 private:
0509    unsigned k;
0510    T z_k;
0511    T z;
0512 };
0513
0514 template <class T, class Policy>
0515 BOOST_MATH_GPU_ENABLED T expint_i_as_series(T z, const Policy& pol)
0516 {
0517    BOOST_MATH_STD_USING
0518    T result = log(z); // (log(z) - log(1 / z)) / 2;
0519    result += constants::euler<T>();
0520    expint_i_series<T> s(z);
0521    boost::math::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
0522    result = tools::sum_series(s, policies::get_epsilon<T, Policy>(), max_iter, result);
0523    policies::check_series_iterations<T>("boost::math::expint_i_series<%1%>(%1%)", max_iter, pol);
0524    return result;
0525 }
0526
0527 template <class T, class Policy, class Tag>
0528 BOOST_MATH_GPU_ENABLED T expint_i_imp(T z, const Policy& pol, const Tag& tag)
0529 {
0530    constexpr auto function = "boost::math::expint<%1%>(%1%)";
0531    if(z < 0)
0532       return -expint_imp(1, T(-z), pol, tag);
0533    if(z == 0)
0534       return -policies::raise_overflow_error<T>(function, nullptr, pol);
0535    return expint_i_as_series(z, pol);
0536 }
0537
0538 template <class T, class Policy>
0539 BOOST_MATH_GPU_ENABLED T expint_i_imp(T z, const Policy& pol, const boost::math::integral_constant<int, 53>& tag)
0540 {
0541    BOOST_MATH_STD_USING
0542    constexpr auto function = "boost::math::expint<%1%>(%1%)";
0543    if(z < 0)
0544       return -expint_imp(1, T(-z), pol, tag);
0545    if(z == 0)
0546       return -policies::raise_overflow_error<T>(function, nullptr, pol);
0547
0548    T result;
0549
0550    if(z <= 6)
0551    {
0552       // Maximum Deviation Found:                     2.852e-18
0553       // Expected Error Term:                         2.852e-18
0554       // Max Error found at double precision =        Poly: 2.636335e-16   Cheb: 4.187027e-16
0555       BOOST_MATH_STATIC const T P[10] = {
0556          BOOST_MATH_BIG_CONSTANT(T, 53, 2.98677224343598593013),
0557          BOOST_MATH_BIG_CONSTANT(T, 53, 0.356343618769377415068),
0558          BOOST_MATH_BIG_CONSTANT(T, 53, 0.780836076283730801839),
0559          BOOST_MATH_BIG_CONSTANT(T, 53, 0.114670926327032002811),
0560          BOOST_MATH_BIG_CONSTANT(T, 53, 0.0499434773576515260534),
0561          BOOST_MATH_BIG_CONSTANT(T, 53, 0.00726224593341228159561),
0562          BOOST_MATH_BIG_CONSTANT(T, 53, 0.00115478237227804306827),
0563          BOOST_MATH_BIG_CONSTANT(T, 53, 0.000116419523609765200999),
0564          BOOST_MATH_BIG_CONSTANT(T, 53, 0.798296365679269702435e-5),
0565          BOOST_MATH_BIG_CONSTANT(T, 53, 0.2777056254402008721e-6)
0566       };
0567       BOOST_MATH_STATIC const T Q[8] = {
0568          BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),
0569          BOOST_MATH_BIG_CONSTANT(T, 53, -1.17090412365413911947),
0570          BOOST_MATH_BIG_CONSTANT(T, 53, 0.62215109846016746276),
0571          BOOST_MATH_BIG_CONSTANT(T, 53, -0.195114782069495403315),
0572          BOOST_MATH_BIG_CONSTANT(T, 53, 0.0391523431392967238166),
0573          BOOST_MATH_BIG_CONSTANT(T, 53, -0.00504800158663705747345),
0574          BOOST_MATH_BIG_CONSTANT(T, 53, 0.000389034007436065401822),
0575          BOOST_MATH_BIG_CONSTANT(T, 53, -0.138972589601781706598e-4)
0576       };
0577
0578       BOOST_MATH_STATIC_LOCAL_VARIABLE const T c1 = BOOST_MATH_BIG_CONSTANT(T, 53, 1677624236387711.0);
0579       BOOST_MATH_STATIC_LOCAL_VARIABLE const T c2 = BOOST_MATH_BIG_CONSTANT(T, 53, 4503599627370496.0);
0580       BOOST_MATH_STATIC_LOCAL_VARIABLE const T r1 = static_cast<T>(c1 / c2);
0581       BOOST_MATH_STATIC_LOCAL_VARIABLE const T r2 = BOOST_MATH_BIG_CONSTANT(T, 53, 0.131401834143860282009280387409357165515556574352422001206362e-16);
0582       BOOST_MATH_STATIC_LOCAL_VARIABLE const T r = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 53, 0.372507410781366634461991866580119133535689497771654051555657435242200120636201854384926049951548942392));
0583       T t = (z / 3) - 1;
0584       result = tools::evaluate_polynomial(P, t)
0585          / tools::evaluate_polynomial(Q, t);
0586       t = (z - r1) - r2;
0587       result *= t;
0588       if(fabs(t) < T(0.1))
0589       {
0590          result += boost::math::log1p(t / r, pol);
0591       }
0592       else
0593       {
0594          result += log(z / r);
0595       }
0596    }
0597    else if (z <= 10)
0598    {
0599       // Maximum Deviation Found:                     6.546e-17
0600       // Expected Error Term:                         6.546e-17
0601       // Max Error found at double precision =        Poly: 6.890169e-17   Cheb: 6.772128e-17
0602       BOOST_MATH_STATIC_LOCAL_VARIABLE const T Y = 1.158985137939453125F;
0603       BOOST_MATH_STATIC const T P[8] = {
0604          BOOST_MATH_BIG_CONSTANT(T, 53, 0.00139324086199402804173),
0605          BOOST_MATH_BIG_CONSTANT(T, 53, -0.0349921221823888744966),
0606          BOOST_MATH_BIG_CONSTANT(T, 53, -0.0264095520754134848538),
0607          BOOST_MATH_BIG_CONSTANT(T, 53, -0.00761224003005476438412),
0608          BOOST_MATH_BIG_CONSTANT(T, 53, -0.00247496209592143627977),
0609          BOOST_MATH_BIG_CONSTANT(T, 53, -0.000374885917942100256775),
0610          BOOST_MATH_BIG_CONSTANT(T, 53, -0.554086272024881826253e-4),
0611          BOOST_MATH_BIG_CONSTANT(T, 53, -0.396487648924804510056e-5)
0612       };
0613       BOOST_MATH_STATIC const T Q[8] = {
0614          BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),
0615          BOOST_MATH_BIG_CONSTANT(T, 53, 0.744625566823272107711),
0616          BOOST_MATH_BIG_CONSTANT(T, 53, 0.329061095011767059236),
0617          BOOST_MATH_BIG_CONSTANT(T, 53, 0.100128624977313872323),
0618          BOOST_MATH_BIG_CONSTANT(T, 53, 0.0223851099128506347278),
0619          BOOST_MATH_BIG_CONSTANT(T, 53, 0.00365334190742316650106),
0620          BOOST_MATH_BIG_CONSTANT(T, 53, 0.000402453408512476836472),
0621          BOOST_MATH_BIG_CONSTANT(T, 53, 0.263649630720255691787e-4)
0622       };
0623       T t = z / 2 - 4;
0624       result = Y + tools::evaluate_polynomial(P, t)
0625          / tools::evaluate_polynomial(Q, t);
0626       result *= exp(z) / z;
0627       result += z;
0628    }
0629    else if(z <= 20)
0630    {
0631       // Maximum Deviation Found:                     1.843e-17
0632       // Expected Error Term:                         -1.842e-17
0633       // Max Error found at double precision =        Poly: 4.375868e-17   Cheb: 5.860967e-17
0634
0635       BOOST_MATH_STATIC_LOCAL_VARIABLE const T Y = 1.0869731903076171875F;
0636       BOOST_MATH_STATIC const T P[9] = {
0637          BOOST_MATH_BIG_CONSTANT(T, 53, -0.00893891094356945667451),
0638          BOOST_MATH_BIG_CONSTANT(T, 53, -0.0484607730127134045806),
0639          BOOST_MATH_BIG_CONSTANT(T, 53, -0.0652810444222236895772),
0640          BOOST_MATH_BIG_CONSTANT(T, 53, -0.0478447572647309671455),
0641          BOOST_MATH_BIG_CONSTANT(T, 53, -0.0226059218923777094596),
0642          BOOST_MATH_BIG_CONSTANT(T, 53, -0.00720603636917482065907),
0643          BOOST_MATH_BIG_CONSTANT(T, 53, -0.00155941947035972031334),
0644          BOOST_MATH_BIG_CONSTANT(T, 53, -0.000209750022660200888349),
0645          BOOST_MATH_BIG_CONSTANT(T, 53, -0.138652200349182596186e-4)
0646       };
0647       BOOST_MATH_STATIC const T Q[9] = {
0648          BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),
0649          BOOST_MATH_BIG_CONSTANT(T, 53, 1.97017214039061194971),
0650          BOOST_MATH_BIG_CONSTANT(T, 53, 1.86232465043073157508),
0651          BOOST_MATH_BIG_CONSTANT(T, 53, 1.09601437090337519977),
0652          BOOST_MATH_BIG_CONSTANT(T, 53, 0.438873285773088870812),
0653          BOOST_MATH_BIG_CONSTANT(T, 53, 0.122537731979686102756),
0654          BOOST_MATH_BIG_CONSTANT(T, 53, 0.0233458478275769288159),
0655          BOOST_MATH_BIG_CONSTANT(T, 53, 0.00278170769163303669021),
0656          BOOST_MATH_BIG_CONSTANT(T, 53, 0.000159150281166108755531)
0657       };
0658       T t = z / 5 - 3;
0659       result = Y + tools::evaluate_polynomial(P, t)
0660          / tools::evaluate_polynomial(Q, t);
0661       result *= exp(z) / z;
0662       result += z;
0663    }
0664    else if(z <= 40)
0665    {
0666       // Maximum Deviation Found:                     5.102e-18
0667       // Expected Error Term:                         5.101e-18
0668       // Max Error found at double precision =        Poly: 1.441088e-16   Cheb: 1.864792e-16
0669
0670
0671       BOOST_MATH_STATIC_LOCAL_VARIABLE const T Y = 1.03937530517578125F;
0672       BOOST_MATH_STATIC const T P[9] = {
0673          BOOST_MATH_BIG_CONSTANT(T, 53, -0.00356165148914447597995),
0674          BOOST_MATH_BIG_CONSTANT(T, 53, -0.0229930320357982333406),
0675          BOOST_MATH_BIG_CONSTANT(T, 53, -0.0449814350482277917716),
0676          BOOST_MATH_BIG_CONSTANT(T, 53, -0.0453759383048193402336),
0677          BOOST_MATH_BIG_CONSTANT(T, 53, -0.0272050837209380717069),
0678          BOOST_MATH_BIG_CONSTANT(T, 53, -0.00994403059883350813295),
0679          BOOST_MATH_BIG_CONSTANT(T, 53, -0.00207592267812291726961),
0680          BOOST_MATH_BIG_CONSTANT(T, 53, -0.000192178045857733706044),
0681          BOOST_MATH_BIG_CONSTANT(T, 53, -0.113161784705911400295e-9)
0682       };
0683       BOOST_MATH_STATIC const T Q[9] = {
0684          BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),
0685          BOOST_MATH_BIG_CONSTANT(T, 53, 2.84354408840148561131),
0686          BOOST_MATH_BIG_CONSTANT(T, 53, 3.6599610090072393012),
0687          BOOST_MATH_BIG_CONSTANT(T, 53, 2.75088464344293083595),
0688          BOOST_MATH_BIG_CONSTANT(T, 53, 1.2985244073998398643),
0689          BOOST_MATH_BIG_CONSTANT(T, 53, 0.383213198510794507409),
0690          BOOST_MATH_BIG_CONSTANT(T, 53, 0.0651165455496281337831),
0691          BOOST_MATH_BIG_CONSTANT(T, 53, 0.00488071077519227853585)
0692       };
0693       T t = z / 10 - 3;
0694       result = Y + tools::evaluate_polynomial(P, t)
0695          / tools::evaluate_polynomial(Q, t);
0696       result *= exp(z) / z;
0697       result += z;
0698    }
0699    else
0700    {
0701       // Max Error found at double precision =        3.381886e-17
0702       BOOST_MATH_STATIC_LOCAL_VARIABLE const T exp40 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 53, 2.35385266837019985407899910749034804508871617254555467236651e17));
0703       BOOST_MATH_STATIC_LOCAL_VARIABLE const T Y= 1.013065338134765625F;
0704       BOOST_MATH_STATIC const T P[6] = {
0705          BOOST_MATH_BIG_CONSTANT(T, 53, -0.0130653381347656243849),
0706          BOOST_MATH_BIG_CONSTANT(T, 53, 0.19029710559486576682),
0707          BOOST_MATH_BIG_CONSTANT(T, 53, 94.7365094537197236011),
0708          BOOST_MATH_BIG_CONSTANT(T, 53, -2516.35323679844256203),
0709          BOOST_MATH_BIG_CONSTANT(T, 53, 18932.0850014925993025),
0710          BOOST_MATH_BIG_CONSTANT(T, 53, -38703.1431362056714134)
0711       };
0712       BOOST_MATH_STATIC const T Q[7] = {
0713          BOOST_MATH_BIG_CONSTANT(T, 53, 1.0),
0714          BOOST_MATH_BIG_CONSTANT(T, 53, 61.9733592849439884145),
0715          BOOST_MATH_BIG_CONSTANT(T, 53, -2354.56211323420194283),
0716          BOOST_MATH_BIG_CONSTANT(T, 53, 22329.1459489893079041),
0717          BOOST_MATH_BIG_CONSTANT(T, 53, -70126.245140396567133),
0718          BOOST_MATH_BIG_CONSTANT(T, 53, 54738.2833147775537106),
0719          BOOST_MATH_BIG_CONSTANT(T, 53, 8297.16296356518409347)
0720       };
0721       T t = 1 / z;
0722       result = Y + tools::evaluate_polynomial(P, t)
0723          / tools::evaluate_polynomial(Q, t);
0724       if(z < 41)
0725          result *= exp(z) / z;
0726       else
0727       {
0728          // Avoid premature overflow if we can:
0729          t = z - 40;
0730          if(t > tools::log_max_value<T>())
0731          {
0732             result = policies::raise_overflow_error<T>(function, nullptr, pol);
0733          }
0734          else
0735          {
0736             result *= exp(z - 40) / z;
0737             if(result > tools::max_value<T>() / exp40)
0738             {
0739                result = policies::raise_overflow_error<T>(function, nullptr, pol);
0740             }
0741             else
0742             {
0743                result *= exp40;
0744             }
0745          }
0746       }
0747       result += z;
0748    }
0749    return result;
0750 }
0751
0752 template <class T, class Policy>
0753 BOOST_MATH_GPU_ENABLED T expint_i_imp(T z, const Policy& pol, const boost::math::integral_constant<int, 64>& tag)
0754 {
0755    BOOST_MATH_STD_USING
0756    constexpr auto function = "boost::math::expint<%1%>(%1%)";
0757    if(z < 0)
0758       return -expint_imp(1, T(-z), pol, tag);
0759    if(z == 0)
0760       return -policies::raise_overflow_error<T>(function, nullptr, pol);
0761
0762    T result;
0763
0764    if(z <= 6)
0765    {
0766       // Maximum Deviation Found:                     3.883e-21
0767       // Expected Error Term:                         3.883e-21
0768       // Max Error found at long double precision =   Poly: 3.344801e-19   Cheb: 4.989937e-19
0769
0770       static const T P[11] = {
0771          BOOST_MATH_BIG_CONSTANT(T, 64, 2.98677224343598593764),
0772          BOOST_MATH_BIG_CONSTANT(T, 64, 0.25891613550886736592),
0773          BOOST_MATH_BIG_CONSTANT(T, 64, 0.789323584998672832285),
0774          BOOST_MATH_BIG_CONSTANT(T, 64, 0.092432587824602399339),
0775          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0514236978728625906656),
0776          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00658477469745132977921),
0777          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00124914538197086254233),
0778          BOOST_MATH_BIG_CONSTANT(T, 64, 0.000131429679565472408551),
0779          BOOST_MATH_BIG_CONSTANT(T, 64, 0.11293331317982763165e-4),
0780          BOOST_MATH_BIG_CONSTANT(T, 64, 0.629499283139417444244e-6),
0781          BOOST_MATH_BIG_CONSTANT(T, 64, 0.177833045143692498221e-7)
0782       };
0783       static const T Q[9] = {
0784          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0785          BOOST_MATH_BIG_CONSTANT(T, 64, -1.20352377969742325748),
0786          BOOST_MATH_BIG_CONSTANT(T, 64, 0.66707904942606479811),
0787          BOOST_MATH_BIG_CONSTANT(T, 64, -0.223014531629140771914),
0788          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0493340022262908008636),
0789          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00741934273050807310677),
0790          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00074353567782087939294),
0791          BOOST_MATH_BIG_CONSTANT(T, 64, -0.455861727069603367656e-4),
0792          BOOST_MATH_BIG_CONSTANT(T, 64, 0.131515429329812837701e-5)
0793       };
0794
0795       static const T c1 = BOOST_MATH_BIG_CONSTANT(T, 64, 1677624236387711.0);
0796       static const T c2 = BOOST_MATH_BIG_CONSTANT(T, 64, 4503599627370496.0);
0797       static const T r1 = c1 / c2;
0798       static const T r2 = BOOST_MATH_BIG_CONSTANT(T, 64, 0.131401834143860282009280387409357165515556574352422001206362e-16);
0799       static const T r = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.372507410781366634461991866580119133535689497771654051555657435242200120636201854384926049951548942392));
0800       T t = (z / 3) - 1;
0801       result = tools::evaluate_polynomial(P, t)
0802          / tools::evaluate_polynomial(Q, t);
0803       t = (z - r1) - r2;
0804       result *= t;
0805       if(fabs(t) < T(0.1))
0806       {
0807          result += boost::math::log1p(t / r, pol);
0808       }
0809       else
0810       {
0811          result += log(z / r);
0812       }
0813    }
0814    else if (z <= 10)
0815    {
0816       // Maximum Deviation Found:                     2.622e-21
0817       // Expected Error Term:                         -2.622e-21
0818       // Max Error found at long double precision =   Poly: 1.208328e-20   Cheb: 1.073723e-20
0819
0820       static const T Y = 1.158985137939453125F;
0821       static const T P[9] = {
0822          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00139324086199409049399),
0823          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0345238388952337563247),
0824          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0382065278072592940767),
0825          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0156117003070560727392),
0826          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00383276012430495387102),
0827          BOOST_MATH_BIG_CONSTANT(T, 64, -0.000697070540945496497992),
0828          BOOST_MATH_BIG_CONSTANT(T, 64, -0.877310384591205930343e-4),
0829          BOOST_MATH_BIG_CONSTANT(T, 64, -0.623067256376494930067e-5),
0830          BOOST_MATH_BIG_CONSTANT(T, 64, -0.377246883283337141444e-6)
0831       };
0832       static const T Q[10] = {
0833          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0834          BOOST_MATH_BIG_CONSTANT(T, 64, 1.08073635708902053767),
0835          BOOST_MATH_BIG_CONSTANT(T, 64, 0.553681133533942532909),
0836          BOOST_MATH_BIG_CONSTANT(T, 64, 0.176763647137553797451),
0837          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0387891748253869928121),
0838          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0060603004848394727017),
0839          BOOST_MATH_BIG_CONSTANT(T, 64, 0.000670519492939992806051),
0840          BOOST_MATH_BIG_CONSTANT(T, 64, 0.4947357050100855646e-4),
0841          BOOST_MATH_BIG_CONSTANT(T, 64, 0.204339282037446434827e-5),
0842          BOOST_MATH_BIG_CONSTANT(T, 64, 0.146951181174930425744e-7)
0843       };
0844       T t = z / 2 - 4;
0845       result = Y + tools::evaluate_polynomial(P, t)
0846          / tools::evaluate_polynomial(Q, t);
0847       result *= exp(z) / z;
0848       result += z;
0849    }
0850    else if(z <= 20)
0851    {
0852       // Maximum Deviation Found:                     3.220e-20
0853       // Expected Error Term:                         3.220e-20
0854       // Max Error found at long double precision =   Poly: 7.696841e-20   Cheb: 6.205163e-20
0855
0856
0857       static const T Y = 1.0869731903076171875F;
0858       static const T P[10] = {
0859          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00893891094356946995368),
0860          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0487562980088748775943),
0861          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0670568657950041926085),
0862          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0509577352851442932713),
0863          BOOST_MATH_BIG_CONSTANT(T, 64, -0.02551800927409034206),
0864          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00892913759760086687083),
0865          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00224469630207344379888),
0866          BOOST_MATH_BIG_CONSTANT(T, 64, -0.000392477245911296982776),
0867          BOOST_MATH_BIG_CONSTANT(T, 64, -0.44424044184395578775e-4),
0868          BOOST_MATH_BIG_CONSTANT(T, 64, -0.252788029251437017959e-5)
0869       };
0870       static const T Q[10] = {
0871          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0872          BOOST_MATH_BIG_CONSTANT(T, 64, 2.00323265503572414261),
0873          BOOST_MATH_BIG_CONSTANT(T, 64, 1.94688958187256383178),
0874          BOOST_MATH_BIG_CONSTANT(T, 64, 1.19733638134417472296),
0875          BOOST_MATH_BIG_CONSTANT(T, 64, 0.513137726038353385661),
0876          BOOST_MATH_BIG_CONSTANT(T, 64, 0.159135395578007264547),
0877          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0358233587351620919881),
0878          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0056716655597009417875),
0879          BOOST_MATH_BIG_CONSTANT(T, 64, 0.000577048986213535829925),
0880          BOOST_MATH_BIG_CONSTANT(T, 64, 0.290976943033493216793e-4)
0881       };
0882       T t = z / 5 - 3;
0883       result = Y + tools::evaluate_polynomial(P, t)
0884          / tools::evaluate_polynomial(Q, t);
0885       result *= exp(z) / z;
0886       result += z;
0887    }
0888    else if(z <= 40)
0889    {
0890       // Maximum Deviation Found:                     2.940e-21
0891       // Expected Error Term:                         -2.938e-21
0892       // Max Error found at long double precision =   Poly: 3.419893e-19   Cheb: 3.359874e-19
0893
0894       static const T Y = 1.03937530517578125F;
0895       static const T P[12] = {
0896          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00356165148914447278177),
0897          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0240235006148610849678),
0898          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0516699967278057976119),
0899          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0586603078706856245674),
0900          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0409960120868776180825),
0901          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0185485073689590665153),
0902          BOOST_MATH_BIG_CONSTANT(T, 64, -0.00537842101034123222417),
0903          BOOST_MATH_BIG_CONSTANT(T, 64, -0.000920988084778273760609),
0904          BOOST_MATH_BIG_CONSTANT(T, 64, -0.716742618812210980263e-4),
0905          BOOST_MATH_BIG_CONSTANT(T, 64, -0.504623302166487346677e-9),
0906          BOOST_MATH_BIG_CONSTANT(T, 64, 0.712662196671896837736e-10),
0907          BOOST_MATH_BIG_CONSTANT(T, 64, -0.533769629702262072175e-11)
0908       };
0909       static const T Q[9] = {
0910          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0911          BOOST_MATH_BIG_CONSTANT(T, 64, 3.13286733695729715455),
0912          BOOST_MATH_BIG_CONSTANT(T, 64, 4.49281223045653491929),
0913          BOOST_MATH_BIG_CONSTANT(T, 64, 3.84900294427622911374),
0914          BOOST_MATH_BIG_CONSTANT(T, 64, 2.15205199043580378211),
0915          BOOST_MATH_BIG_CONSTANT(T, 64, 0.802912186540269232424),
0916          BOOST_MATH_BIG_CONSTANT(T, 64, 0.194793170017818925388),
0917          BOOST_MATH_BIG_CONSTANT(T, 64, 0.0280128013584653182994),
0918          BOOST_MATH_BIG_CONSTANT(T, 64, 0.00182034930799902922549)
0919       };
0920       T t = z / 10 - 3;
0921       result = Y + tools::evaluate_polynomial(P, t)
0922          / tools::evaluate_polynomial(Q, t);
0923       BOOST_MATH_INSTRUMENT_VARIABLE(result)
0924       result *= exp(z) / z;
0925       BOOST_MATH_INSTRUMENT_VARIABLE(result)
0926       result += z;
0927       BOOST_MATH_INSTRUMENT_VARIABLE(result)
0928    }
0929    else
0930    {
0931       // Maximum Deviation Found:                     3.536e-20
0932       // Max Error found at long double precision =   Poly: 1.310671e-19   Cheb: 8.630943e-11
0933
0934       static const T exp40 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.35385266837019985407899910749034804508871617254555467236651e17));
0935       static const T Y= 1.013065338134765625F;
0936       static const T P[9] = {
0937          BOOST_MATH_BIG_CONSTANT(T, 64, -0.0130653381347656250004),
0938          BOOST_MATH_BIG_CONSTANT(T, 64, 0.644487780349757303739),
0939          BOOST_MATH_BIG_CONSTANT(T, 64, 143.995670348227433964),
0940          BOOST_MATH_BIG_CONSTANT(T, 64, -13918.9322758014173709),
0941          BOOST_MATH_BIG_CONSTANT(T, 64, 476260.975133624194484),
0942          BOOST_MATH_BIG_CONSTANT(T, 64, -7437102.15135982802122),
0943          BOOST_MATH_BIG_CONSTANT(T, 64, 53732298.8764767916542),
0944          BOOST_MATH_BIG_CONSTANT(T, 64, -160695051.957997452509),
0945          BOOST_MATH_BIG_CONSTANT(T, 64, 137839271.592778020028)
0946       };
0947       static const T Q[9] = {
0948          BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
0949          BOOST_MATH_BIG_CONSTANT(T, 64, 27.2103343964943718802),
0950          BOOST_MATH_BIG_CONSTANT(T, 64, -8785.48528692879413676),
0951          BOOST_MATH_BIG_CONSTANT(T, 64, 397530.290000322626766),
0952          BOOST_MATH_BIG_CONSTANT(T, 64, -7356441.34957799368252),
0953          BOOST_MATH_BIG_CONSTANT(T, 64, 63050914.5343400957524),
0954          BOOST_MATH_BIG_CONSTANT(T, 64, -246143779.638307701369),
0955          BOOST_MATH_BIG_CONSTANT(T, 64, 384647824.678554961174),
0956          BOOST_MATH_BIG_CONSTANT(T, 64, -166288297.874583961493)
0957       };
0958       T t = 1 / z;
0959       result = Y + tools::evaluate_polynomial(P, t)
0960          / tools::evaluate_polynomial(Q, t);
0961       if(z < 41)
0962          result *= exp(z) / z;
0963       else
0964       {
0965          // Avoid premature overflow if we can:
0966          t = z - 40;
0967          if(t > tools::log_max_value<T>())
0968          {
0969             result = policies::raise_overflow_error<T>(function, nullptr, pol);
0970          }
0971          else
0972          {
0973             result *= exp(z - 40) / z;
0974             if(result > tools::max_value<T>() / exp40)
0975             {
0976                result = policies::raise_overflow_error<T>(function, nullptr, pol);
0977             }
0978             else
0979             {
0980                result *= exp40;
0981             }
0982          }
0983       }
0984       result += z;
0985    }
0986    return result;
0987 }
0988
0989 template <class T, class Policy>
0990 BOOST_MATH_GPU_ENABLED void expint_i_imp_113a(T& result, const T& z, const Policy& pol)
0991 {
0992    BOOST_MATH_STD_USING
0993    // Maximum Deviation Found:                     1.230e-36
0994    // Expected Error Term:                         -1.230e-36
0995    // Max Error found at long double precision =   Poly: 4.355299e-34   Cheb: 7.512581e-34
0996
0997
0998    static const T P[15] = {
0999       BOOST_MATH_BIG_CONSTANT(T, 113, 2.98677224343598593765287235997328555),
1000       BOOST_MATH_BIG_CONSTANT(T, 113, -0.333256034674702967028780537349334037),
1001       BOOST_MATH_BIG_CONSTANT(T, 113, 0.851831522798101228384971644036708463),
1002       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0657854833494646206186773614110374948),
1003       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0630065662557284456000060708977935073),
1004       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00311759191425309373327784154659649232),
1005       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00176213568201493949664478471656026771),
1006       BOOST_MATH_BIG_CONSTANT(T, 113, -0.491548660404172089488535218163952295e-4),
1007       BOOST_MATH_BIG_CONSTANT(T, 113, 0.207764227621061706075562107748176592e-4),
1008       BOOST_MATH_BIG_CONSTANT(T, 113, -0.225445398156913584846374273379402765e-6),
1009       BOOST_MATH_BIG_CONSTANT(T, 113, 0.996939977231410319761273881672601592e-7),
1010       BOOST_MATH_BIG_CONSTANT(T, 113, 0.212546902052178643330520878928100847e-9),
1011       BOOST_MATH_BIG_CONSTANT(T, 113, 0.154646053060262871360159325115980023e-9),
1012       BOOST_MATH_BIG_CONSTANT(T, 113, 0.143971277122049197323415503594302307e-11),
1013       BOOST_MATH_BIG_CONSTANT(T, 113, 0.306243138978114692252817805327426657e-13)
1014    };
1015    static const T Q[15] = {
1016       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
1017       BOOST_MATH_BIG_CONSTANT(T, 113, -1.40178870313943798705491944989231793),
1018       BOOST_MATH_BIG_CONSTANT(T, 113, 0.943810968269701047641218856758605284),
1019       BOOST_MATH_BIG_CONSTANT(T, 113, -0.405026631534345064600850391026113165),
1020       BOOST_MATH_BIG_CONSTANT(T, 113, 0.123924153524614086482627660399122762),
1021       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0286364505373369439591132549624317707),
1022       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00516148845910606985396596845494015963),
1023       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000738330799456364820380739850924783649),
1024       BOOST_MATH_BIG_CONSTANT(T, 113, 0.843737760991856114061953265870882637e-4),
1025       BOOST_MATH_BIG_CONSTANT(T, 113, -0.767957673431982543213661388914587589e-5),
1026       BOOST_MATH_BIG_CONSTANT(T, 113, 0.549136847313854595809952100614840031e-6),
1027       BOOST_MATH_BIG_CONSTANT(T, 113, -0.299801381513743676764008325949325404e-7),
1028       BOOST_MATH_BIG_CONSTANT(T, 113, 0.118419479055346106118129130945423483e-8),
1029       BOOST_MATH_BIG_CONSTANT(T, 113, -0.30372295663095470359211949045344607e-10),
1030       BOOST_MATH_BIG_CONSTANT(T, 113, 0.382742953753485333207877784720070523e-12)
1031    };
1032
1033    static const T c1 = BOOST_MATH_BIG_CONSTANT(T, 113, 1677624236387711.0);
1034    static const T c2 = BOOST_MATH_BIG_CONSTANT(T, 113, 4503599627370496.0);
1035    static const T c3 = BOOST_MATH_BIG_CONSTANT(T, 113, 266514582277687.0);
1036    static const T c4 = BOOST_MATH_BIG_CONSTANT(T, 113, 4503599627370496.0);
1037    static const T c5 = BOOST_MATH_BIG_CONSTANT(T, 113, 4503599627370496.0);
1038    static const T r1 = c1 / c2;
1039    static const T r2 = c3 / c4 / c5;
1040    static const T r3 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, 0.283806480836357377069325311780969887585024578164571984232357e-31));
1041    static const T r = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, 0.372507410781366634461991866580119133535689497771654051555657435242200120636201854384926049951548942392));
1042    T t = (z / 3) - 1;
1043    result = tools::evaluate_polynomial(P, t)
1044       / tools::evaluate_polynomial(Q, t);
1045    t = ((z - r1) - r2) - r3;
1046    result *= t;
1047    if(fabs(t) < 0.1)
1048    {
1049       result += boost::math::log1p(t / r, pol);
1050    }
1051    else
1052    {
1053       result += log(z / r);
1054    }
1055 }
1056
1057 template <class T>
1058 BOOST_MATH_GPU_ENABLED void expint_i_113b(T& result, const T& z)
1059 {
1060    BOOST_MATH_STD_USING
1061    // Maximum Deviation Found:                     7.779e-36
1062    // Expected Error Term:                         -7.779e-36
1063    // Max Error found at long double precision =   Poly: 2.576723e-35   Cheb: 1.236001e-34
1064
1065    static const T Y = 1.158985137939453125F;
1066    static const T P[15] = {
1067       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00139324086199409049282472239613554817),
1068       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0338173111691991289178779840307998955),
1069       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0555972290794371306259684845277620556),
1070       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0378677976003456171563136909186202177),
1071       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0152221583517528358782902783914356667),
1072       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00428283334203873035104248217403126905),
1073       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000922782631491644846511553601323435286),
1074       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000155513428088853161562660696055496696),
1075       BOOST_MATH_BIG_CONSTANT(T, 113, -0.205756580255359882813545261519317096e-4),
1076       BOOST_MATH_BIG_CONSTANT(T, 113, -0.220327406578552089820753181821115181e-5),
1077       BOOST_MATH_BIG_CONSTANT(T, 113, -0.189483157545587592043421445645377439e-6),
1078       BOOST_MATH_BIG_CONSTANT(T, 113, -0.122426571518570587750898968123803867e-7),
1079       BOOST_MATH_BIG_CONSTANT(T, 113, -0.635187358949437991465353268374523944e-9),
1080       BOOST_MATH_BIG_CONSTANT(T, 113, -0.203015132965870311935118337194860863e-10),
1081       BOOST_MATH_BIG_CONSTANT(T, 113, -0.384276705503357655108096065452950822e-12)
1082    };
1083    static const T Q[15] = {
1084       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
1085       BOOST_MATH_BIG_CONSTANT(T, 113, 1.58784732785354597996617046880946257),
1086       BOOST_MATH_BIG_CONSTANT(T, 113, 1.18550755302279446339364262338114098),
1087       BOOST_MATH_BIG_CONSTANT(T, 113, 0.55598993549661368604527040349702836),
1088       BOOST_MATH_BIG_CONSTANT(T, 113, 0.184290888380564236919107835030984453),
1089       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0459658051803613282360464632326866113),
1090       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0089505064268613225167835599456014705),
1091       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00139042673882987693424772855926289077),
1092       BOOST_MATH_BIG_CONSTANT(T, 113, 0.000174210708041584097450805790176479012),
1093       BOOST_MATH_BIG_CONSTANT(T, 113, 0.176324034009707558089086875136647376e-4),
1094       BOOST_MATH_BIG_CONSTANT(T, 113, 0.142935845999505649273084545313710581e-5),
1095       BOOST_MATH_BIG_CONSTANT(T, 113, 0.907502324487057260675816233312747784e-7),
1096       BOOST_MATH_BIG_CONSTANT(T, 113, 0.431044337808893270797934621235918418e-8),
1097       BOOST_MATH_BIG_CONSTANT(T, 113, 0.139007266881450521776529705677086902e-9),
1098       BOOST_MATH_BIG_CONSTANT(T, 113, 0.234715286125516430792452741830364672e-11)
1099    };
1100    T t = z / 2 - 4;
1101    result = Y + tools::evaluate_polynomial(P, t)
1102       / tools::evaluate_polynomial(Q, t);
1103    result *= exp(z) / z;
1104    result += z;
1105 }
1106
1107 template <class T>
1108 BOOST_MATH_GPU_ENABLED void expint_i_113c(T& result, const T& z)
1109 {
1110    BOOST_MATH_STD_USING
1111    // Maximum Deviation Found:                     1.082e-34
1112    // Expected Error Term:                         1.080e-34
1113    // Max Error found at long double precision =   Poly: 1.958294e-34   Cheb: 2.472261e-34
1114
1115
1116    static const T Y = 1.091579437255859375F;
1117    static const T P[17] = {
1118       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00685089599550151282724924894258520532),
1119       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0443313550253580053324487059748497467),
1120       BOOST_MATH_BIG_CONSTANT(T, 113, -0.071538561252424027443296958795814874),
1121       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0622923153354102682285444067843300583),
1122       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0361631270264607478205393775461208794),
1123       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0153192826839624850298106509601033261),
1124       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00496967904961260031539602977748408242),
1125       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00126989079663425780800919171538920589),
1126       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000258933143097125199914724875206326698),
1127       BOOST_MATH_BIG_CONSTANT(T, 113, -0.422110326689204794443002330541441956e-4),
1128       BOOST_MATH_BIG_CONSTANT(T, 113, -0.546004547590412661451073996127115221e-5),
1129       BOOST_MATH_BIG_CONSTANT(T, 113, -0.546775260262202177131068692199272241e-6),
1130       BOOST_MATH_BIG_CONSTANT(T, 113, -0.404157632825805803833379568956559215e-7),
1131       BOOST_MATH_BIG_CONSTANT(T, 113, -0.200612596196561323832327013027419284e-8),
1132       BOOST_MATH_BIG_CONSTANT(T, 113, -0.502538501472133913417609379765434153e-10),
1133       BOOST_MATH_BIG_CONSTANT(T, 113, -0.326283053716799774936661568391296584e-13),
1134       BOOST_MATH_BIG_CONSTANT(T, 113, 0.869226483473172853557775877908693647e-15)
1135    };
1136    static const T Q[15] = {
1137       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
1138       BOOST_MATH_BIG_CONSTANT(T, 113, 2.23227220874479061894038229141871087),
1139       BOOST_MATH_BIG_CONSTANT(T, 113, 2.40221000361027971895657505660959863),
1140       BOOST_MATH_BIG_CONSTANT(T, 113, 1.65476320985936174728238416007084214),
1141       BOOST_MATH_BIG_CONSTANT(T, 113, 0.816828602963895720369875535001248227),
1142       BOOST_MATH_BIG_CONSTANT(T, 113, 0.306337922909446903672123418670921066),
1143       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0902400121654409267774593230720600752),
1144       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0212708882169429206498765100993228086),
1145       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00404442626252467471957713495828165491),
1146       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0006195601618842253612635241404054589),
1147       BOOST_MATH_BIG_CONSTANT(T, 113, 0.755930932686543009521454653994321843e-4),
1148       BOOST_MATH_BIG_CONSTANT(T, 113, 0.716004532773778954193609582677482803e-5),
1149       BOOST_MATH_BIG_CONSTANT(T, 113, 0.500881663076471627699290821742924233e-6),
1150       BOOST_MATH_BIG_CONSTANT(T, 113, 0.233593219218823384508105943657387644e-7),
1151       BOOST_MATH_BIG_CONSTANT(T, 113, 0.554900353169148897444104962034267682e-9)
1152    };
1153    T t = z / 4 - 3.5;
1154    result = Y + tools::evaluate_polynomial(P, t)
1155       / tools::evaluate_polynomial(Q, t);
1156    result *= exp(z) / z;
1157    result += z;
1158 }
1159
1160 template <class T>
1161 BOOST_MATH_GPU_ENABLED void expint_i_113d(T& result, const T& z)
1162 {
1163    BOOST_MATH_STD_USING
1164    // Maximum Deviation Found:                     3.163e-35
1165    // Expected Error Term:                         3.163e-35
1166    // Max Error found at long double precision =   Poly: 4.158110e-35   Cheb: 5.385532e-35
1167
1168    static const T Y = 1.051731109619140625F;
1169    static const T P[14] = {
1170       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00144552494420652573815404828020593565),
1171       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0126747451594545338365684731262912741),
1172       BOOST_MATH_BIG_CONSTANT(T, 113, -0.01757394877502366717526779263438073),
1173       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0126838952395506921945756139424722588),
1174       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0060045057928894974954756789352443522),
1175       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00205349237147226126653803455793107903),
1176       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000532606040579654887676082220195624207),
1177       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000107344687098019891474772069139014662),
1178       BOOST_MATH_BIG_CONSTANT(T, 113, -0.169536802705805811859089949943435152e-4),
1179       BOOST_MATH_BIG_CONSTANT(T, 113, -0.20863311729206543881826553010120078e-5),
1180       BOOST_MATH_BIG_CONSTANT(T, 113, -0.195670358542116256713560296776654385e-6),
1181       BOOST_MATH_BIG_CONSTANT(T, 113, -0.133291168587253145439184028259772437e-7),
1182       BOOST_MATH_BIG_CONSTANT(T, 113, -0.595500337089495614285777067722823397e-9),
1183       BOOST_MATH_BIG_CONSTANT(T, 113, -0.133141358866324100955927979606981328e-10)
1184    };
1185    static const T Q[14] = {
1186       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
1187       BOOST_MATH_BIG_CONSTANT(T, 113, 1.72490783907582654629537013560044682),
1188       BOOST_MATH_BIG_CONSTANT(T, 113, 1.44524329516800613088375685659759765),
1189       BOOST_MATH_BIG_CONSTANT(T, 113, 0.778241785539308257585068744978050181),
1190       BOOST_MATH_BIG_CONSTANT(T, 113, 0.300520486589206605184097270225725584),
1191       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0879346899691339661394537806057953957),
1192       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0200802415843802892793583043470125006),
1193       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00362842049172586254520256100538273214),
1194       BOOST_MATH_BIG_CONSTANT(T, 113, 0.000519731362862955132062751246769469957),
1195       BOOST_MATH_BIG_CONSTANT(T, 113, 0.584092147914050999895178697392282665e-4),
1196       BOOST_MATH_BIG_CONSTANT(T, 113, 0.501851497707855358002773398333542337e-5),
1197       BOOST_MATH_BIG_CONSTANT(T, 113, 0.313085677467921096644895738538865537e-6),
1198       BOOST_MATH_BIG_CONSTANT(T, 113, 0.127552010539733113371132321521204458e-7),
1199       BOOST_MATH_BIG_CONSTANT(T, 113, 0.25737310826983451144405899970774587e-9)
1200    };
1201    T t = z / 4 - 5.5;
1202    result = Y + tools::evaluate_polynomial(P, t)
1203       / tools::evaluate_polynomial(Q, t);
1204    BOOST_MATH_INSTRUMENT_VARIABLE(result)
1205    result *= exp(z) / z;
1206    BOOST_MATH_INSTRUMENT_VARIABLE(result)
1207    result += z;
1208    BOOST_MATH_INSTRUMENT_VARIABLE(result)
1209 }
1210
1211 template <class T>
1212 BOOST_MATH_GPU_ENABLED void expint_i_113e(T& result, const T& z)
1213 {
1214    BOOST_MATH_STD_USING
1215    // Maximum Deviation Found:                     7.972e-36
1216    // Expected Error Term:                         7.962e-36
1217    // Max Error found at long double precision =   Poly: 1.711721e-34   Cheb: 3.100018e-34
1218
1219    static const T Y = 1.032726287841796875F;
1220    static const T P[15] = {
1221       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00141056919297307534690895009969373233),
1222       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0123384175302540291339020257071411437),
1223       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0298127270706864057791526083667396115),
1224       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0390686759471630584626293670260768098),
1225       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0338226792912607409822059922949035589),
1226       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0211659736179834946452561197559654582),
1227       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0100428887460879377373158821400070313),
1228       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00370717396015165148484022792801682932),
1229       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0010768667551001624764329000496561659),
1230       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000246127328761027039347584096573123531),
1231       BOOST_MATH_BIG_CONSTANT(T, 113, -0.437318110527818613580613051861991198e-4),
1232       BOOST_MATH_BIG_CONSTANT(T, 113, -0.587532682329299591501065482317771497e-5),
1233       BOOST_MATH_BIG_CONSTANT(T, 113, -0.565697065670893984610852937110819467e-6),
1234       BOOST_MATH_BIG_CONSTANT(T, 113, -0.350233957364028523971768887437839573e-7),
1235       BOOST_MATH_BIG_CONSTANT(T, 113, -0.105428907085424234504608142258423505e-8)
1236    };
1237    static const T Q[16] = {
1238       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
1239       BOOST_MATH_BIG_CONSTANT(T, 113, 3.17261315255467581204685605414005525),
1240       BOOST_MATH_BIG_CONSTANT(T, 113, 4.85267952971640525245338392887217426),
1241       BOOST_MATH_BIG_CONSTANT(T, 113, 4.74341914912439861451492872946725151),
1242       BOOST_MATH_BIG_CONSTANT(T, 113, 3.31108463283559911602405970817931801),
1243       BOOST_MATH_BIG_CONSTANT(T, 113, 1.74657006336994649386607925179848899),
1244       BOOST_MATH_BIG_CONSTANT(T, 113, 0.718255607416072737965933040353653244),
1245       BOOST_MATH_BIG_CONSTANT(T, 113, 0.234037553177354542791975767960643864),
1246       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0607470145906491602476833515412605389),
1247       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0125048143774226921434854172947548724),
1248       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00201034366420433762935768458656609163),
1249       BOOST_MATH_BIG_CONSTANT(T, 113, 0.000244823338417452367656368849303165721),
1250       BOOST_MATH_BIG_CONSTANT(T, 113, 0.213511655166983177960471085462540807e-4),
1251       BOOST_MATH_BIG_CONSTANT(T, 113, 0.119323998465870686327170541547982932e-5),
1252       BOOST_MATH_BIG_CONSTANT(T, 113, 0.322153582559488797803027773591727565e-7),
1253       BOOST_MATH_BIG_CONSTANT(T, 113, -0.161635525318683508633792845159942312e-16)
1254    };
1255    T t = z / 8 - 4.25;
1256    result = Y + tools::evaluate_polynomial(P, t)
1257       / tools::evaluate_polynomial(Q, t);
1258    BOOST_MATH_INSTRUMENT_VARIABLE(result)
1259    result *= exp(z) / z;
1260    BOOST_MATH_INSTRUMENT_VARIABLE(result)
1261    result += z;
1262    BOOST_MATH_INSTRUMENT_VARIABLE(result)
1263 }
1264
1265 template <class T>
1266 BOOST_MATH_GPU_ENABLED void expint_i_113f(T& result, const T& z)
1267 {
1268    BOOST_MATH_STD_USING
1269    // Maximum Deviation Found:                     4.469e-36
1270    // Expected Error Term:                         4.468e-36
1271    // Max Error found at long double precision =   Poly: 1.288958e-35   Cheb: 2.304586e-35
1272
1273    static const T Y = 1.0216197967529296875F;
1274    static const T P[12] = {
1275       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000322999116096627043476023926572650045),
1276       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00385606067447365187909164609294113346),
1277       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00686514524727568176735949971985244415),
1278       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00606260649593050194602676772589601799),
1279       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00334382362017147544335054575436194357),
1280       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00126108534260253075708625583630318043),
1281       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000337881489347846058951220431209276776),
1282       BOOST_MATH_BIG_CONSTANT(T, 113, -0.648480902304640018785370650254018022e-4),
1283       BOOST_MATH_BIG_CONSTANT(T, 113, -0.87652644082970492211455290209092766e-5),
1284       BOOST_MATH_BIG_CONSTANT(T, 113, -0.794712243338068631557849449519994144e-6),
1285       BOOST_MATH_BIG_CONSTANT(T, 113, -0.434084023639508143975983454830954835e-7),
1286       BOOST_MATH_BIG_CONSTANT(T, 113, -0.107839681938752337160494412638656696e-8)
1287    };
1288    static const T Q[12] = {
1289       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
1290       BOOST_MATH_BIG_CONSTANT(T, 113, 2.09913805456661084097134805151524958),
1291       BOOST_MATH_BIG_CONSTANT(T, 113, 2.07041755535439919593503171320431849),
1292       BOOST_MATH_BIG_CONSTANT(T, 113, 1.26406517226052371320416108604874734),
1293       BOOST_MATH_BIG_CONSTANT(T, 113, 0.529689923703770353961553223973435569),
1294       BOOST_MATH_BIG_CONSTANT(T, 113, 0.159578150879536711042269658656115746),
1295       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0351720877642000691155202082629857131),
1296       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00565313621289648752407123620997063122),
1297       BOOST_MATH_BIG_CONSTANT(T, 113, 0.000646920278540515480093843570291218295),
1298       BOOST_MATH_BIG_CONSTANT(T, 113, 0.499904084850091676776993523323213591e-4),
1299       BOOST_MATH_BIG_CONSTANT(T, 113, 0.233740058688179614344680531486267142e-5),
1300       BOOST_MATH_BIG_CONSTANT(T, 113, 0.498800627828842754845418576305379469e-7)
1301    };
1302    T t = z / 7 - 7;
1303    result = Y + tools::evaluate_polynomial(P, t)
1304       / tools::evaluate_polynomial(Q, t);
1305    BOOST_MATH_INSTRUMENT_VARIABLE(result)
1306    result *= exp(z) / z;
1307    BOOST_MATH_INSTRUMENT_VARIABLE(result)
1308    result += z;
1309    BOOST_MATH_INSTRUMENT_VARIABLE(result)
1310 }
1311
1312 template <class T>
1313 BOOST_MATH_GPU_ENABLED void expint_i_113g(T& result, const T& z)
1314 {
1315    BOOST_MATH_STD_USING
1316    // Maximum Deviation Found:                     5.588e-35
1317    // Expected Error Term:                         -5.566e-35
1318    // Max Error found at long double precision =   Poly: 9.976345e-35   Cheb: 8.358865e-35
1319
1320    static const T Y = 1.015148162841796875F;
1321    static const T P[11] = {
1322       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000435714784725086961464589957142615216),
1323       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00432114324353830636009453048419094314),
1324       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0100740363285526177522819204820582424),
1325       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0116744115827059174392383504427640362),
1326       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00816145387784261141360062395898644652),
1327       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00371380272673500791322744465394211508),
1328       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00112958263488611536502153195005736563),
1329       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000228316462389404645183269923754256664),
1330       BOOST_MATH_BIG_CONSTANT(T, 113, -0.29462181955852860250359064291292577e-4),
1331       BOOST_MATH_BIG_CONSTANT(T, 113, -0.21972450610957417963227028788460299e-5),
1332       BOOST_MATH_BIG_CONSTANT(T, 113, -0.720558173805289167524715527536874694e-7)
1333    };
1334    static const T Q[11] = {
1335       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
1336       BOOST_MATH_BIG_CONSTANT(T, 113, 2.95918362458402597039366979529287095),
1337       BOOST_MATH_BIG_CONSTANT(T, 113, 3.96472247520659077944638411856748924),
1338       BOOST_MATH_BIG_CONSTANT(T, 113, 3.15563251550528513747923714884142131),
1339       BOOST_MATH_BIG_CONSTANT(T, 113, 1.64674612007093983894215359287448334),
1340       BOOST_MATH_BIG_CONSTANT(T, 113, 0.58695020129846594405856226787156424),
1341       BOOST_MATH_BIG_CONSTANT(T, 113, 0.144358385319329396231755457772362793),
1342       BOOST_MATH_BIG_CONSTANT(T, 113, 0.024146911506411684815134916238348063),
1343       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0026257132337460784266874572001650153),
1344       BOOST_MATH_BIG_CONSTANT(T, 113, 0.000167479843750859222348869769094711093),
1345       BOOST_MATH_BIG_CONSTANT(T, 113, 0.475673638665358075556452220192497036e-5)
1346    };
1347    T t = z / 14 - 5;
1348    result = Y + tools::evaluate_polynomial(P, t)
1349       / tools::evaluate_polynomial(Q, t);
1350    BOOST_MATH_INSTRUMENT_VARIABLE(result)
1351    result *= exp(z) / z;
1352    BOOST_MATH_INSTRUMENT_VARIABLE(result)
1353    result += z;
1354    BOOST_MATH_INSTRUMENT_VARIABLE(result)
1355 }
1356
1357 template <class T>
1358 BOOST_MATH_GPU_ENABLED void expint_i_113h(T& result, const T& z)
1359 {
1360    BOOST_MATH_STD_USING
1361    // Maximum Deviation Found:                     4.448e-36
1362    // Expected Error Term:                         4.445e-36
1363    // Max Error found at long double precision =   Poly: 2.058532e-35   Cheb: 2.165465e-27
1364
1365    static const T Y= 1.00849151611328125F;
1366    static const T P[9] = {
1367       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0084915161132812500000001440233607358),
1368       BOOST_MATH_BIG_CONSTANT(T, 113, 1.84479378737716028341394223076147872),
1369       BOOST_MATH_BIG_CONSTANT(T, 113, -130.431146923726715674081563022115568),
1370       BOOST_MATH_BIG_CONSTANT(T, 113, 4336.26945491571504885214176203512015),
1371       BOOST_MATH_BIG_CONSTANT(T, 113, -76279.0031974974730095170437591004177),
1372       BOOST_MATH_BIG_CONSTANT(T, 113, 729577.956271997673695191455111727774),
1373       BOOST_MATH_BIG_CONSTANT(T, 113, -3661928.69330208734947103004900349266),
1374       BOOST_MATH_BIG_CONSTANT(T, 113, 8570600.041606912735872059184527855),
1375       BOOST_MATH_BIG_CONSTANT(T, 113, -6758379.93672362080947905580906028645)
1376    };
1377    static const T Q[10] = {
1378       BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
1379       BOOST_MATH_BIG_CONSTANT(T, 113, -99.4868026047611434569541483506091713),
1380       BOOST_MATH_BIG_CONSTANT(T, 113, 3879.67753690517114249705089803055473),
1381       BOOST_MATH_BIG_CONSTANT(T, 113, -76495.82413252517165830203774900806),
1382       BOOST_MATH_BIG_CONSTANT(T, 113, 820773.726408311894342553758526282667),
1383       BOOST_MATH_BIG_CONSTANT(T, 113, -4803087.64956923577571031564909646579),
1384       BOOST_MATH_BIG_CONSTANT(T, 113, 14521246.227703545012713173740895477),
1385       BOOST_MATH_BIG_CONSTANT(T, 113, -19762752.0196769712258527849159393044),
1386       BOOST_MATH_BIG_CONSTANT(T, 113, 8354144.67882768405803322344185185517),
1387       BOOST_MATH_BIG_CONSTANT(T, 113, 355076.853106511136734454134915432571)
1388    };
1389    T t = 1 / z;
1390    result = Y + tools::evaluate_polynomial(P, t)
1391       / tools::evaluate_polynomial(Q, t);
1392    result *= exp(z) / z;
1393    result += z;
1394 }
1395
1396 template <class T, class Policy>
1397 BOOST_MATH_GPU_ENABLED T expint_i_imp(T z, const Policy& pol, const boost::math::integral_constant<int, 113>& tag)
1398 {
1399    BOOST_MATH_STD_USING
1400    constexpr auto function = "boost::math::expint<%1%>(%1%)";
1401    if(z < 0)
1402       return -expint_imp(1, T(-z), pol, tag);
1403    if(z == 0)
1404       return -policies::raise_overflow_error<T>(function, nullptr, pol);
1405
1406    T result;
1407
1408    if(z <= 6)
1409    {
1410       expint_i_imp_113a(result, z, pol);
1411    }
1412    else if (z <= 10)
1413    {
1414       expint_i_113b(result, z);
1415    }
1416    else if(z <= 18)
1417    {
1418       expint_i_113c(result, z);
1419    }
1420    else if(z <= 26)
1421    {
1422       expint_i_113d(result, z);
1423    }
1424    else if(z <= 42)
1425    {
1426       expint_i_113e(result, z);
1427    }
1428    else if(z <= 56)
1429    {
1430       expint_i_113f(result, z);
1431    }
1432    else if(z <= 84)
1433    {
1434       expint_i_113g(result, z);
1435    }
1436    else if(z <= 210)
1437    {
1438       expint_i_113h(result, z);
1439    }
1440    else // z > 210
1441    {
1442       // Maximum Deviation Found:                     3.963e-37
1443       // Expected Error Term:                         3.963e-37
1444       // Max Error found at long double precision =   Poly: 1.248049e-36   Cheb: 2.843486e-29
1445
1446       static const T exp40 = static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, 2.35385266837019985407899910749034804508871617254555467236651e17));
1447       static const T Y= 1.00252532958984375F;
1448       static const T P[8] = {
1449          BOOST_MATH_BIG_CONSTANT(T, 113, -0.00252532958984375000000000000000000085),
1450          BOOST_MATH_BIG_CONSTANT(T, 113, 1.16591386866059087390621952073890359),
1451          BOOST_MATH_BIG_CONSTANT(T, 113, -67.8483431314018462417456828499277579),
1452          BOOST_MATH_BIG_CONSTANT(T, 113, 1567.68688154683822956359536287575892),
1453          BOOST_MATH_BIG_CONSTANT(T, 113, -17335.4683325819116482498725687644986),
1454          BOOST_MATH_BIG_CONSTANT(T, 113, 93632.6567462673524739954389166550069),
1455          BOOST_MATH_BIG_CONSTANT(T, 113, -225025.189335919133214440347510936787),
1456          BOOST_MATH_BIG_CONSTANT(T, 113, 175864.614717440010942804684741336853)
1457       };
1458       static const T Q[9] = {
1459          BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
1460          BOOST_MATH_BIG_CONSTANT(T, 113, -65.6998869881600212224652719706425129),
1461          BOOST_MATH_BIG_CONSTANT(T, 113, 1642.73850032324014781607859416890077),
1462          BOOST_MATH_BIG_CONSTANT(T, 113, -19937.2610222467322481947237312818575),
1463          BOOST_MATH_BIG_CONSTANT(T, 113, 124136.267326632742667972126625064538),
1464          BOOST_MATH_BIG_CONSTANT(T, 113, -384614.251466704550678760562965502293),
1465          BOOST_MATH_BIG_CONSTANT(T, 113, 523355.035910385688578278384032026998),
1466          BOOST_MATH_BIG_CONSTANT(T, 113, -217809.552260834025885677791936351294),
1467          BOOST_MATH_BIG_CONSTANT(T, 113, -8555.81719551123640677261226549550872)
1468       };
1469       T t = 1 / z;
1470       result = Y + tools::evaluate_polynomial(P, t)
1471          / tools::evaluate_polynomial(Q, t);
1472       if(z < 41)
1473          result *= exp(z) / z;
1474       else
1475       {
1476          // Avoid premature overflow if we can:
1477          t = z - 40;
1478          if(t > tools::log_max_value<T>())
1479          {
1480             result = policies::raise_overflow_error<T>(function, nullptr, pol);
1481          }
1482          else
1483          {
1484             result *= exp(z - 40) / z;
1485             if(result > tools::max_value<T>() / exp40)
1486             {
1487                result = policies::raise_overflow_error<T>(function, nullptr, pol);
1488             }
1489             else
1490             {
1491                result *= exp40;
1492             }
1493          }
1494       }
1495       result += z;
1496    }
1497    return result;
1498 }
1499
1500 template <class T, class Policy, class tag>
1501 struct expint_i_initializer
1502 {
1503    struct init
1504    {
1505       BOOST_MATH_GPU_ENABLED init()
1506       {
1507          do_init(tag());
1508       }
1509       BOOST_MATH_GPU_ENABLED static void do_init(const boost::math::integral_constant<int, 0>&){}
1510       BOOST_MATH_GPU_ENABLED static void do_init(const boost::math::integral_constant<int, 53>&)
1511       {
1512          boost::math::expint(T(5), Policy());
1513          boost::math::expint(T(7), Policy());
1514          boost::math::expint(T(18), Policy());
1515          boost::math::expint(T(38), Policy());
1516          boost::math::expint(T(45), Policy());
1517       }
1518       BOOST_MATH_GPU_ENABLED static void do_init(const boost::math::integral_constant<int, 64>&)
1519       {
1520          boost::math::expint(T(5), Policy());
1521          boost::math::expint(T(7), Policy());
1522          boost::math::expint(T(18), Policy());
1523          boost::math::expint(T(38), Policy());
1524          boost::math::expint(T(45), Policy());
1525       }
1526       BOOST_MATH_GPU_ENABLED static void do_init(const boost::math::integral_constant<int, 113>&)
1527       {
1528          boost::math::expint(T(5), Policy());
1529          boost::math::expint(T(7), Policy());
1530          boost::math::expint(T(17), Policy());
1531          boost::math::expint(T(25), Policy());
1532          boost::math::expint(T(40), Policy());
1533          boost::math::expint(T(50), Policy());
1534          boost::math::expint(T(80), Policy());
1535          boost::math::expint(T(200), Policy());
1536          boost::math::expint(T(220), Policy());
1537       }
1538       BOOST_MATH_GPU_ENABLED void force_instantiate()const{}
1539    };
1540    static const init initializer;
1541    BOOST_MATH_GPU_ENABLED static void force_instantiate()
1542    {
1543       #ifndef BOOST_MATH_HAS_GPU_SUPPORT
1544       initializer.force_instantiate();
1545       #endif
1546    }
1547 };
1548
1549 template <class T, class Policy, class tag>
1550 const typename expint_i_initializer<T, Policy, tag>::init expint_i_initializer<T, Policy, tag>::initializer;
1551
1552 template <class T, class Policy, class tag>
1553 struct expint_1_initializer
1554 {
1555    struct init
1556    {
1557       BOOST_MATH_GPU_ENABLED init()
1558       {
1559          do_init(tag());
1560       }
1561       BOOST_MATH_GPU_ENABLED static void do_init(const boost::math::integral_constant<int, 0>&){}
1562       BOOST_MATH_GPU_ENABLED static void do_init(const boost::math::integral_constant<int, 53>&)
1563       {
1564          boost::math::expint(1, T(0.5), Policy());
1565          boost::math::expint(1, T(2), Policy());
1566       }
1567       BOOST_MATH_GPU_ENABLED static void do_init(const boost::math::integral_constant<int, 64>&)
1568       {
1569          boost::math::expint(1, T(0.5), Policy());
1570          boost::math::expint(1, T(2), Policy());
1571       }
1572       BOOST_MATH_GPU_ENABLED static void do_init(const boost::math::integral_constant<int, 113>&)
1573       {
1574          boost::math::expint(1, T(0.5), Policy());
1575          boost::math::expint(1, T(2), Policy());
1576          boost::math::expint(1, T(6), Policy());
1577       }
1578       BOOST_MATH_GPU_ENABLED void force_instantiate()const{}
1579    };
1580    static const init initializer;
1581    BOOST_MATH_GPU_ENABLED static void force_instantiate()
1582    {
1583       #ifndef BOOST_MATH_HAS_GPU_SUPPORT
1584       initializer.force_instantiate();
1585       #endif
1586    }
1587 };
1588
1589 template <class T, class Policy, class tag>
1590 const typename expint_1_initializer<T, Policy, tag>::init expint_1_initializer<T, Policy, tag>::initializer;
1591
1592 template <class T, class Policy>
1593 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type
1594    expint_forwarder(T z, const Policy& /*pol*/, boost::math::true_type const&)
1595 {
1596    typedef typename tools::promote_args<T>::type result_type;
1597    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1598    typedef typename policies::precision<result_type, Policy>::type precision_type;
1599    typedef typename policies::normalise<
1600       Policy,
1601       policies::promote_float<false>,
1602       policies::promote_double<false>,
1603       policies::discrete_quantile<>,
1604       policies::assert_undefined<> >::type forwarding_policy;
1605    typedef boost::math::integral_constant<int,
1606       precision_type::value <= 0 ? 0 :
1607       precision_type::value <= 53 ? 53 :
1608       precision_type::value <= 64 ? 64 :
1609       precision_type::value <= 113 ? 113 : 0
1610    > tag_type;
1611
1612    expint_i_initializer<value_type, forwarding_policy, tag_type>::force_instantiate();
1613
1614    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expint_i_imp(
1615       static_cast<value_type>(z),
1616       forwarding_policy(),
1617       tag_type()), "boost::math::expint<%1%>(%1%)");
1618 }
1619
1620 template <class T>
1621 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type
1622 expint_forwarder(unsigned n, T z, const boost::math::false_type&)
1623 {
1624    return boost::math::expint(n, z, policies::policy<>());
1625 }
1626
1627 } // namespace detail
1628
1629 template <class T, class Policy>
1630 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type
1631    expint(unsigned n, T z, const Policy& /*pol*/)
1632 {
1633    typedef typename tools::promote_args<T>::type result_type;
1634    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1635    typedef typename policies::precision<result_type, Policy>::type precision_type;
1636    typedef typename policies::normalise<
1637       Policy,
1638       policies::promote_float<false>,
1639       policies::promote_double<false>,
1640       policies::discrete_quantile<>,
1641       policies::assert_undefined<> >::type forwarding_policy;
1642    typedef boost::math::integral_constant<int,
1643       precision_type::value <= 0 ? 0 :
1644       precision_type::value <= 53 ? 53 :
1645       precision_type::value <= 64 ? 64 :
1646       precision_type::value <= 113 ? 113 : 0
1647    > tag_type;
1648
1649    detail::expint_1_initializer<value_type, forwarding_policy, tag_type>::force_instantiate();
1650
1651    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expint_imp(
1652       n,
1653       static_cast<value_type>(z),
1654       forwarding_policy(),
1655       tag_type()), "boost::math::expint<%1%>(unsigned, %1%)");
1656 }
1657
1658 template <class T, class U>
1659 BOOST_MATH_GPU_ENABLED inline typename detail::expint_result<T, U>::type
1660    expint(T const z, U const u)
1661 {
1662    typedef typename policies::is_policy<U>::type tag_type;
1663    return detail::expint_forwarder(z, u, tag_type());
1664 }
1665
1666 template <class T>
1667 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type
1668    expint(T z)
1669 {
1670    return expint(z, policies::policy<>());
1671 }
1672
1673 }} // namespaces
1674
1675 #ifdef _MSC_VER
1676 #pragma warning(pop)
1677 #endif
1678
1679 #endif // BOOST_MATH_EXPINT_HPP
1680
1681