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

   1 // Special functions -*- C++ -*-
   2
   3 // Copyright (C) 2006-2015 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 namespace tr1
  52 {
  53   // [5.2] Special functions
  54
  55   // Implementation-space details.
  56   namespace __detail
  57   {
  58   _GLIBCXX_BEGIN_NAMESPACE_VERSION
  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         {
  81           __bet = std::tr1::tgamma(__x)
  82                 / std::tr1::tgamma(__x + __y);
  83           __bet *= std::tr1::tgamma(__y);
  84         }
  85       else
  86         {
  87           __bet = std::tr1::tgamma(__y)
  88                 / std::tr1::tgamma(__x + __y);
  89           __bet *= std::tr1::tgamma(__x);
  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
 125       _Tp __bet = std::tr1::lgamma(__x)
 126                 + std::tr1::lgamma(__y)
 127                 - std::tr1::lgamma(__x + __y);
 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
 192   _GLIBCXX_END_NAMESPACE_VERSION
 193   } // namespace std::tr1::__detail
 194 }
 195 }
 196
 197 #endif // _GLIBCXX_TR1_BETA_FUNCTION_TCC