1// Special functions -*- C++ -*-
2
3// Copyright (C) 2006-2021 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