1// Special functions -*- C++ -*- 
2 
3// Copyright (C) 2006-2019 Free Software Foundation, Inc. 
4// 
5// This file is part of the GNU ISO C++ Library. This library is free 
6// software; you can redistribute it and/or modify it under the 
7// terms of the GNU General Public License as published by the 
8// Free Software Foundation; either version 3, or (at your option) 
9// any later version. 
10// 
11// This library is distributed in the hope that it will be useful, 
12// but WITHOUT ANY WARRANTY; without even the implied warranty of 
13// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
14// GNU General Public License for more details. 
15// 
16// Under Section 7 of GPL version 3, you are granted additional 
17// permissions described in the GCC Runtime Library Exception, version 
18// 3.1, as published by the Free Software Foundation. 
19 
20// You should have received a copy of the GNU General Public License and 
21// a copy of the GCC Runtime Library Exception along with this program; 
22// see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 
23// <http://www.gnu.org/licenses/>. 
24 
25/** @file tr1/exp_integral.tcc 
26 * This is an internal header file, included by other library headers. 
27 * Do not attempt to use it directly. @headername{tr1/cmath} 
28 */ 
29 
30// 
31// ISO C++ 14882 TR1: 5.2 Special functions 
32// 
33 
34// Written by Edward Smith-Rowland based on: 
35// 
36// (1) Handbook of Mathematical Functions, 
37// Ed. by Milton Abramowitz and Irene A. Stegun, 
38// Dover Publications, New-York, Section 5, pp. 228-251. 
39// (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl 
40// (3) Numerical Recipes in C, by W. H. Press, S. A. Teukolsky, 
41// W. T. Vetterling, B. P. Flannery, Cambridge University Press (1992), 
42// 2nd ed, pp. 222-225. 
43// 
44 
45#ifndef _GLIBCXX_TR1_EXP_INTEGRAL_TCC 
46#define _GLIBCXX_TR1_EXP_INTEGRAL_TCC 1 
47 
48#include <tr1/special_function_util.h> 
49 
50namespace std _GLIBCXX_VISIBILITY(default) 
51{ 
52_GLIBCXX_BEGIN_NAMESPACE_VERSION 
53 
54#if _GLIBCXX_USE_STD_SPEC_FUNCS 
55#elif defined(_GLIBCXX_TR1_CMATH) 
56namespace tr1 
57{ 
58#else 
59# error do not include this header directly, use <cmath> or <tr1/cmath> 
60#endif 
61 // [5.2] Special functions 
62 
63 // Implementation-space details. 
64 namespace __detail 
65 { 
66 template<typename _Tp> _Tp __expint_E1(_Tp); 
67 
68 /** 
69 * @brief Return the exponential integral @f$ E_1(x) @f$ 
70 * by series summation. This should be good 
71 * for @f$ x < 1 @f$. 
72 *  
73 * The exponential integral is given by 
74 * \f[ 
75 * E_1(x) = \int_{1}^{\infty} \frac{e^{-xt}}{t} dt 
76 * \f] 
77 *  
78 * @param __x The argument of the exponential integral function. 
79 * @return The exponential integral. 
80 */ 
81 template<typename _Tp> 
82 _Tp 
83 __expint_E1_series(_Tp __x) 
84 { 
85 const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 
86 _Tp __term = _Tp(1); 
87 _Tp __esum = _Tp(0); 
88 _Tp __osum = _Tp(0); 
89 const unsigned int __max_iter = 1000; 
90 for (unsigned int __i = 1; __i < __max_iter; ++__i) 
91 { 
92 __term *= - __x / __i; 
93 if (std::abs(__term) < __eps) 
94 break; 
95 if (__term >= _Tp(0)) 
96 __esum += __term / __i; 
97 else 
98 __osum += __term / __i; 
99 } 
100 
101 return - __esum - __osum 
102 - __numeric_constants<_Tp>::__gamma_e() - std::log(__x); 
103 } 
104 
105 
106 /** 
107 * @brief Return the exponential integral @f$ E_1(x) @f$ 
108 * by asymptotic expansion. 
109 *  
110 * The exponential integral is given by 
111 * \f[ 
112 * E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt 
113 * \f] 
114 *  
115 * @param __x The argument of the exponential integral function. 
116 * @return The exponential integral. 
117 */ 
118 template<typename _Tp> 
119 _Tp 
120 __expint_E1_asymp(_Tp __x) 
121 { 
122 _Tp __term = _Tp(1); 
123 _Tp __esum = _Tp(1); 
124 _Tp __osum = _Tp(0); 
125 const unsigned int __max_iter = 1000; 
126 for (unsigned int __i = 1; __i < __max_iter; ++__i) 
127 { 
128 _Tp __prev = __term; 
129 __term *= - __i / __x; 
130 if (std::abs(__term) > std::abs(__prev)) 
131 break; 
132 if (__term >= _Tp(0)) 
133 __esum += __term; 
134 else 
135 __osum += __term; 
136 } 
137 
138 return std::exp(- __x) * (__esum + __osum) / __x; 
139 } 
140 
141 
142 /** 
143 * @brief Return the exponential integral @f$ E_n(x) @f$ 
144 * by series summation. 
145 *  
146 * The exponential integral is given by 
147 * \f[ 
148 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 
149 * \f] 
150 *  
151 * @param __n The order of the exponential integral function. 
152 * @param __x The argument of the exponential integral function. 
153 * @return The exponential integral. 
154 */ 
155 template<typename _Tp> 
156 _Tp 
157 __expint_En_series(unsigned int __n, _Tp __x) 
158 { 
159 const unsigned int __max_iter = 1000; 
160 const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 
161 const int __nm1 = __n - 1; 
162 _Tp __ans = (__nm1 != 0 
163 ? _Tp(1) / __nm1 : -std::log(__x) 
164 - __numeric_constants<_Tp>::__gamma_e()); 
165 _Tp __fact = _Tp(1); 
166 for (int __i = 1; __i <= __max_iter; ++__i) 
167 { 
168 __fact *= -__x / _Tp(__i); 
169 _Tp __del; 
170 if ( __i != __nm1 ) 
171 __del = -__fact / _Tp(__i - __nm1); 
172 else 
173 { 
174 _Tp __psi = -__numeric_constants<_Tp>::gamma_e(); 
175 for (int __ii = 1; __ii <= __nm1; ++__ii) 
176 __psi += _Tp(1) / _Tp(__ii); 
177 __del = __fact * (__psi - std::log(__x));  
178 } 
179 __ans += __del; 
180 if (std::abs(__del) < __eps * std::abs(__ans)) 
181 return __ans; 
182 } 
183 std::__throw_runtime_error(__N("Series summation failed " 
184 "in __expint_En_series.")); 
185 } 
186 
187 
188 /** 
189 * @brief Return the exponential integral @f$ E_n(x) @f$ 
190 * by continued fractions. 
191 *  
192 * The exponential integral is given by 
193 * \f[ 
194 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 
195 * \f] 
196 *  
197 * @param __n The order of the exponential integral function. 
198 * @param __x The argument of the exponential integral function. 
199 * @return The exponential integral. 
200 */ 
201 template<typename _Tp> 
202 _Tp 
203 __expint_En_cont_frac(unsigned int __n, _Tp __x) 
204 { 
205 const unsigned int __max_iter = 1000; 
206 const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 
207 const _Tp __fp_min = std::numeric_limits<_Tp>::min(); 
208 const int __nm1 = __n - 1; 
209 _Tp __b = __x + _Tp(__n); 
210 _Tp __c = _Tp(1) / __fp_min; 
211 _Tp __d = _Tp(1) / __b; 
212 _Tp __h = __d; 
213 for ( unsigned int __i = 1; __i <= __max_iter; ++__i ) 
214 { 
215 _Tp __a = -_Tp(__i * (__nm1 + __i)); 
216 __b += _Tp(2); 
217 __d = _Tp(1) / (__a * __d + __b); 
218 __c = __b + __a / __c; 
219 const _Tp __del = __c * __d; 
220 __h *= __del; 
221 if (std::abs(__del - _Tp(1)) < __eps) 
222 { 
223 const _Tp __ans = __h * std::exp(-__x); 
224 return __ans; 
225 } 
226 } 
227 std::__throw_runtime_error(__N("Continued fraction failed " 
228 "in __expint_En_cont_frac.")); 
229 } 
230 
231 
232 /** 
233 * @brief Return the exponential integral @f$ E_n(x) @f$ 
234 * by recursion. Use upward recursion for @f$ x < n @f$ 
235 * and downward recursion (Miller's algorithm) otherwise. 
236 *  
237 * The exponential integral is given by 
238 * \f[ 
239 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 
240 * \f] 
241 *  
242 * @param __n The order of the exponential integral function. 
243 * @param __x The argument of the exponential integral function. 
244 * @return The exponential integral. 
245 */ 
246 template<typename _Tp> 
247 _Tp 
248 __expint_En_recursion(unsigned int __n, _Tp __x) 
249 { 
250 _Tp __En; 
251 _Tp __E1 = __expint_E1(__x); 
252 if (__x < _Tp(__n)) 
253 { 
254 // Forward recursion is stable only for n < x. 
255 __En = __E1; 
256 for (unsigned int __j = 2; __j < __n; ++__j) 
257 __En = (std::exp(-__x) - __x * __En) / _Tp(__j - 1); 
258 } 
259 else 
260 { 
261 // Backward recursion is stable only for n >= x. 
262 __En = _Tp(1); 
263 const int __N = __n + 20; // TODO: Check this starting number. 
264 _Tp __save = _Tp(0); 
265 for (int __j = __N; __j > 0; --__j) 
266 { 
267 __En = (std::exp(-__x) - __j * __En) / __x; 
268 if (__j == __n) 
269 __save = __En; 
270 } 
271 _Tp __norm = __En / __E1; 
272 __En /= __norm; 
273 } 
274 
275 return __En; 
276 } 
277 
278 /** 
279 * @brief Return the exponential integral @f$ Ei(x) @f$ 
280 * by series summation. 
281 *  
282 * The exponential integral is given by 
283 * \f[ 
284 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt 
285 * \f] 
286 *  
287 * @param __x The argument of the exponential integral function. 
288 * @return The exponential integral. 
289 */ 
290 template<typename _Tp> 
291 _Tp 
292 __expint_Ei_series(_Tp __x) 
293 { 
294 _Tp __term = _Tp(1); 
295 _Tp __sum = _Tp(0); 
296 const unsigned int __max_iter = 1000; 
297 for (unsigned int __i = 1; __i < __max_iter; ++__i) 
298 { 
299 __term *= __x / __i; 
300 __sum += __term / __i; 
301 if (__term < std::numeric_limits<_Tp>::epsilon() * __sum) 
302 break; 
303 } 
304 
305 return __numeric_constants<_Tp>::__gamma_e() + __sum + std::log(__x); 
306 } 
307 
308 
309 /** 
310 * @brief Return the exponential integral @f$ Ei(x) @f$ 
311 * by asymptotic expansion. 
312 *  
313 * The exponential integral is given by 
314 * \f[ 
315 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt 
316 * \f] 
317 *  
318 * @param __x The argument of the exponential integral function. 
319 * @return The exponential integral. 
320 */ 
321 template<typename _Tp> 
322 _Tp 
323 __expint_Ei_asymp(_Tp __x) 
324 { 
325 _Tp __term = _Tp(1); 
326 _Tp __sum = _Tp(1); 
327 const unsigned int __max_iter = 1000; 
328 for (unsigned int __i = 1; __i < __max_iter; ++__i) 
329 { 
330 _Tp __prev = __term; 
331 __term *= __i / __x; 
332 if (__term < std::numeric_limits<_Tp>::epsilon()) 
333 break; 
334 if (__term >= __prev) 
335 break; 
336 __sum += __term; 
337 } 
338 
339 return std::exp(__x) * __sum / __x; 
340 } 
341 
342 
343 /** 
344 * @brief Return the exponential integral @f$ Ei(x) @f$. 
345 *  
346 * The exponential integral is given by 
347 * \f[ 
348 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt 
349 * \f] 
350 *  
351 * @param __x The argument of the exponential integral function. 
352 * @return The exponential integral. 
353 */ 
354 template<typename _Tp> 
355 _Tp 
356 __expint_Ei(_Tp __x) 
357 { 
358 if (__x < _Tp(0)) 
359 return -__expint_E1(-__x); 
360 else if (__x < -std::log(std::numeric_limits<_Tp>::epsilon())) 
361 return __expint_Ei_series(__x); 
362 else 
363 return __expint_Ei_asymp(__x); 
364 } 
365 
366 
367 /** 
368 * @brief Return the exponential integral @f$ E_1(x) @f$. 
369 *  
370 * The exponential integral is given by 
371 * \f[ 
372 * E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt 
373 * \f] 
374 *  
375 * @param __x The argument of the exponential integral function. 
376 * @return The exponential integral. 
377 */ 
378 template<typename _Tp> 
379 _Tp 
380 __expint_E1(_Tp __x) 
381 { 
382 if (__x < _Tp(0)) 
383 return -__expint_Ei(-__x); 
384 else if (__x < _Tp(1)) 
385 return __expint_E1_series(__x); 
386 else if (__x < _Tp(100)) // TODO: Find a good asymptotic switch point. 
387 return __expint_En_cont_frac(1, __x); 
388 else 
389 return __expint_E1_asymp(__x); 
390 } 
391 
392 
393 /** 
394 * @brief Return the exponential integral @f$ E_n(x) @f$ 
395 * for large argument. 
396 *  
397 * The exponential integral is given by 
398 * \f[ 
399 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 
400 * \f] 
401 *  
402 * This is something of an extension. 
403 *  
404 * @param __n The order of the exponential integral function. 
405 * @param __x The argument of the exponential integral function. 
406 * @return The exponential integral. 
407 */ 
408 template<typename _Tp> 
409 _Tp 
410 __expint_asymp(unsigned int __n, _Tp __x) 
411 { 
412 _Tp __term = _Tp(1); 
413 _Tp __sum = _Tp(1); 
414 for (unsigned int __i = 1; __i <= __n; ++__i) 
415 { 
416 _Tp __prev = __term; 
417 __term *= -(__n - __i + 1) / __x; 
418 if (std::abs(__term) > std::abs(__prev)) 
419 break; 
420 __sum += __term; 
421 } 
422 
423 return std::exp(-__x) * __sum / __x; 
424 } 
425 
426 
427 /** 
428 * @brief Return the exponential integral @f$ E_n(x) @f$ 
429 * for large order. 
430 *  
431 * The exponential integral is given by 
432 * \f[ 
433 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 
434 * \f] 
435 *  
436 * This is something of an extension. 
437 *  
438 * @param __n The order of the exponential integral function. 
439 * @param __x The argument of the exponential integral function. 
440 * @return The exponential integral. 
441 */ 
442 template<typename _Tp> 
443 _Tp 
444 __expint_large_n(unsigned int __n, _Tp __x) 
445 { 
446 const _Tp __xpn = __x + __n; 
447 const _Tp __xpn2 = __xpn * __xpn; 
448 _Tp __term = _Tp(1); 
449 _Tp __sum = _Tp(1); 
450 for (unsigned int __i = 1; __i <= __n; ++__i) 
451 { 
452 _Tp __prev = __term; 
453 __term *= (__n - 2 * (__i - 1) * __x) / __xpn2; 
454 if (std::abs(__term) < std::numeric_limits<_Tp>::epsilon()) 
455 break; 
456 __sum += __term; 
457 } 
458 
459 return std::exp(-__x) * __sum / __xpn; 
460 } 
461 
462 
463 /** 
464 * @brief Return the exponential integral @f$ E_n(x) @f$. 
465 *  
466 * The exponential integral is given by 
467 * \f[ 
468 * E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt 
469 * \f] 
470 * This is something of an extension. 
471 *  
472 * @param __n The order of the exponential integral function. 
473 * @param __x The argument of the exponential integral function. 
474 * @return The exponential integral. 
475 */ 
476 template<typename _Tp> 
477 _Tp 
478 __expint(unsigned int __n, _Tp __x) 
479 { 
480 // Return NaN on NaN input. 
481 if (__isnan(__x)) 
482 return std::numeric_limits<_Tp>::quiet_NaN(); 
483 else if (__n <= 1 && __x == _Tp(0)) 
484 return std::numeric_limits<_Tp>::infinity(); 
485 else 
486 { 
487 _Tp __E0 = std::exp(__x) / __x; 
488 if (__n == 0) 
489 return __E0; 
490 
491 _Tp __E1 = __expint_E1(__x); 
492 if (__n == 1) 
493 return __E1; 
494 
495 if (__x == _Tp(0)) 
496 return _Tp(1) / static_cast<_Tp>(__n - 1); 
497 
498 _Tp __En = __expint_En_recursion(__n, __x); 
499 
500 return __En; 
501 } 
502 } 
503 
504 
505 /** 
506 * @brief Return the exponential integral @f$ Ei(x) @f$. 
507 *  
508 * The exponential integral is given by 
509 * \f[ 
510 * Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt 
511 * \f] 
512 *  
513 * @param __x The argument of the exponential integral function. 
514 * @return The exponential integral. 
515 */ 
516 template<typename _Tp> 
517 inline _Tp 
518 __expint(_Tp __x) 
519 { 
520 if (__isnan(__x)) 
521 return std::numeric_limits<_Tp>::quiet_NaN(); 
522 else 
523 return __expint_Ei(__x); 
524 } 
525 } // namespace __detail 
526#if ! _GLIBCXX_USE_STD_SPEC_FUNCS && defined(_GLIBCXX_TR1_CMATH) 
527} // namespace tr1 
528#endif 
529 
530_GLIBCXX_END_NAMESPACE_VERSION 
531} 
532 
533#endif // _GLIBCXX_TR1_EXP_INTEGRAL_TCC 
534