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