ceph/src/boost/libs/math/include/boost/math/special_functions/expint.hpp

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