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