exp_integral.tcc

Go to the documentation of this file.
00001 // Special functions -*- C++ -*-
00002 
00003 // Copyright (C) 2006, 2007, 2008
00004 // Free Software Foundation, Inc.
00005 //
00006 // This file is part of the GNU ISO C++ Library.  This library is free
00007 // software; you can redistribute it and/or modify it under the
00008 // terms of the GNU General Public License as published by the
00009 // Free Software Foundation; either version 2, or (at your option)
00010 // any later version.
00011 //
00012 // This library is distributed in the hope that it will be useful,
00013 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00014 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015 // GNU General Public License for more details.
00016 //
00017 // You should have received a copy of the GNU General Public License along
00018 // with this library; see the file COPYING.  If not, write to the Free
00019 // Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,
00020 // USA.
00021 //
00022 // As a special exception, you may use this file as part of a free software
00023 // library without restriction.  Specifically, if other files instantiate
00024 // templates or use macros or inline functions from this file, or you compile
00025 // this file and link it with other files to produce an executable, this
00026 // file does not by itself cause the resulting executable to be covered by
00027 // the GNU General Public License.  This exception does not however
00028 // invalidate any other reasons why the executable file might be covered by
00029 // the GNU General Public License.
00030 
00031 /** @file tr1/exp_integral.tcc
00032  *  This is an internal header file, included by other library headers.
00033  *  You should not attempt to use it directly.
00034  */
00035 
00036 //
00037 // ISO C++ 14882 TR1: 5.2  Special functions
00038 //
00039 
00040 //  Written by Edward Smith-Rowland based on:
00041 //
00042 //   (1) Handbook of Mathematical Functions,
00043 //       Ed. by Milton Abramowitz and Irene A. Stegun,
00044 //       Dover Publications, New-York, Section 5, pp. 228-251.
00045 //   (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl
00046 //   (3) Numerical Recipes in C, by W. H. Press, S. A. Teukolsky,
00047 //       W. T. Vetterling, B. P. Flannery, Cambridge University Press (1992),
00048 //       2nd ed, pp. 222-225.
00049 //
00050 
00051 #ifndef _GLIBCXX_TR1_EXP_INTEGRAL_TCC
00052 #define _GLIBCXX_TR1_EXP_INTEGRAL_TCC 1
00053 
00054 #include "special_function_util.h"
00055 
00056 namespace std
00057 {
00058 namespace tr1
00059 {
00060 
00061   // [5.2] Special functions
00062 
00063   // Implementation-space details.
00064   namespace __detail
00065   {
00066 
00067     /**
00068      *   @brief Return the exponential integral @f$ E_1(x) @f$
00069      *          by series summation.  This should be good
00070      *          for @f$ x < 1 @f$.
00071      * 
00072      *   The exponential integral is given by
00073      *          \f[
00074      *            E_1(x) = \int_{1}^{\infty} \frac{e^{-xt}}{t} dt
00075      *          \f]
00076      * 
00077      *   @param  __x  The argument of the exponential integral function.
00078      *   @return  The exponential integral.
00079      */
00080     template<typename _Tp>
00081     _Tp
00082     __expint_E1_series(const _Tp __x)
00083     {
00084       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
00085       _Tp __term = _Tp(1);
00086       _Tp __esum = _Tp(0);
00087       _Tp __osum = _Tp(0);
00088       const unsigned int __max_iter = 100;
00089       for (unsigned int __i = 1; __i < __max_iter; ++__i)
00090         {
00091           __term *= - __x / __i;
00092           if (std::abs(__term) < __eps)
00093             break;
00094           if (__term >= _Tp(0))
00095             __esum += __term / __i;
00096           else
00097             __osum += __term / __i;
00098         }
00099 
00100       return - __esum - __osum
00101              - __numeric_constants<_Tp>::__gamma_e() - std::log(__x);
00102     }
00103 
00104 
00105     /**
00106      *   @brief Return the exponential integral @f$ E_1(x) @f$
00107      *          by asymptotic expansion.
00108      * 
00109      *   The exponential integral is given by
00110      *          \f[
00111      *            E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt
00112      *          \f]
00113      * 
00114      *   @param  __x  The argument of the exponential integral function.
00115      *   @return  The exponential integral.
00116      */
00117     template<typename _Tp>
00118     _Tp
00119     __expint_E1_asymp(const _Tp __x)
00120     {
00121       _Tp __term = _Tp(1);
00122       _Tp __esum = _Tp(1);
00123       _Tp __osum = _Tp(0);
00124       const unsigned int __max_iter = 1000;
00125       for (unsigned int __i = 1; __i < __max_iter; ++__i)
00126         {
00127           _Tp __prev = __term;
00128           __term *= - __i / __x;
00129           if (std::abs(__term) > std::abs(__prev))
00130             break;
00131           if (__term >= _Tp(0))
00132             __esum += __term;
00133           else
00134             __osum += __term;
00135         }
00136 
00137       return std::exp(- __x) * (__esum + __osum) / __x;
00138     }
00139 
00140 
00141     /**
00142      *   @brief Return the exponential integral @f$ E_n(x) @f$
00143      *          by series summation.
00144      * 
00145      *   The exponential integral is given by
00146      *          \f[
00147      *            E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
00148      *          \f]
00149      * 
00150      *   @param  __n  The order of the exponential integral function.
00151      *   @param  __x  The argument of the exponential integral function.
00152      *   @return  The exponential integral.
00153      */
00154     template<typename _Tp>
00155     _Tp
00156     __expint_En_series(const unsigned int __n, const _Tp __x)
00157     {
00158       const unsigned int __max_iter = 100;
00159       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
00160       const int __nm1 = __n - 1;
00161       _Tp __ans = (__nm1 != 0
00162                 ? _Tp(1) / __nm1 : -std::log(__x)
00163                                    - __numeric_constants<_Tp>::__gamma_e());
00164       _Tp __fact = _Tp(1);
00165       for (int __i = 1; __i <= __max_iter; ++__i)
00166         {
00167           __fact *= -__x / _Tp(__i);
00168           _Tp __del;
00169           if ( __i != __nm1 )
00170             __del = -__fact / _Tp(__i - __nm1);
00171           else
00172             {
00173               _Tp __psi = -_TR1_GAMMA_TCC;
00174               for (int __ii = 1; __ii <= __nm1; ++__ii)
00175                 __psi += _Tp(1) / _Tp(__ii);
00176               __del = __fact * (__psi - std::log(__x)); 
00177             }
00178           __ans += __del;
00179           if (std::abs(__del) < __eps * std::abs(__ans))
00180             return __ans;
00181         }
00182       std::__throw_runtime_error(__N("Series summation failed "
00183                                      "in __expint_En_series."));
00184     }
00185 
00186 
00187     /**
00188      *   @brief Return the exponential integral @f$ E_n(x) @f$
00189      *          by continued fractions.
00190      * 
00191      *   The exponential integral is given by
00192      *          \f[
00193      *            E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
00194      *          \f]
00195      * 
00196      *   @param  __n  The order of the exponential integral function.
00197      *   @param  __x  The argument of the exponential integral function.
00198      *   @return  The exponential integral.
00199      */
00200     template<typename _Tp>
00201     _Tp
00202     __expint_En_cont_frac(const unsigned int __n, const _Tp __x)
00203     {
00204       const unsigned int __max_iter = 100;
00205       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
00206       const _Tp __fp_min = std::numeric_limits<_Tp>::min();
00207       const int __nm1 = __n - 1;
00208       _Tp __b = __x + _Tp(__n);
00209       _Tp __c = _Tp(1) / __fp_min;
00210       _Tp __d = _Tp(1) / __b;
00211       _Tp __h = __d;
00212       for ( unsigned int __i = 1; __i <= __max_iter; ++__i )
00213         {
00214           _Tp __a = -_Tp(__i * (__nm1 + __i));
00215           __b += _Tp(2);
00216           __d = _Tp(1) / (__a * __d + __b);
00217           __c = __b + __a / __c;
00218           const _Tp __del = __c * __d;
00219           __h *= __del;
00220           if (std::abs(__del - _Tp(1)) < __eps)
00221             {
00222               const _Tp __ans = __h * std::exp(-__x);
00223               return __ans;
00224             }
00225         }
00226       std::__throw_runtime_error(__N("Continued fraction failed "
00227                                      "in __expint_En_cont_frac."));
00228     }
00229 
00230 
00231     /**
00232      *   @brief Return the exponential integral @f$ E_n(x) @f$
00233      *          by recursion.  Use upward recursion for @f$ x < n @f$
00234      *          and downward recursion (Miller's algorithm) otherwise.
00235      * 
00236      *   The exponential integral is given by
00237      *          \f[
00238      *            E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
00239      *          \f]
00240      * 
00241      *   @param  __n  The order of the exponential integral function.
00242      *   @param  __x  The argument of the exponential integral function.
00243      *   @return  The exponential integral.
00244      */
00245     template<typename _Tp>
00246     _Tp
00247     __expint_En_recursion(const unsigned int __n, const _Tp __x)
00248     {
00249       _Tp __En;
00250       _Tp __E1 = __expint_E1(__x);
00251       if (__x < _Tp(__n))
00252         {
00253           //  Forward recursion is stable only for n < x.
00254           __En = __E1;
00255           for (unsigned int __j = 2; __j < __n; ++__j)
00256             __En = (std::exp(-__x) - __x * __En) / _Tp(__j - 1);
00257         }
00258       else
00259         {
00260           //  Backward recursion is stable only for n >= x.
00261           __En = _Tp(1);
00262           const int __N = __n + 20;  //  TODO: Check this starting number.
00263           _Tp __save = _Tp(0);
00264           for (int __j = __N; __j > 0; --__j)
00265             {
00266               __En = (std::exp(-__x) - __j * __En) / __x;
00267               if (__j == __n)
00268                 __save = __En;
00269             }
00270             _Tp __norm = __En / __E1;
00271             __En /= __norm;
00272         }
00273 
00274       return __En;
00275     }
00276 
00277     /**
00278      *   @brief Return the exponential integral @f$ Ei(x) @f$
00279      *          by series summation.
00280      * 
00281      *   The exponential integral is given by
00282      *          \f[
00283      *            Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
00284      *          \f]
00285      * 
00286      *   @param  __x  The argument of the exponential integral function.
00287      *   @return  The exponential integral.
00288      */
00289     template<typename _Tp>
00290     _Tp
00291     __expint_Ei_series(const _Tp __x)
00292     {
00293       _Tp __term = _Tp(1);
00294       _Tp __sum = _Tp(0);
00295       const unsigned int __max_iter = 1000;
00296       for (unsigned int __i = 1; __i < __max_iter; ++__i)
00297         {
00298           __term *= __x / __i;
00299           __sum += __term / __i;
00300           if (__term < std::numeric_limits<_Tp>::epsilon() * __sum)
00301             break;
00302         }
00303 
00304       return __numeric_constants<_Tp>::__gamma_e() + __sum + std::log(__x);
00305     }
00306 
00307 
00308     /**
00309      *   @brief Return the exponential integral @f$ Ei(x) @f$
00310      *          by asymptotic expansion.
00311      * 
00312      *   The exponential integral is given by
00313      *          \f[
00314      *            Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
00315      *          \f]
00316      * 
00317      *   @param  __x  The argument of the exponential integral function.
00318      *   @return  The exponential integral.
00319      */
00320     template<typename _Tp>
00321     _Tp
00322     __expint_Ei_asymp(const _Tp __x)
00323     {
00324       _Tp __term = _Tp(1);
00325       _Tp __sum = _Tp(1);
00326       const unsigned int __max_iter = 1000;
00327       for (unsigned int __i = 1; __i < __max_iter; ++__i)
00328         {
00329           _Tp __prev = __term;
00330           __term *= __i / __x;
00331           if (__term < std::numeric_limits<_Tp>::epsilon())
00332             break;
00333           if (__term >= __prev)
00334             break;
00335           __sum += __term;
00336         }
00337 
00338       return std::exp(__x) * __sum / __x;
00339     }
00340 
00341 
00342     /**
00343      *   @brief Return the exponential integral @f$ Ei(x) @f$.
00344      * 
00345      *   The exponential integral is given by
00346      *          \f[
00347      *            Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
00348      *          \f]
00349      * 
00350      *   @param  __x  The argument of the exponential integral function.
00351      *   @return  The exponential integral.
00352      */
00353     template<typename _Tp>
00354     _Tp
00355     __expint_Ei(const _Tp __x)
00356     {
00357       if (__x < _Tp(0))
00358         return -__expint_E1(-__x);
00359       else if (__x < -std::log(std::numeric_limits<_Tp>::epsilon()))
00360         return __expint_Ei_series(__x);
00361       else
00362         return __expint_Ei_asymp(__x);
00363     }
00364 
00365 
00366     /**
00367      *   @brief Return the exponential integral @f$ E_1(x) @f$.
00368      * 
00369      *   The exponential integral is given by
00370      *          \f[
00371      *            E_1(x) = \int_{1}^\infty \frac{e^{-xt}}{t} dt
00372      *          \f]
00373      * 
00374      *   @param  __x  The argument of the exponential integral function.
00375      *   @return  The exponential integral.
00376      */
00377     template<typename _Tp>
00378     _Tp
00379     __expint_E1(const _Tp __x)
00380     {
00381       if (__x < _Tp(0))
00382         return -__expint_Ei(-__x);
00383       else if (__x < _Tp(1))
00384         return __expint_E1_series(__x);
00385       else if (__x < _Tp(100))  //  TODO: Find a good asymptotic switch point.
00386         return __expint_En_cont_frac(1, __x);
00387       else
00388         return __expint_E1_asymp(__x);
00389     }
00390 
00391 
00392     /**
00393      *   @brief Return the exponential integral @f$ E_n(x) @f$
00394      *          for large argument.
00395      * 
00396      *   The exponential integral is given by
00397      *          \f[
00398      *            E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
00399      *          \f]
00400      * 
00401      *   This is something of an extension.
00402      * 
00403      *   @param  __n  The order of the exponential integral function.
00404      *   @param  __x  The argument of the exponential integral function.
00405      *   @return  The exponential integral.
00406      */
00407     template<typename _Tp>
00408     _Tp
00409     __expint_asymp(const unsigned int __n, const _Tp __x)
00410     {
00411       _Tp __term = _Tp(1);
00412       _Tp __sum = _Tp(1);
00413       for (unsigned int __i = 1; __i <= __n; ++__i)
00414         {
00415           _Tp __prev = __term;
00416           __term *= -(__n - __i + 1) / __x;
00417           if (std::abs(__term) > std::abs(__prev))
00418             break;
00419           __sum += __term;
00420         }
00421 
00422       return std::exp(-__x) * __sum / __x;
00423     }
00424 
00425 
00426     /**
00427      *   @brief Return the exponential integral @f$ E_n(x) @f$
00428      *          for large order.
00429      * 
00430      *   The exponential integral is given by
00431      *          \f[
00432      *            E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
00433      *          \f]
00434      *        
00435      *   This is something of an extension.
00436      * 
00437      *   @param  __n  The order of the exponential integral function.
00438      *   @param  __x  The argument of the exponential integral function.
00439      *   @return  The exponential integral.
00440      */
00441     template<typename _Tp>
00442     _Tp
00443     __expint_large_n(const unsigned int __n, const _Tp __x)
00444     {
00445       const _Tp __xpn = __x + __n;
00446       const _Tp __xpn2 = __xpn * __xpn;
00447       _Tp __term = _Tp(1);
00448       _Tp __sum = _Tp(1);
00449       for (unsigned int __i = 1; __i <= __n; ++__i)
00450         {
00451           _Tp __prev = __term;
00452           __term *= (__n - 2 * (__i - 1) * __x) / __xpn2;
00453           if (std::abs(__term) < std::numeric_limits<_Tp>::epsilon())
00454             break;
00455           __sum += __term;
00456         }
00457 
00458       return std::exp(-__x) * __sum / __xpn;
00459     }
00460 
00461 
00462     /**
00463      *   @brief Return the exponential integral @f$ E_n(x) @f$.
00464      * 
00465      *   The exponential integral is given by
00466      *          \f[
00467      *            E_n(x) = \int_{1}^\infty \frac{e^{-xt}}{t^n} dt
00468      *          \f]
00469      *   This is something of an extension.
00470      * 
00471      *   @param  __n  The order of the exponential integral function.
00472      *   @param  __x  The argument of the exponential integral function.
00473      *   @return  The exponential integral.
00474      */
00475     template<typename _Tp>
00476     _Tp
00477     __expint(const unsigned int __n, const _Tp __x)
00478     {
00479       //  Return NaN on NaN input.
00480       if (__isnan(__x))
00481         return std::numeric_limits<_Tp>::quiet_NaN();
00482       else if (__n <= 1 && __x == _Tp(0))
00483         return std::numeric_limits<_Tp>::infinity();
00484       else
00485         {
00486           _Tp __E0 = std::exp(__x) / __x;
00487           if (__n == 0)
00488             return __E0;
00489 
00490           _Tp __E1 = __expint_E1(__x);
00491           if (__n == 1)
00492             return __E1;
00493 
00494           if (__x == _Tp(0))
00495             return _Tp(1) / static_cast<_Tp>(__n - 1);
00496 
00497           _Tp __En = __expint_En_recursion(__n, __x);
00498 
00499           return __En;
00500         }
00501     }
00502 
00503 
00504     /**
00505      *   @brief Return the exponential integral @f$ Ei(x) @f$.
00506      * 
00507      *   The exponential integral is given by
00508      *   \f[
00509      *     Ei(x) = -\int_{-x}^\infty \frac{e^t}{t} dt
00510      *   \f]
00511      * 
00512      *   @param  __x  The argument of the exponential integral function.
00513      *   @return  The exponential integral.
00514      */
00515     template<typename _Tp>
00516     inline _Tp
00517     __expint(const _Tp __x)
00518     {
00519       if (__isnan(__x))
00520         return std::numeric_limits<_Tp>::quiet_NaN();
00521       else
00522         return __expint_Ei(__x);
00523     }
00524 
00525   } // namespace std::tr1::__detail
00526 }
00527 }
00528 
00529 #endif // _GLIBCXX_TR1_EXP_INTEGRAL_TCC

Generated on Wed Dec 31 12:48:54 2008 for libstdc++ by  doxygen 1.5.6