l4/pkg/libstdc++-v3/contrib/libstdc++-v3-4.3.3/include/tr1/exp_integral.tcc

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