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