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
50namespace std
51{
52namespace 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))
188        return std::numeric_limits<_Tp>::quiet_NaN();
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
198