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

   1 // Special functions -*- C++ -*-
   2
   3 // Copyright (C) 2006-2007
   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 2, 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 // You should have received a copy of the GNU General Public License along
  18 // with this library; see the file COPYING.  If not, write to the Free
  19 // Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,
  20 // USA.
  21 //
  22 // As a special exception, you may use this file as part of a free software
  23 // library without restriction.  Specifically, if other files instantiate
  24 // templates or use macros or inline functions from this file, or you compile
  25 // this file and link it with other files to produce an executable, this
  26 // file does not by itself cause the resulting executable to be covered by
  27 // the GNU General Public License.  This exception does not however
  28 // invalidate any other reasons why the executable file might be covered by
  29 // the GNU General Public License.
  30
  31 /** @file tr1/beta_function.tcc
  32  *  This is an internal header file, included by other library headers.
  33  *  You should not attempt to use it directly.
  34  */
  35
  36 //
  37 // ISO C++ 14882 TR1: 5.2  Special functions
  38 //
  39
  40 // Written by Edward Smith-Rowland based on:
  41 //   (1) Handbook of Mathematical Functions,
  42 //       ed. Milton Abramowitz and Irene A. Stegun,
  43 //       Dover Publications,
  44 //       Section 6, pp. 253-266
  45 //   (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl
  46 //   (3) Numerical Recipes in C, by W. H. Press, S. A. Teukolsky,
  47 //       W. T. Vetterling, B. P. Flannery, Cambridge University Press (1992),
  48 //       2nd ed, pp. 213-216
  49 //   (4) Gamma, Exploring Euler's Constant, Julian Havil,
  50 //       Princeton, 2003.
  51
  52 #ifndef _TR1_BETA_FUNCTION_TCC
  53 #define _TR1_BETA_FUNCTION_TCC 1
  54
  55 namespace std
  56 {
  57 _GLIBCXX_BEGIN_NAMESPACE(_GLIBCXX_TR1)
  58
  59   // [5.2] Special functions
  60
  61   /**
  62    * @ingroup tr1_math_spec_func
  63    * @{
  64    */
  65
  66   //
  67   // Implementation-space details.
  68   //
  69   namespace __detail
  70   {
  71
  72     /**
  73      *   @brief  Return the beta function: \f$B(x,y)\f$.
  74      *
  75      *   The beta function is defined by
  76      *   @f[
  77      *     B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
  78      *   @f]
  79      *
  80      *   @param __x The first argument of the beta function.
  81      *   @param __y The second argument of the beta function.
  82      *   @return  The beta function.
  83      */
  84     template<typename _Tp>
  85     _Tp
  86     __beta_gamma(_Tp __x, _Tp __y)
  87     {
  88
  89       _Tp __bet;
  90 #if _GLIBCXX_USE_C99_MATH_TR1
  91       if (__x > __y)
  92         {
  93           __bet = std::_GLIBCXX_TR1::tgamma(__x)
  94                 / std::_GLIBCXX_TR1::tgamma(__x + __y);
  95           __bet *= std::_GLIBCXX_TR1::tgamma(__y);
  96         }
  97       else
  98         {
  99           __bet = std::_GLIBCXX_TR1::tgamma(__y)
 100                 / std::_GLIBCXX_TR1::tgamma(__x + __y);
 101           __bet *= std::_GLIBCXX_TR1::tgamma(__x);
 102         }
 103 #else
 104       if (__x > __y)
 105         {
 106           __bet = __gamma(__x) / __gamma(__x + __y);
 107           __bet *= __gamma(__y);
 108         }
 109       else
 110         {
 111           __bet = __gamma(__y) / __gamma(__x + __y);
 112           __bet *= __gamma(__x);
 113         }
 114 #endif
 115
 116       return __bet;
 117     }
 118
 119     /**
 120      *   @brief  Return the beta function \f$B(x,y)\f$ using
 121      *           the log gamma functions.
 122      *
 123      *   The beta function is defined by
 124      *   @f[
 125      *     B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
 126      *   @f]
 127      *
 128      *   @param __x The first argument of the beta function.
 129      *   @param __y The second argument of the beta function.
 130      *   @return  The beta function.
 131      */
 132     template<typename _Tp>
 133     _Tp
 134     __beta_lgamma(_Tp __x, _Tp __y)
 135     {
 136 #if _GLIBCXX_USE_C99_MATH_TR1
 137       _Tp __bet = std::_GLIBCXX_TR1::lgamma(__x)
 138                 + std::_GLIBCXX_TR1::lgamma(__y)
 139                 - std::_GLIBCXX_TR1::lgamma(__x + __y);
 140 #else
 141       _Tp __bet = __log_gamma(__x)
 142                 + __log_gamma(__y)
 143                 - __log_gamma(__x + __y);
 144 #endif
 145       __bet = std::exp(__bet);
 146       return __bet;
 147     }
 148
 149
 150     /**
 151      *   @brief  Return the beta function \f$B(x,y)\f$ using
 152      *           the product form.
 153      *
 154      *   The beta function is defined by
 155      *   @f[
 156      *     B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
 157      *   @f]
 158      *
 159      *   @param __x The first argument of the beta function.
 160      *   @param __y The second argument of the beta function.
 161      *   @return  The beta function.
 162      */
 163     template<typename _Tp>
 164     _Tp
 165     __beta_product(_Tp __x, _Tp __y)
 166     {
 167
 168       _Tp __bet = (__x + __y) / (__x * __y);
 169
 170       unsigned int __max_iter = 1000000;
 171       for (unsigned int __k = 1; __k < __max_iter; ++__k)
 172         {
 173           _Tp __term = (_Tp(1) + (__x + __y) / __k)
 174                      / ((_Tp(1) + __x / __k) * (_Tp(1) + __y / __k));
 175           __bet *= __term;
 176         }
 177
 178       return __bet;
 179     }
 180
 181
 182     /**
 183      *   @brief  Return the beta function \f$ B(x,y) \f$.
 184      *
 185      *   The beta function is defined by
 186      *   @f[
 187      *     B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
 188      *   @f]
 189      *
 190      *   @param __x The first argument of the beta function.
 191      *   @param __y The second argument of the beta function.
 192      *   @return  The beta function.
 193      */
 194     template<typename _Tp>
 195     inline _Tp
 196     __beta(_Tp __x, _Tp __y)
 197     {
 198       if (__isnan(__x) || __isnan(__y))
 199         return std::numeric_limits<_Tp>::quiet_NaN();
 200       else
 201         return __beta_lgamma(__x, __y);
 202     }
 203
 204   } // namespace std::tr1::__detail
 205
 206   /* @} */ // group tr1_math_spec_func
 207
 208 _GLIBCXX_END_NAMESPACE
 209 }
 210
 211 #endif // _TR1_BETA_FUNCTION_TCC