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