]>
Commit | Line | Data |
---|---|---|
7c62b943 BK |
1 | // Special functions -*- C++ -*- |
2 | ||
aa118a03 | 3 | // Copyright (C) 2006-2014 Free Software Foundation, Inc. |
7c62b943 BK |
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 | |
748086b7 | 8 | // Free Software Foundation; either version 3, or (at your option) |
7c62b943 BK |
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 | // | |
748086b7 JJ |
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/>. | |
7c62b943 BK |
24 | |
25 | /** @file tr1/beta_function.tcc | |
26 | * This is an internal header file, included by other library headers. | |
f910786b | 27 | * Do not attempt to use it directly. @headername{tr1/cmath} |
7c62b943 BK |
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 | ||
e133ace8 PC |
46 | #ifndef _GLIBCXX_TR1_BETA_FUNCTION_TCC |
47 | #define _GLIBCXX_TR1_BETA_FUNCTION_TCC 1 | |
7c62b943 | 48 | |
12ffa228 | 49 | namespace std _GLIBCXX_VISIBILITY(default) |
7c62b943 | 50 | { |
e133ace8 PC |
51 | namespace tr1 |
52 | { | |
7c62b943 BK |
53 | // [5.2] Special functions |
54 | ||
7c62b943 | 55 | // Implementation-space details. |
7c62b943 BK |
56 | namespace __detail |
57 | { | |
12ffa228 | 58 | _GLIBCXX_BEGIN_NAMESPACE_VERSION |
7c62b943 BK |
59 | |
60 | /** | |
61 | * @brief Return the beta function: \f$B(x,y)\f$. | |
62 | * | |
63 | * The beta function is defined by | |
64 | * @f[ | |
65 | * B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)} | |
66 | * @f] | |
67 | * | |
68 | * @param __x The first argument of the beta function. | |
69 | * @param __y The second argument of the beta function. | |
70 | * @return The beta function. | |
71 | */ | |
72 | template<typename _Tp> | |
73 | _Tp | |
74 | __beta_gamma(_Tp __x, _Tp __y) | |
75 | { | |
76 | ||
77 | _Tp __bet; | |
78 | #if _GLIBCXX_USE_C99_MATH_TR1 | |
79 | if (__x > __y) | |
80 | { | |
e133ace8 PC |
81 | __bet = std::tr1::tgamma(__x) |
82 | / std::tr1::tgamma(__x + __y); | |
83 | __bet *= std::tr1::tgamma(__y); | |
7c62b943 BK |
84 | } |
85 | else | |
86 | { | |
e133ace8 PC |
87 | __bet = std::tr1::tgamma(__y) |
88 | / std::tr1::tgamma(__x + __y); | |
89 | __bet *= std::tr1::tgamma(__x); | |
7c62b943 BK |
90 | } |
91 | #else | |
92 | if (__x > __y) | |
93 | { | |
94 | __bet = __gamma(__x) / __gamma(__x + __y); | |
95 | __bet *= __gamma(__y); | |
96 | } | |
97 | else | |
98 | { | |
99 | __bet = __gamma(__y) / __gamma(__x + __y); | |
100 | __bet *= __gamma(__x); | |
101 | } | |
102 | #endif | |
103 | ||
104 | return __bet; | |
105 | } | |
106 | ||
107 | /** | |
108 | * @brief Return the beta function \f$B(x,y)\f$ using | |
109 | * the log gamma functions. | |
110 | * | |
111 | * The beta function is defined by | |
112 | * @f[ | |
113 | * B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)} | |
114 | * @f] | |
115 | * | |
116 | * @param __x The first argument of the beta function. | |
117 | * @param __y The second argument of the beta function. | |
118 | * @return The beta function. | |
119 | */ | |
120 | template<typename _Tp> | |
121 | _Tp | |
122 | __beta_lgamma(_Tp __x, _Tp __y) | |
123 | { | |
124 | #if _GLIBCXX_USE_C99_MATH_TR1 | |
e133ace8 PC |
125 | _Tp __bet = std::tr1::lgamma(__x) |
126 | + std::tr1::lgamma(__y) | |
127 | - std::tr1::lgamma(__x + __y); | |
7c62b943 BK |
128 | #else |
129 | _Tp __bet = __log_gamma(__x) | |
130 | + __log_gamma(__y) | |
131 | - __log_gamma(__x + __y); | |
132 | #endif | |
133 | __bet = std::exp(__bet); | |
134 | return __bet; | |
135 | } | |
136 | ||
137 | ||
138 | /** | |
139 | * @brief Return the beta function \f$B(x,y)\f$ using | |
140 | * the product form. | |
141 | * | |
142 | * The beta function is defined by | |
143 | * @f[ | |
144 | * B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)} | |
145 | * @f] | |
146 | * | |
147 | * @param __x The first argument of the beta function. | |
148 | * @param __y The second argument of the beta function. | |
149 | * @return The beta function. | |
150 | */ | |
151 | template<typename _Tp> | |
152 | _Tp | |
153 | __beta_product(_Tp __x, _Tp __y) | |
154 | { | |
155 | ||
156 | _Tp __bet = (__x + __y) / (__x * __y); | |
157 | ||
158 | unsigned int __max_iter = 1000000; | |
159 | for (unsigned int __k = 1; __k < __max_iter; ++__k) | |
160 | { | |
161 | _Tp __term = (_Tp(1) + (__x + __y) / __k) | |
162 | / ((_Tp(1) + __x / __k) * (_Tp(1) + __y / __k)); | |
163 | __bet *= __term; | |
164 | } | |
165 | ||
166 | return __bet; | |
167 | } | |
168 | ||
169 | ||
170 | /** | |
171 | * @brief Return the beta function \f$ B(x,y) \f$. | |
172 | * | |
173 | * The beta function is defined by | |
174 | * @f[ | |
175 | * B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)} | |
176 | * @f] | |
177 | * | |
178 | * @param __x The first argument of the beta function. | |
179 | * @param __y The second argument of the beta function. | |
180 | * @return The beta function. | |
181 | */ | |
182 | template<typename _Tp> | |
183 | inline _Tp | |
184 | __beta(_Tp __x, _Tp __y) | |
185 | { | |
186 | if (__isnan(__x) || __isnan(__y)) | |
187 | return std::numeric_limits<_Tp>::quiet_NaN(); | |
188 | else | |
189 | return __beta_lgamma(__x, __y); | |
190 | } | |
191 | ||
12ffa228 | 192 | _GLIBCXX_END_NAMESPACE_VERSION |
7c62b943 | 193 | } // namespace std::tr1::__detail |
e133ace8 | 194 | } |
7c62b943 BK |
195 | } |
196 | ||
e133ace8 | 197 | #endif // __GLIBCXX_TR1_BETA_FUNCTION_TCC |