1// Special functions -*- C++ -*-
2
3// Copyright (C) 2006-2018 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