exp_integral.tcc

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

Generated on Thu Jul 23 21:16:04 2009 for libstdc++ by  doxygen 1.5.8