i686-linux-gnu-4.7/usr/include/c++/4.7/tr1/exp_integral.tcc

   1 // Special functions -*- C++ -*-
   2
   3 // Copyright (C) 2006, 2007, 2008, 2009, 2010
   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  *  Do not attempt to use it directly. @headername{tr1/cmath}
  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 _GLIBCXX_VISIBILITY(default)
  52 {
  53 namespace tr1
  54 {
  55   // [5.2] Special functions
  56
  57   // Implementation-space details.
  58   namespace __detail
  59   {
  60   _GLIBCXX_BEGIN_NAMESPACE_VERSION
  61
  62     template<typename _Tp> _Tp __expint_E1(const _Tp);
  63
  64     /**
  65      *   @brief Return the exponential integral @f$ E_1(x) @f$
  66      *          by series summation.  This should be good
  67      *          for @f$ x < 1 @f$.
  68      *
  69      *   The exponential integral is given by
  70      *          \f[
  71      *            E_1(x) = \int_{1}^{\infty} \frac{e^{-xt}}{t} dt
  72      *          \f]
  73      *
  74      *   @param  __x  The argument of the exponential integral function.
  75      *   @return  The exponential integral.
  76      */
  77     template<typename _Tp>
  78     _Tp
  79     __expint_E1_series(const _Tp __x)
  80     {
  81       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
  82       _Tp __term = _Tp(1);
  83       _Tp __esum = _Tp(0);
  84       _Tp __osum = _Tp(0);
  85       const unsigned int __max_iter = 100;
  86       for (unsigned int __i = 1; __i < __max_iter; ++__i)
  87         {
  88           __term *= - __x / __i;
  89           if (std::abs(__term) < __eps)
  90             break;
  91           if (__term >= _Tp(0))
  92             __esum += __term / __i;
  93           else
  94             __osum += __term / __i;
  95         }
  96
  97       return - __esum - __osum
  98              - __numeric_constants<_Tp>::__gamma_e() - std::log(__x);
  99     }
 100
 101
 102     /**
 103      *   @brief Return the exponential integral @f$ E_1(x) @f$
 104      *          by asymptotic expansion.
 105      *
 106      *   The exponential integral is given by
 107      *          \f[
 108      *            E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt
 109      *          \f]
 110      *
 111      *   @param  __x  The argument of the exponential integral function.
 112      *   @return  The exponential integral.
 113      */
 114     template<typename _Tp>
 115     _Tp
 116     __expint_E1_asymp(const _Tp __x)
 117     {
 118       _Tp __term = _Tp(1);
 119       _Tp __esum = _Tp(1);
 120       _Tp __osum = _Tp(0);
 121       const unsigned int __max_iter = 1000;
 122       for (unsigned int __i = 1; __i < __max_iter; ++__i)
 123         {
 124           _Tp __prev = __term;
 125           __term *= - __i / __x;
 126           if (std::abs(__term) > std::abs(__prev))
 127             break;
 128           if (__term >= _Tp(0))
 129             __esum += __term;
 130           else
 131             __osum += __term;
 132         }
 133
 134       return std::exp(- __x) * (__esum + __osum) / __x;
 135     }
 136
 137
 138     /**
 139      *   @brief Return the exponential integral @f$ E_n(x) @f$
 140      *          by series summation.
 141      *
 142      *   The exponential integral is given by
 143      *          \f[
 144      *            E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
 145      *          \f]
 146      *
 147      *   @param  __n  The order of the exponential integral function.
 148      *   @param  __x  The argument of the exponential integral function.
 149      *   @return  The exponential integral.
 150      */
 151     template<typename _Tp>
 152     _Tp
 153     __expint_En_series(const unsigned int __n, const _Tp __x)
 154     {
 155       const unsigned int __max_iter = 100;
 156       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
 157       const int __nm1 = __n - 1;
 158       _Tp __ans = (__nm1 != 0
 159                 ? _Tp(1) / __nm1 : -std::log(__x)
 160                                    - __numeric_constants<_Tp>::__gamma_e());
 161       _Tp __fact = _Tp(1);
 162       for (int __i = 1; __i <= __max_iter; ++__i)
 163         {
 164           __fact *= -__x / _Tp(__i);
 165           _Tp __del;
 166           if ( __i != __nm1 )
 167             __del = -__fact / _Tp(__i - __nm1);
 168           else
 169             {
 170               _Tp __psi = -__numeric_constants<_Tp>::gamma_e();
 171               for (int __ii = 1; __ii <= __nm1; ++__ii)
 172                 __psi += _Tp(1) / _Tp(__ii);
 173               __del = __fact * (__psi - std::log(__x));
 174             }
 175           __ans += __del;
 176           if (std::abs(__del) < __eps * std::abs(__ans))
 177             return __ans;
 178         }
 179       std::__throw_runtime_error(__N("Series summation failed "
 180                                      "in __expint_En_series."));
 181     }
 182
 183
 184     /**
 185      *   @brief Return the exponential integral @f$ E_n(x) @f$
 186      *          by continued fractions.
 187      *
 188      *   The exponential integral is given by
 189      *          \f[
 190      *            E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
 191      *          \f]
 192      *
 193      *   @param  __n  The order of the exponential integral function.
 194      *   @param  __x  The argument of the exponential integral function.
 195      *   @return  The exponential integral.
 196      */
 197     template<typename _Tp>
 198     _Tp
 199     __expint_En_cont_frac(const unsigned int __n, const _Tp __x)
 200     {
 201       const unsigned int __max_iter = 100;
 202       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
 203       const _Tp __fp_min = std::numeric_limits<_Tp>::min();
 204       const int __nm1 = __n - 1;
 205       _Tp __b = __x + _Tp(__n);
 206       _Tp __c = _Tp(1) / __fp_min;
 207       _Tp __d = _Tp(1) / __b;
 208       _Tp __h = __d;
 209       for ( unsigned int __i = 1; __i <= __max_iter; ++__i )
 210         {
 211           _Tp __a = -_Tp(__i * (__nm1 + __i));
 212           __b += _Tp(2);
 213           __d = _Tp(1) / (__a * __d + __b);
 214           __c = __b + __a / __c;
 215           const _Tp __del = __c * __d;
 216           __h *= __del;
 217           if (std::abs(__del - _Tp(1)) < __eps)
 218             {
 219               const _Tp __ans = __h * std::exp(-__x);
 220               return __ans;
 221             }
 222         }
 223       std::__throw_runtime_error(__N("Continued fraction failed "
 224                                      "in __expint_En_cont_frac."));
 225     }
 226
 227
 228     /**
 229      *   @brief Return the exponential integral @f$ E_n(x) @f$
 230      *          by recursion.  Use upward recursion for @f$ x < n @f$
 231      *          and downward recursion (Miller's algorithm) otherwise.
 232      *
 233      *   The exponential integral is given by
 234      *          \f[
 235      *            E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
 236      *          \f]
 237      *
 238      *   @param  __n  The order of the exponential integral function.
 239      *   @param  __x  The argument of the exponential integral function.
 240      *   @return  The exponential integral.
 241      */
 242     template<typename _Tp>
 243     _Tp
 244     __expint_En_recursion(const unsigned int __n, const _Tp __x)
 245     {
 246       _Tp __En;
 247       _Tp __E1 = __expint_E1(__x);
 248       if (__x < _Tp(__n))
 249         {
 250           //  Forward recursion is stable only for n < x.
 251           __En = __E1;
 252           for (unsigned int __j = 2; __j < __n; ++__j)
 253             __En = (std::exp(-__x) - __x * __En) / _Tp(__j - 1);
 254         }
 255       else
 256         {
 257           //  Backward recursion is stable only for n >= x.
 258           __En = _Tp(1);
 259           const int __N = __n + 20;  //  TODO: Check this starting number.
 260           _Tp __save = _Tp(0);
 261           for (int __j = __N; __j > 0; --__j)
 262             {
 263               __En = (std::exp(-__x) - __j * __En) / __x;
 264               if (__j == __n)
 265                 __save = __En;
 266             }
 267             _Tp __norm = __En / __E1;
 268             __En /= __norm;
 269         }
 270
 271       return __En;
 272     }
 273
 274     /**
 275      *   @brief Return the exponential integral @f$ Ei(x) @f$
 276      *          by series summation.
 277      *
 278      *   The exponential integral is given by
 279      *          \f[
 280      *            Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
 281      *          \f]
 282      *
 283      *   @param  __x  The argument of the exponential integral function.
 284      *   @return  The exponential integral.
 285      */
 286     template<typename _Tp>
 287     _Tp
 288     __expint_Ei_series(const _Tp __x)
 289     {
 290       _Tp __term = _Tp(1);
 291       _Tp __sum = _Tp(0);
 292       const unsigned int __max_iter = 1000;
 293       for (unsigned int __i = 1; __i < __max_iter; ++__i)
 294         {
 295           __term *= __x / __i;
 296           __sum += __term / __i;
 297           if (__term < std::numeric_limits<_Tp>::epsilon() * __sum)
 298             break;
 299         }
 300
 301       return __numeric_constants<_Tp>::__gamma_e() + __sum + std::log(__x);
 302     }
 303
 304
 305     /**
 306      *   @brief Return the exponential integral @f$ Ei(x) @f$
 307      *          by asymptotic expansion.
 308      *
 309      *   The exponential integral is given by
 310      *          \f[
 311      *            Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
 312      *          \f]
 313      *
 314      *   @param  __x  The argument of the exponential integral function.
 315      *   @return  The exponential integral.
 316      */
 317     template<typename _Tp>
 318     _Tp
 319     __expint_Ei_asymp(const _Tp __x)
 320     {
 321       _Tp __term = _Tp(1);
 322       _Tp __sum = _Tp(1);
 323       const unsigned int __max_iter = 1000;
 324       for (unsigned int __i = 1; __i < __max_iter; ++__i)
 325         {
 326           _Tp __prev = __term;
 327           __term *= __i / __x;
 328           if (__term < std::numeric_limits<_Tp>::epsilon())
 329             break;
 330           if (__term >= __prev)
 331             break;
 332           __sum += __term;
 333         }
 334
 335       return std::exp(__x) * __sum / __x;
 336     }
 337
 338
 339     /**
 340      *   @brief Return the exponential integral @f$ Ei(x) @f$.
 341      *
 342      *   The exponential integral is given by
 343      *          \f[
 344      *            Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
 345      *          \f]
 346      *
 347      *   @param  __x  The argument of the exponential integral function.
 348      *   @return  The exponential integral.
 349      */
 350     template<typename _Tp>
 351     _Tp
 352     __expint_Ei(const _Tp __x)
 353     {
 354       if (__x < _Tp(0))
 355         return -__expint_E1(-__x);
 356       else if (__x < -std::log(std::numeric_limits<_Tp>::epsilon()))
 357         return __expint_Ei_series(__x);
 358       else
 359         return __expint_Ei_asymp(__x);
 360     }
 361
 362
 363     /**
 364      *   @brief Return the exponential integral @f$ E_1(x) @f$.
 365      *
 366      *   The exponential integral is given by
 367      *          \f[
 368      *            E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt
 369      *          \f]
 370      *
 371      *   @param  __x  The argument of the exponential integral function.
 372      *   @return  The exponential integral.
 373      */
 374     template<typename _Tp>
 375     _Tp
 376     __expint_E1(const _Tp __x)
 377     {
 378       if (__x < _Tp(0))
 379         return -__expint_Ei(-__x);
 380       else if (__x < _Tp(1))
 381         return __expint_E1_series(__x);
 382       else if (__x < _Tp(100))  //  TODO: Find a good asymptotic switch point.
 383         return __expint_En_cont_frac(1, __x);
 384       else
 385         return __expint_E1_asymp(__x);
 386     }
 387
 388
 389     /**
 390      *   @brief Return the exponential integral @f$ E_n(x) @f$
 391      *          for large argument.
 392      *
 393      *   The exponential integral is given by
 394      *          \f[
 395      *            E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
 396      *          \f]
 397      *
 398      *   This is something of an extension.
 399      *
 400      *   @param  __n  The order of the exponential integral function.
 401      *   @param  __x  The argument of the exponential integral function.
 402      *   @return  The exponential integral.
 403      */
 404     template<typename _Tp>
 405     _Tp
 406     __expint_asymp(const unsigned int __n, const _Tp __x)
 407     {
 408       _Tp __term = _Tp(1);
 409       _Tp __sum = _Tp(1);
 410       for (unsigned int __i = 1; __i <= __n; ++__i)
 411         {
 412           _Tp __prev = __term;
 413           __term *= -(__n - __i + 1) / __x;
 414           if (std::abs(__term) > std::abs(__prev))
 415             break;
 416           __sum += __term;
 417         }
 418
 419       return std::exp(-__x) * __sum / __x;
 420     }
 421
 422
 423     /**
 424      *   @brief Return the exponential integral @f$ E_n(x) @f$
 425      *          for large order.
 426      *
 427      *   The exponential integral is given by
 428      *          \f[
 429      *            E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
 430      *          \f]
 431      *
 432      *   This is something of an extension.
 433      *
 434      *   @param  __n  The order of the exponential integral function.
 435      *   @param  __x  The argument of the exponential integral function.
 436      *   @return  The exponential integral.
 437      */
 438     template<typename _Tp>
 439     _Tp
 440     __expint_large_n(const unsigned int __n, const _Tp __x)
 441     {
 442       const _Tp __xpn = __x + __n;
 443       const _Tp __xpn2 = __xpn * __xpn;
 444       _Tp __term = _Tp(1);
 445       _Tp __sum = _Tp(1);
 446       for (unsigned int __i = 1; __i <= __n; ++__i)
 447         {
 448           _Tp __prev = __term;
 449           __term *= (__n - 2 * (__i - 1) * __x) / __xpn2;
 450           if (std::abs(__term) < std::numeric_limits<_Tp>::epsilon())
 451             break;
 452           __sum += __term;
 453         }
 454
 455       return std::exp(-__x) * __sum / __xpn;
 456     }
 457
 458
 459     /**
 460      *   @brief Return the exponential integral @f$ E_n(x) @f$.
 461      *
 462      *   The exponential integral is given by
 463      *          \f[
 464      *            E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
 465      *          \f]
 466      *   This is something of an extension.
 467      *
 468      *   @param  __n  The order of the exponential integral function.
 469      *   @param  __x  The argument of the exponential integral function.
 470      *   @return  The exponential integral.
 471      */
 472     template<typename _Tp>
 473     _Tp
 474     __expint(const unsigned int __n, const _Tp __x)
 475     {
 476       //  Return NaN on NaN input.
 477       if (__isnan(__x))
 478         return std::numeric_limits<_Tp>::quiet_NaN();
 479       else if (__n <= 1 && __x == _Tp(0))
 480         return std::numeric_limits<_Tp>::infinity();
 481       else
 482         {
 483           _Tp __E0 = std::exp(__x) / __x;
 484           if (__n == 0)
 485             return __E0;
 486
 487           _Tp __E1 = __expint_E1(__x);
 488           if (__n == 1)
 489             return __E1;
 490
 491           if (__x == _Tp(0))
 492             return _Tp(1) / static_cast<_Tp>(__n - 1);
 493
 494           _Tp __En = __expint_En_recursion(__n, __x);
 495
 496           return __En;
 497         }
 498     }
 499
 500
 501     /**
 502      *   @brief Return the exponential integral @f$ Ei(x) @f$.
 503      *
 504      *   The exponential integral is given by
 505      *   \f[
 506      *     Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
 507      *   \f]
 508      *
 509      *   @param  __x  The argument of the exponential integral function.
 510      *   @return  The exponential integral.
 511      */
 512     template<typename _Tp>
 513     inline _Tp
 514     __expint(const _Tp __x)
 515     {
 516       if (__isnan(__x))
 517         return std::numeric_limits<_Tp>::quiet_NaN();
 518       else
 519         return __expint_Ei(__x);
 520     }
 521
 522   _GLIBCXX_END_NAMESPACE_VERSION
 523   } // namespace std::tr1::__detail
 524 }
 525 }
 526
 527 #endif // _GLIBCXX_TR1_EXP_INTEGRAL_TCC