Mercurial > hg > vamp-build-and-test

comparison DEPENDENCIES/generic/include/boost/math/special_functions/expint.hpp @ 16:2665513ce2d3

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