contrib/gcc-4.4/libstdc++-v3/include/tr1/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