beta_function.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/beta_function.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 //   (1) Handbook of Mathematical Functions,
00037 //       ed. Milton Abramowitz and Irene A. Stegun,
00038 //       Dover Publications,
00039 //       Section 6, pp. 253-266
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. 213-216
00044 //   (4) Gamma, Exploring Euler's Constant, Julian Havil,
00045 //       Princeton, 2003.
00046 
00047 #ifndef _GLIBCXX_TR1_BETA_FUNCTION_TCC
00048 #define _GLIBCXX_TR1_BETA_FUNCTION_TCC 1
00049 
00050 namespace std
00051 {
00052 namespace tr1
00053 {
00054 
00055   // [5.2] Special functions
00056 
00057   // Implementation-space details.
00058   namespace __detail
00059   {
00060 
00061     /**
00062      *   @brief  Return the beta function: \f$B(x,y)\f$.
00063      * 
00064      *   The beta function is defined by
00065      *   @f[
00066      *     B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
00067      *   @f]
00068      *
00069      *   @param __x The first argument of the beta function.
00070      *   @param __y The second argument of the beta function.
00071      *   @return  The beta function.
00072      */
00073     template<typename _Tp>
00074     _Tp
00075     __beta_gamma(_Tp __x, _Tp __y)
00076     {
00077 
00078       _Tp __bet;
00079 #if _GLIBCXX_USE_C99_MATH_TR1
00080       if (__x > __y)
00081         {
00082           __bet = std::tr1::tgamma(__x)
00083                 / std::tr1::tgamma(__x + __y);
00084           __bet *= std::tr1::tgamma(__y);
00085         }
00086       else
00087         {
00088           __bet = std::tr1::tgamma(__y)
00089                 / std::tr1::tgamma(__x + __y);
00090           __bet *= std::tr1::tgamma(__x);
00091         }
00092 #else
00093       if (__x > __y)
00094         {
00095           __bet = __gamma(__x) / __gamma(__x + __y);
00096           __bet *= __gamma(__y);
00097         }
00098       else
00099         {
00100           __bet = __gamma(__y) / __gamma(__x + __y);
00101           __bet *= __gamma(__x);
00102         }
00103 #endif
00104 
00105       return __bet;
00106     }
00107 
00108     /**
00109      *   @brief  Return the beta function \f$B(x,y)\f$ using
00110      *           the log gamma functions.
00111      * 
00112      *   The beta function is defined by
00113      *   @f[
00114      *     B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
00115      *   @f]
00116      *
00117      *   @param __x The first argument of the beta function.
00118      *   @param __y The second argument of the beta function.
00119      *   @return  The beta function.
00120      */
00121     template<typename _Tp>
00122     _Tp
00123     __beta_lgamma(_Tp __x, _Tp __y)
00124     {
00125 #if _GLIBCXX_USE_C99_MATH_TR1
00126       _Tp __bet = std::tr1::lgamma(__x)
00127                 + std::tr1::lgamma(__y)
00128                 - std::tr1::lgamma(__x + __y);
00129 #else
00130       _Tp __bet = __log_gamma(__x)
00131                 + __log_gamma(__y)
00132                 - __log_gamma(__x + __y);
00133 #endif
00134       __bet = std::exp(__bet);
00135       return __bet;
00136     }
00137 
00138 
00139     /**
00140      *   @brief  Return the beta function \f$B(x,y)\f$ using
00141      *           the product form.
00142      * 
00143      *   The beta function is defined by
00144      *   @f[
00145      *     B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
00146      *   @f]
00147      *
00148      *   @param __x The first argument of the beta function.
00149      *   @param __y The second argument of the beta function.
00150      *   @return  The beta function.
00151      */
00152     template<typename _Tp>
00153     _Tp
00154     __beta_product(_Tp __x, _Tp __y)
00155     {
00156 
00157       _Tp __bet = (__x + __y) / (__x * __y);
00158 
00159       unsigned int __max_iter = 1000000;
00160       for (unsigned int __k = 1; __k < __max_iter; ++__k)
00161         {
00162           _Tp __term = (_Tp(1) + (__x + __y) / __k)
00163                      / ((_Tp(1) + __x / __k) * (_Tp(1) + __y / __k));
00164           __bet *= __term;
00165         }
00166 
00167       return __bet;
00168     }
00169 
00170 
00171     /**
00172      *   @brief  Return the beta function \f$ B(x,y) \f$.
00173      * 
00174      *   The beta function is defined by
00175      *   @f[
00176      *     B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
00177      *   @f]
00178      *
00179      *   @param __x The first argument of the beta function.
00180      *   @param __y The second argument of the beta function.
00181      *   @return  The beta function.
00182      */
00183     template<typename _Tp>
00184     inline _Tp
00185     __beta(_Tp __x, _Tp __y)
00186     {
00187       if (__isnan(__x) || __isnan(__y))
00188         return std::numeric_limits<_Tp>::quiet_NaN();
00189       else
00190         return __beta_lgamma(__x, __y);
00191     }
00192 
00193   } // namespace std::tr1::__detail
00194 }
00195 }
00196 
00197 #endif // __GLIBCXX_TR1_BETA_FUNCTION_TCC

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