libstdc++
beta_function.tcc
Go to the documentation of this file.
1 // Special functions -*- C++ -*-
2 
3 // Copyright (C) 2006, 2007, 2008, 2009
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 3, 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 // 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.
20 
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/>.
25 
26 /** @file tr1/beta_function.tcc
27  * This is an internal header file, included by other library headers.
28  * You should not attempt to use it directly.
29  */
30 
31 //
32 // ISO C++ 14882 TR1: 5.2 Special functions
33 //
34 
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,
45 // Princeton, 2003.
46 
47 #ifndef _GLIBCXX_TR1_BETA_FUNCTION_TCC
48 #define _GLIBCXX_TR1_BETA_FUNCTION_TCC 1
49 
50 namespace std
51 {
52 namespace tr1
53 {
54 
55  // [5.2] Special functions
56 
57  // Implementation-space details.
58  namespace __detail
59  {
60 
61  /**
62  * @brief Return the beta function: \f$B(x,y)\f$.
63  *
64  * The beta function is defined by
65  * @f[
66  * B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
67  * @f]
68  *
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.
72  */
73  template<typename _Tp>
74  _Tp
75  __beta_gamma(_Tp __x, _Tp __y)
76  {
77 
78  _Tp __bet;
79 #if _GLIBCXX_USE_C99_MATH_TR1
80  if (__x > __y)
81  {
82  __bet = std::tr1::tgamma(__x)
83  / std::tr1::tgamma(__x + __y);
84  __bet *= std::tr1::tgamma(__y);
85  }
86  else
87  {
88  __bet = std::tr1::tgamma(__y)
89  / std::tr1::tgamma(__x + __y);
90  __bet *= std::tr1::tgamma(__x);
91  }
92 #else
93  if (__x > __y)
94  {
95  __bet = __gamma(__x) / __gamma(__x + __y);
96  __bet *= __gamma(__y);
97  }
98  else
99  {
100  __bet = __gamma(__y) / __gamma(__x + __y);
101  __bet *= __gamma(__x);
102  }
103 #endif
104 
105  return __bet;
106  }
107 
108  /**
109  * @brief Return the beta function \f$B(x,y)\f$ using
110  * the log gamma functions.
111  *
112  * The beta function is defined by
113  * @f[
114  * B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
115  * @f]
116  *
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.
120  */
121  template<typename _Tp>
122  _Tp
123  __beta_lgamma(_Tp __x, _Tp __y)
124  {
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);
129 #else
130  _Tp __bet = __log_gamma(__x)
131  + __log_gamma(__y)
132  - __log_gamma(__x + __y);
133 #endif
134  __bet = std::exp(__bet);
135  return __bet;
136  }
137 
138 
139  /**
140  * @brief Return the beta function \f$B(x,y)\f$ using
141  * the product form.
142  *
143  * The beta function is defined by
144  * @f[
145  * B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
146  * @f]
147  *
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.
151  */
152  template<typename _Tp>
153  _Tp
154  __beta_product(_Tp __x, _Tp __y)
155  {
156 
157  _Tp __bet = (__x + __y) / (__x * __y);
158 
159  unsigned int __max_iter = 1000000;
160  for (unsigned int __k = 1; __k < __max_iter; ++__k)
161  {
162  _Tp __term = (_Tp(1) + (__x + __y) / __k)
163  / ((_Tp(1) + __x / __k) * (_Tp(1) + __y / __k));
164  __bet *= __term;
165  }
166 
167  return __bet;
168  }
169 
170 
171  /**
172  * @brief Return the beta function \f$ B(x,y) \f$.
173  *
174  * The beta function is defined by
175  * @f[
176  * B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
177  * @f]
178  *
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.
182  */
183  template<typename _Tp>
184  inline _Tp
185  __beta(_Tp __x, _Tp __y)
186  {
187  if (__isnan(__x) || __isnan(__y))
189  else
190  return __beta_lgamma(__x, __y);
191  }
192 
193  } // namespace std::tr1::__detail
194 }
195 }
196 
197 #endif // __GLIBCXX_TR1_BETA_FUNCTION_TCC