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 50namespace std 51{ 52namespace 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