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