contrib/gcc-5.0/libstdc++-v3/include/tr1/exp_integral.tcc

   1 // Special functions -*- C++ -*-
   2
   3 // Copyright (C) 2006-2015 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 "special_function_util.h"
  49
  50 namespace std _GLIBCXX_VISIBILITY(default)
  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
  63     /**
  64      *   @brief Return the exponential integral @f$ E_1(x) @f$
  65      *          by series summation.  This should be good
  66      *          for @f$ x < 1 @f$.
  67      *
  68      *   The exponential integral is given by
  69      *          \f[
  70      *            E_1(x) = \int_{1}^{\infty} \frac{e^{-xt}}{t} dt
  71      *          \f]
  72      *
  73      *   @param  __x  The argument of the exponential integral function.
  74      *   @return  The exponential integral.
  75      */
  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
 101     /**
 102      *   @brief Return the exponential integral @f$ E_1(x) @f$
 103      *          by asymptotic expansion.
 104      *
 105      *   The exponential integral is given by
 106      *          \f[
 107      *            E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt
 108      *          \f]
 109      *
 110      *   @param  __x  The argument of the exponential integral function.
 111      *   @return  The exponential integral.
 112      */
 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
 137     /**
 138      *   @brief Return the exponential integral @f$ E_n(x) @f$
 139      *          by series summation.
 140      *
 141      *   The exponential integral is given by
 142      *          \f[
 143      *            E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
 144      *          \f]
 145      *
 146      *   @param  __n  The order of the exponential integral function.
 147      *   @param  __x  The argument of the exponential integral function.
 148      *   @return  The exponential integral.
 149      */
 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
 183     /**
 184      *   @brief Return the exponential integral @f$ E_n(x) @f$
 185      *          by continued fractions.
 186      *
 187      *   The exponential integral is given by
 188      *          \f[
 189      *            E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
 190      *          \f]
 191      *
 192      *   @param  __n  The order of the exponential integral function.
 193      *   @param  __x  The argument of the exponential integral function.
 194      *   @return  The exponential integral.
 195      */
 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
 227     /**
 228      *   @brief Return the exponential integral @f$ E_n(x) @f$
 229      *          by recursion.  Use upward recursion for @f$ x < n @f$
 230      *          and downward recursion (Miller's algorithm) otherwise.
 231      *
 232      *   The exponential integral is given by
 233      *          \f[
 234      *            E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
 235      *          \f]
 236      *
 237      *   @param  __n  The order of the exponential integral function.
 238      *   @param  __x  The argument of the exponential integral function.
 239      *   @return  The exponential integral.
 240      */
 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
 273     /**
 274      *   @brief Return the exponential integral @f$ Ei(x) @f$
 275      *          by series summation.
 276      *
 277      *   The exponential integral is given by
 278      *          \f[
 279      *            Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
 280      *          \f]
 281      *
 282      *   @param  __x  The argument of the exponential integral function.
 283      *   @return  The exponential integral.
 284      */
 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
 304     /**
 305      *   @brief Return the exponential integral @f$ Ei(x) @f$
 306      *          by asymptotic expansion.
 307      *
 308      *   The exponential integral is given by
 309      *          \f[
 310      *            Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
 311      *          \f]
 312      *
 313      *   @param  __x  The argument of the exponential integral function.
 314      *   @return  The exponential integral.
 315      */
 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
 338     /**
 339      *   @brief Return the exponential integral @f$ Ei(x) @f$.
 340      *
 341      *   The exponential integral is given by
 342      *          \f[
 343      *            Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
 344      *          \f]
 345      *
 346      *   @param  __x  The argument of the exponential integral function.
 347      *   @return  The exponential integral.
 348      */
 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
 362     /**
 363      *   @brief Return the exponential integral @f$ E_1(x) @f$.
 364      *
 365      *   The exponential integral is given by
 366      *          \f[
 367      *            E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt
 368      *          \f]
 369      *
 370      *   @param  __x  The argument of the exponential integral function.
 371      *   @return  The exponential integral.
 372      */
 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
 388     /**
 389      *   @brief Return the exponential integral @f$ E_n(x) @f$
 390      *          for large argument.
 391      *
 392      *   The exponential integral is given by
 393      *          \f[
 394      *            E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
 395      *          \f]
 396      *
 397      *   This is something of an extension.
 398      *
 399      *   @param  __n  The order of the exponential integral function.
 400      *   @param  __x  The argument of the exponential integral function.
 401      *   @return  The exponential integral.
 402      */
 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
 422     /**
 423      *   @brief Return the exponential integral @f$ E_n(x) @f$
 424      *          for large order.
 425      *
 426      *   The exponential integral is given by
 427      *          \f[
 428      *            E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
 429      *          \f]
 430      *
 431      *   This is something of an extension.
 432      *
 433      *   @param  __n  The order of the exponential integral function.
 434      *   @param  __x  The argument of the exponential integral function.
 435      *   @return  The exponential integral.
 436      */
 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
 458     /**
 459      *   @brief Return the exponential integral @f$ E_n(x) @f$.
 460      *
 461      *   The exponential integral is given by
 462      *          \f[
 463      *            E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
 464      *          \f]
 465      *   This is something of an extension.
 466      *
 467      *   @param  __n  The order of the exponential integral function.
 468      *   @param  __x  The argument of the exponential integral function.
 469      *   @return  The exponential integral.
 470      */
 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
 500     /**
 501      *   @brief Return the exponential integral @f$ Ei(x) @f$.
 502      *
 503      *   The exponential integral is given by
 504      *   \f[
 505      *     Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
 506      *   \f]
 507      *
 508      *   @param  __x  The argument of the exponential integral function.
 509      *   @return  The exponential integral.
 510      */
 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 }
 525
 526 #endif // _GLIBCXX_TR1_EXP_INTEGRAL_TCC