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