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