1 // Special functions -*- C++ -*-
3 // Copyright (C) 2006, 2007, 2008, 2009, 2010
4 // Free Software Foundation, Inc.
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)
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.
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.
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/>.
26 /** @file tr1/beta_function.tcc
27 * This is an internal header file, included by other library headers.
28 * Do not attempt to use it directly. @headername{tr1/cmath}
32 // ISO C++ 14882 TR1: 5.2 Special functions
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,
47 #ifndef _GLIBCXX_TR1_BETA_FUNCTION_TCC
48 #define _GLIBCXX_TR1_BETA_FUNCTION_TCC 1
50 namespace std _GLIBCXX_VISIBILITY(default)
54 // [5.2] Special functions
56 // Implementation-space details.
59 _GLIBCXX_BEGIN_NAMESPACE_VERSION
62 * @brief Return the beta function: \f$B(x,y)\f$.
64 * The beta function is defined by
66 * B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
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.
73 template<typename _Tp>
75 __beta_gamma(_Tp __x, _Tp __y)
79 #if _GLIBCXX_USE_C99_MATH_TR1
82 __bet = std::tr1::tgamma(__x)
83 / std::tr1::tgamma(__x + __y);
84 __bet *= std::tr1::tgamma(__y);
88 __bet = std::tr1::tgamma(__y)
89 / std::tr1::tgamma(__x + __y);
90 __bet *= std::tr1::tgamma(__x);
95 __bet = __gamma(__x) / __gamma(__x + __y);
96 __bet *= __gamma(__y);
100 __bet = __gamma(__y) / __gamma(__x + __y);
101 __bet *= __gamma(__x);
109 * @brief Return the beta function \f$B(x,y)\f$ using
110 * the log gamma functions.
112 * The beta function is defined by
114 * B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
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.
121 template<typename _Tp>
123 __beta_lgamma(_Tp __x, _Tp __y)
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);
130 _Tp __bet = __log_gamma(__x)
132 - __log_gamma(__x + __y);
134 __bet = std::exp(__bet);
140 * @brief Return the beta function \f$B(x,y)\f$ using
143 * The beta function is defined by
145 * B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
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.
152 template<typename _Tp>
154 __beta_product(_Tp __x, _Tp __y)
157 _Tp __bet = (__x + __y) / (__x * __y);
159 unsigned int __max_iter = 1000000;
160 for (unsigned int __k = 1; __k < __max_iter; ++__k)
162 _Tp __term = (_Tp(1) + (__x + __y) / __k)
163 / ((_Tp(1) + __x / __k) * (_Tp(1) + __y / __k));
172 * @brief Return the beta function \f$ B(x,y) \f$.
174 * The beta function is defined by
176 * B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
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.
183 template<typename _Tp>
185 __beta(_Tp __x, _Tp __y)
187 if (__isnan(__x) || __isnan(__y))
188 return std::numeric_limits<_Tp>::quiet_NaN();
190 return __beta_lgamma(__x, __y);
193 _GLIBCXX_END_NAMESPACE_VERSION
194 } // namespace std::tr1::__detail
198 #endif // __GLIBCXX_TR1_BETA_FUNCTION_TCC