1 // Special functions -*- C++ -*-
2 
3 // Copyright (C) 2006-2020 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 
50 namespace std _GLIBCXX_VISIBILITY(default)
51 {
52 _GLIBCXX_BEGIN_NAMESPACE_VERSION
53 
54 #if _GLIBCXX_USE_STD_SPEC_FUNCS
55 #elif defined(_GLIBCXX_TR1_CMATH)
56 namespace 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
__expint_E1_series(_Tp __x)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
__expint_E1_asymp(_Tp __x)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
__expint_En_series(unsigned int __n,_Tp __x)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
__expint_En_cont_frac(unsigned int __n,_Tp __x)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
__expint_En_recursion(unsigned int __n,_Tp __x)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
__expint_Ei_series(_Tp __x)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
__expint_Ei_asymp(_Tp __x)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
__expint_Ei(_Tp __x)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
__expint_E1(_Tp __x)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
__expint_asymp(unsigned int __n,_Tp __x)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
__expint_large_n(unsigned int __n,_Tp __x)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
__expint(unsigned int __n,_Tp __x)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
__expint(_Tp __x)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