vamp-build-and-test: DEPENDENCIES/generic/include/boost/math/special

annotate DEPENDENCIES/generic/include/boost/math/special_functions/expint.hpp @ 125:34e428693f5d vext

Vext -> Repoint

author	Chris Cannam
date	Thu, 14 Jun 2018 11:15:39 +0100
parents	c530137014c0
children

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

Mercurial > hg > vamp-build-and-test

annotate DEPENDENCIES/generic/include/boost/math/special_functions/expint.hpp @ 125:34e428693f5d vext