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