STLdoc
STLdocumentation
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Macros | Functions
exp_integral.tcc File Reference
#include "special_function_util.h"

Macros

#define _GLIBCXX_TR1_EXP_INTEGRAL_TCC   1
 

Functions

namespace std _GLIBCXX_VISIBILITY (default)
 

Detailed Description

This is an internal header file, included by other library headers. Do not attempt to use it directly. {tr1/cmath}

Macro Definition Documentation

#define _GLIBCXX_TR1_EXP_INTEGRAL_TCC   1

Function Documentation

namespace std _GLIBCXX_VISIBILITY ( default  )

Return the exponential integral $ E_1(x) $ by series summation. This should be good for $ x < 1 $.

The exponential integral is given by

\[ E_1(x) = \int_{1}^{\infty} \frac{e^{-xt}}{t} dt \]

Parameters
__xThe argument of the exponential integral function.
Returns
The exponential integral.

Return the exponential integral $ E_1(x) $ by asymptotic expansion.

The exponential integral is given by

\[ E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt \]

Parameters
__xThe argument of the exponential integral function.
Returns
The exponential integral.

Return the exponential integral $ E_n(x) $ by series summation.

The exponential integral is given by

\[ E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt \]

Parameters
__nThe order of the exponential integral function.
__xThe argument of the exponential integral function.
Returns
The exponential integral.

Return the exponential integral $ E_n(x) $ by continued fractions.

The exponential integral is given by

\[ E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt \]

Parameters
__nThe order of the exponential integral function.
__xThe argument of the exponential integral function.
Returns
The exponential integral.

