1 // Special functions -*- C++ -*- 2 3 // Copyright (C) 2006, 2007, 2008, 2009 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 3, 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 // Under Section 7 of GPL version 3, you are granted additional 18 // permissions described in the GCC Runtime Library Exception, version 19 // 3.1, as published by the Free Software Foundation. 20 21 // You should have received a copy of the GNU General Public License and 22 // a copy of the GCC Runtime Library Exception along with this program; 23 // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 24 // <http://www.gnu.org/licenses/>. 25 26 /** @file tr1/beta_function.tcc 27 * This is an internal header file, included by other library headers. 28 * You should not attempt to use it directly. 29 */ 30 31 // 32 // ISO C++ 14882 TR1: 5.2 Special functions 33 // 34 35 // Written by Edward Smith-Rowland based on: 36 // (1) Handbook of Mathematical Functions, 37 // ed. Milton Abramowitz and Irene A. Stegun, 38 // Dover Publications, 39 // Section 6, pp. 253-266 40 // (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl 41 // (3) Numerical Recipes in C, by W. H. Press, S. A. Teukolsky, 42 // W. T. Vetterling, B. P. Flannery, Cambridge University Press (1992), 43 // 2nd ed, pp. 213-216 44 // (4) Gamma, Exploring Euler's Constant, Julian Havil, 45 // Princeton, 2003. 46 47 #ifndef _GLIBCXX_TR1_BETA_FUNCTION_TCC 48 #define _GLIBCXX_TR1_BETA_FUNCTION_TCC 1 49 50 namespace std 51 { 52 namespace tr1 53 { 54 55 // [5.2] Special functions 56 57 // Implementation-space details. 58 namespace __detail 59 { 60 61 /** 62 * @brief Return the beta function: \f$B(x,y)\f$. 63 * 64 * The beta function is defined by 65 * @f[ 66 * B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)} 67 * @f] 68 * 69 * @param __x The first argument of the beta function. 70 * @param __y The second argument of the beta function. 71 * @return The beta function. 72 */ 73 template<typename _Tp> 74 _Tp 75 __beta_gamma(_Tp __x, _Tp __y) 76 { 77 78 _Tp __bet; 79 #if _GLIBCXX_USE_C99_MATH_TR1 80 if (__x > __y) 81 { 82 __bet = std::tr1::tgamma(__x) 83 / std::tr1::tgamma(__x + __y); 84 __bet *= std::tr1::tgamma(__y); 85 } 86 else 87 { 88 __bet = std::tr1::tgamma(__y) 89 / std::tr1::tgamma(__x + __y); 90 __bet *= std::tr1::tgamma(__x); 91 } 92 #else 93 if (__x > __y) 94 { 95 __bet = __gamma(__x) / __gamma(__x + __y); 96 __bet *= __gamma(__y); 97 } 98 else 99 { 100 __bet = __gamma(__y) / __gamma(__x + __y); 101 __bet *= __gamma(__x); 102 } 103 #endif 104 105 return __bet; 106 } 107 108 /** 109 * @brief Return the beta function \f$B(x,y)\f$ using 110 * the log gamma functions. 111 * 112 * The beta function is defined by 113 * @f[ 114 * B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)} 115 * @f] 116 * 117 * @param __x The first argument of the beta function. 118 * @param __y The second argument of the beta function. 119 * @return The beta function. 120 */ 121 template<typename _Tp> 122 _Tp 123 __beta_lgamma(_Tp __x, _Tp __y) 124 { 125 #if _GLIBCXX_USE_C99_MATH_TR1 126 _Tp __bet = std::tr1::lgamma(__x) 127 + std::tr1::lgamma(__y) 128 - std::tr1::lgamma(__x + __y); 129 #else 130 _Tp __bet = __log_gamma(__x) 131 + __log_gamma(__y) 132 - __log_gamma(__x + __y); 133 #endif 134 __bet = std::exp(__bet); 135 return __bet; 136 } 137 138 139 /** 140 * @brief Return the beta function \f$B(x,y)\f$ using 141 * the product form. 142 * 143 * The beta function is defined by 144 * @f[ 145 * B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)} 146 * @f] 147 * 148 * @param __x The first argument of the beta function. 149 * @param __y The second argument of the beta function. 150 * @return The beta function. 151 */ 152 template<typename _Tp> 153 _Tp 154 __beta_product(_Tp __x, _Tp __y) 155 { 156 157 _Tp __bet = (__x + __y) / (__x * __y); 158 159 unsigned int __max_iter = 1000000; 160 for (unsigned int __k = 1; __k < __max_iter; ++__k) 161 { 162 _Tp __term = (_Tp(1) + (__x + __y) / __k) 163 / ((_Tp(1) + __x / __k) * (_Tp(1) + __y / __k)); 164 __bet *= __term; 165 } 166 167 return __bet; 168 } 169 170 171 /** 172 * @brief Return the beta function \f$ B(x,y) \f$. 173 * 174 * The beta function is defined by 175 * @f[ 176 * B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)} 177 * @f] 178 * 179 * @param __x The first argument of the beta function. 180 * @param __y The second argument of the beta function. 181 * @return The beta function. 182 */ 183 template<typename _Tp> 184 inline _Tp 185 __beta(_Tp __x, _Tp __y) 186 { 187 if (__isnan(__x) || __isnan(__y)) 188 return std::numeric_limits<_Tp>::quiet_NaN(); 189 else 190 return __beta_lgamma(__x, __y); 191 } 192 193 } // namespace std::tr1::__detail 194 } 195 } 196 197 #endif // __GLIBCXX_TR1_BETA_FUNCTION_TCC 198
