libstdc++-v3/include/tr1/exp_integral.tcc

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