Return the exponential integral $ E_n(x) $ by recursion. Use upward recursion for $ x < n $ and downward recursion (Miller's algorithm) otherwise.

The exponential integral is given by

\[ E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt \]

Parameters
__nThe order of the exponential integral function.
__xThe argument of the exponential integral function.
Returns
The exponential integral.

Return the exponential integral $ Ei(x) $ by series summation.

The exponential integral is given by

\[ Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt \]

Parameters
__xThe argument of the exponential integral function.
Returns
The exponential integral.

Return the exponential integral $ Ei(x) $ by asymptotic expansion.

The exponential integral is given by

\[ Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt \]

Parameters
__xThe argument of the exponential integral function.
Returns
The exponential integral.

Return the exponential integral $ Ei(x) $.

The exponential integral is given by

\[ Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt \]

Parameters
__xThe argument of the exponential integral function.
Returns
The exponential integral.

Return the exponential integral $ E_1(x) $.

The exponential integral is given by

\[ E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt \]

Parameters
__xThe argument of the exponential integral function.
Returns
The exponential integral.

Return the exponential integral $ E_n(x) $ for large argument.

The exponential integral is given by

\[ E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt \]

This is something of an extension.

Parameters
__nThe order of the exponential integral function.
__xThe argument of the exponential integral function.
Returns
The exponential integral.

Return the exponential integral $ E_n(x) $ for large order.

The exponential integral is given by

\[ E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt \]

This is something of an extension.

Parameters
__nThe order of the exponential integral function.
__xThe argument of the exponential integral function.
Returns
The exponential integral.

Return the exponential integral $ E_n(x) $.

The exponential integral is given by

\[ E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt \]

This is something of an extension.

Parameters
__nThe order of the exponential integral function.
__xThe argument of the exponential integral function.
Returns
The exponential integral.

Return the exponential integral $ Ei(x) $.

The exponential integral is given by

\[ Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt \]

Parameters
__xThe argument of the exponential integral function.
Returns
The exponential integral.
51 {
52 namespace tr1
53 {
54  // [5.2] Special functions
55 
56  // Implementation-space details.
57  namespace __detail
58  {
59  _GLIBCXX_BEGIN_NAMESPACE_VERSION
60 
61  template<typename _Tp> _Tp __expint_E1(_Tp);
62 
76  template<typename _Tp>
77  _Tp
78  __expint_E1_series(_Tp __x)
79  {
80  const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
81  _Tp __term = _Tp(1);
82  _Tp __esum = _Tp(0);
83  _Tp __osum = _Tp(0);
84  const unsigned int __max_iter = 100;
85  for (unsigned int __i = 1; __i < __max_iter; ++__i)
86  {
87  __term *= - __x / __i;
88  if (std::abs(__term) < __eps)
89  break;
90  if (__term >= _Tp(0))
91  __esum += __term / __i;
92  else
93  __osum += __term / __i;
94  }
95 
96  return - __esum - __osum
97  - __numeric_constants<_Tp>::__gamma_e() - std::log(__x);
98  }
99 
100 
113  template<typename _Tp>
114  _Tp
115  __expint_E1_asymp(_Tp __x)
116  {
117  _Tp __term = _Tp(1);
118  _Tp __esum = _Tp(1);
119  _Tp __osum = _Tp(0);
120  const unsigned int __max_iter = 1000;
121  for (unsigned int __i = 1; __i < __max_iter; ++__i)
122  {
123  _Tp __prev = __term;
124  __term *= - __i / __x;
125  if (std::abs(__term) > std::abs(__prev))
126  break;
127  if (__term >= _Tp(0))
128  __esum += __term;
129  else
130  __osum += __term;
131  }
132 
133  return std::exp(- __x) * (__esum + __osum) / __x;
134  }
135 
136 
150  template<typename _Tp>
151  _Tp
152  __expint_En_series(unsigned int __n, _Tp __x)
153  {
154  const unsigned int __max_iter = 100;
155  const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
156  const int __nm1 = __n - 1;
157  _Tp __ans = (__nm1 != 0
158  ? _Tp(1) / __nm1 : -std::log(__x)
159  - __numeric_constants<_Tp>::__gamma_e());
160  _Tp __fact = _Tp(1);
161  for (int __i = 1; __i <= __max_iter; ++__i)
162  {
163  __fact *= -__x / _Tp(__i);
164  _Tp __del;
165  if ( __i != __nm1 )
166  __del = -__fact / _Tp(__i - __nm1);
167  else
168  {
169  _Tp __psi = -__numeric_constants<_Tp>::gamma_e();
170  for (int __ii = 1; __ii <= __nm1; ++__ii)
171  __psi += _Tp(1) / _Tp(__ii);
172  __del = __fact * (__psi - std::log(__x));
173  }
174  __ans += __del;
175  if (std::abs(__del) < __eps * std::abs(__ans))
176  return __ans;
177  }
178  std::__throw_runtime_error(__N("Series summation failed "
179  "in __expint_En_series."));
180  }
181 
182 
196  template<typename _Tp>
197  _Tp
198  __expint_En_cont_frac(unsigned int __n, _Tp __x)
199  {
200  const unsigned int __max_iter = 100;
201  const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
202  const _Tp __fp_min = std::numeric_limits<_Tp>::min();
203  const int __nm1 = __n - 1;
204  _Tp __b = __x + _Tp(__n);
205  _Tp __c = _Tp(1) / __fp_min;
206  _Tp __d = _Tp(1) / __b;
207  _Tp __h = __d;
208  for ( unsigned int __i = 1; __i <= __max_iter; ++__i )
209  {
210  _Tp __a = -_Tp(__i * (__nm1 + __i));
211  __b += _Tp(2);
212  __d = _Tp(1) / (__a * __d + __b);
213  __c = __b + __a / __c;
214  const _Tp __del = __c * __d;
215  __h *= __del;
216  if (std::abs(__del - _Tp(1)) < __eps)
217  {
218  const _Tp __ans = __h * std::exp(-__x);
219  return __ans;
220  }
221  }
222  std::__throw_runtime_error(__N("Continued fraction failed "
223  "in __expint_En_cont_frac."));
224  }
225 
226 
241  template<typename _Tp>
242  _Tp
243  __expint_En_recursion(unsigned int __n, _Tp __x)
244  {
245  _Tp __En;
246  _Tp __E1 = __expint_E1(__x);
247  if (__x < _Tp(__n))
248  {
249  // Forward recursion is stable only for n < x.
250  __En = __E1;
251  for (unsigned int __j = 2; __j < __n; ++__j)
252  __En = (std::exp(-__x) - __x * __En) / _Tp(__j - 1);
253  }
254  else
255  {
256  // Backward recursion is stable only for n >= x.
257  __En = _Tp(1);
258  const int __N = __n + 20; // TODO: Check this starting number.
259  _Tp __save = _Tp(0);
260  for (int __j = __N; __j > 0; --__j)
261  {
262  __En = (std::exp(-__x) - __j * __En) / __x;
263  if (__j == __n)
264  __save = __En;
265  }
266  _Tp __norm = __En / __E1;
267  __En /= __norm;
268  }
269 
270  return __En;
271  }
272 
285  template<typename _Tp>
286  _Tp
287  __expint_Ei_series(_Tp __x)
288  {
289  _Tp __term = _Tp(1);
290  _Tp __sum = _Tp(0);
291  const unsigned int __max_iter = 1000;
292  for (unsigned int __i = 1; __i < __max_iter; ++__i)
293  {
294  __term *= __x / __i;
295  __sum += __term / __i;
296  if (__term < std::numeric_limits<_Tp>::epsilon() * __sum)
297  break;
298  }
299 
300  return __numeric_constants<_Tp>::__gamma_e() + __sum + std::log(__x);
301  }
302 
303 
316  template<typename _Tp>
317  _Tp
318  __expint_Ei_asymp(_Tp __x)
319  {
320  _Tp __term = _Tp(1);
321  _Tp __sum = _Tp(1);
322  const unsigned int __max_iter = 1000;
323  for (unsigned int __i = 1; __i < __max_iter; ++__i)
324  {
325  _Tp __prev = __term;
326  __term *= __i / __x;
327  if (__term < std::numeric_limits<_Tp>::epsilon())
328  break;
329  if (__term >= __prev)
330  break;
331  __sum += __term;
332  }
333 
334  return std::exp(__x) * __sum / __x;
335  }
336 
337 
349  template<typename _Tp>
350  _Tp
351  __expint_Ei(_Tp __x)
352  {
353  if (__x < _Tp(0))
354  return -__expint_E1(-__x);
355  else if (__x < -std::log(std::numeric_limits<_Tp>::epsilon()))
356  return __expint_Ei_series(__x);
357  else
358  return __expint_Ei_asymp(__x);
359  }
360 
361 
373  template<typename _Tp>
374  _Tp
375  __expint_E1(_Tp __x)
376  {
377  if (__x < _Tp(0))
378  return -__expint_Ei(-__x);
379  else if (__x < _Tp(1))
380  return __expint_E1_series(__x);
381  else if (__x < _Tp(100)) // TODO: Find a good asymptotic switch point.
382  return __expint_En_cont_frac(1, __x);
383  else
384  return __expint_E1_asymp(__x);
385  }
386 
387 
403  template<typename _Tp>
404  _Tp
405  __expint_asymp(unsigned int __n, _Tp __x)
406  {
407  _Tp __term = _Tp(1);
408  _Tp __sum = _Tp(1);
409  for (unsigned int __i = 1; __i <= __n; ++__i)
410  {
411  _Tp __prev = __term;
412  __term *= -(__n - __i + 1) / __x;
413  if (std::abs(__term) > std::abs(__prev))
414  break;
415  __sum += __term;
416  }
417 
418  return std::exp(-__x) * __sum / __x;
419  }
420 
421 
437  template<typename _Tp>
438  _Tp
439  __expint_large_n(unsigned int __n, _Tp __x)
440  {
441  const _Tp __xpn = __x + __n;
442  const _Tp __xpn2 = __xpn * __xpn;
443  _Tp __term = _Tp(1);
444  _Tp __sum = _Tp(1);
445  for (unsigned int __i = 1; __i <= __n; ++__i)
446  {
447  _Tp __prev = __term;
448  __term *= (__n - 2 * (__i - 1) * __x) / __xpn2;
449  if (std::abs(__term) < std::numeric_limits<_Tp>::epsilon())
450  break;
451  __sum += __term;
452  }
453 
454  return std::exp(-__x) * __sum / __xpn;
455  }
456 
457 
471  template<typename _Tp>
472  _Tp
473  __expint(unsigned int __n, _Tp __x)
474  {
475  // Return NaN on NaN input.
476  if (__isnan(__x))
477  return std::numeric_limits<_Tp>::quiet_NaN();
478  else if (__n <= 1 && __x == _Tp(0))
479  return std::numeric_limits<_Tp>::infinity();
480  else
481  {
482  _Tp __E0 = std::exp(__x) / __x;
483  if (__n == 0)
484  return __E0;
485 
486  _Tp __E1 = __expint_E1(__x);
487  if (__n == 1)
488  return __E1;
489 
490  if (__x == _Tp(0))
491  return _Tp(1) / static_cast<_Tp>(__n - 1);
492 
493  _Tp __En = __expint_En_recursion(__n, __x);
494 
495  return __En;
496  }
497  }
498 
499 
511  template<typename _Tp>
512  inline _Tp
513  __expint(_Tp __x)
514  {
515  if (__isnan(__x))
516  return std::numeric_limits<_Tp>::quiet_NaN();
517  else
518  return __expint_Ei(__x);
519  }
520 
521  _GLIBCXX_END_NAMESPACE_VERSION
522  } // namespace std::tr1::__detail
523 }
524 }
__gnu_parallel::min
const _Tp & min(const _Tp &__a, const _Tp &__b)
Equivalent to std::min.
Definition: base.h:144
Generated on Sat Jul 22 2017 14:42:36 for STLdoc by   doxygen 1.8.6