libstdc++-v3/include/tr1/beta_function.tcc

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