mingw/lib/gcc/mingw32/4.3.3/include/c++/tr1/beta_function.tcc

   1 // Special functions -*- C++ -*-
   2
   3 // Copyright (C) 2006, 2007, 2008
   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 _GLIBCXX_TR1_BETA_FUNCTION_TCC
  53 #define _GLIBCXX_TR1_BETA_FUNCTION_TCC 1
  54
  55 namespace std
  56 {
  57 namespace tr1
  58 {
  59
  60   // [5.2] Special functions
  61
  62   // Implementation-space details.
  63   namespace __detail
  64   {
  65
  66     /**
  67      *   @brief  Return the beta function: \f$B(x,y)\f$.
  68      *
  69      *   The beta function is defined by
  70      *   @f[
  71      *     B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
  72      *   @f]
  73      *
  74      *   @param __x The first argument of the beta function.
  75      *   @param __y The second argument of the beta function.
  76      *   @return  The beta function.
  77      */
  78     template<typename _Tp>
  79     _Tp
  80     __beta_gamma(_Tp __x, _Tp __y)
  81     {
  82
  83       _Tp __bet;
  84 #if _GLIBCXX_USE_C99_MATH_TR1
  85       if (__x > __y)
  86         {
  87           __bet = std::tr1::tgamma(__x)
  88                 / std::tr1::tgamma(__x + __y);
  89           __bet *= std::tr1::tgamma(__y);
  90         }
  91       else
  92         {
  93           __bet = std::tr1::tgamma(__y)
  94                 / std::tr1::tgamma(__x + __y);
  95           __bet *= std::tr1::tgamma(__x);
  96         }
  97 #else
  98       if (__x > __y)
  99         {
 100           __bet = __gamma(__x) / __gamma(__x + __y);
 101           __bet *= __gamma(__y);
 102         }
 103       else
 104         {
 105           __bet = __gamma(__y) / __gamma(__x + __y);
 106           __bet *= __gamma(__x);
 107         }
 108 #endif
 109
 110       return __bet;
 111     }
 112
 113     /**
 114      *   @brief  Return the beta function \f$B(x,y)\f$ using
 115      *           the log gamma functions.
 116      *
 117      *   The beta function is defined by
 118      *   @f[
 119      *     B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
 120      *   @f]
 121      *
 122      *   @param __x The first argument of the beta function.
 123      *   @param __y The second argument of the beta function.
 124      *   @return  The beta function.
 125      */
 126     template<typename _Tp>
 127     _Tp
 128     __beta_lgamma(_Tp __x, _Tp __y)
 129     {
 130 #if _GLIBCXX_USE_C99_MATH_TR1
 131       _Tp __bet = std::tr1::lgamma(__x)
 132                 + std::tr1::lgamma(__y)
 133                 - std::tr1::lgamma(__x + __y);
 134 #else
 135       _Tp __bet = __log_gamma(__x)
 136                 + __log_gamma(__y)
 137                 - __log_gamma(__x + __y);
 138 #endif
 139       __bet = std::exp(__bet);
 140       return __bet;
 141     }
 142
 143
 144     /**
 145      *   @brief  Return the beta function \f$B(x,y)\f$ using
 146      *           the product form.
 147      *
 148      *   The beta function is defined by
 149      *   @f[
 150      *     B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
 151      *   @f]
 152      *
 153      *   @param __x The first argument of the beta function.
 154      *   @param __y The second argument of the beta function.
 155      *   @return  The beta function.
 156      */
 157     template<typename _Tp>
 158     _Tp
 159     __beta_product(_Tp __x, _Tp __y)
 160     {
 161
 162       _Tp __bet = (__x + __y) / (__x * __y);
 163
 164       unsigned int __max_iter = 1000000;
 165       for (unsigned int __k = 1; __k < __max_iter; ++__k)
 166         {
 167           _Tp __term = (_Tp(1) + (__x + __y) / __k)
 168                      / ((_Tp(1) + __x / __k) * (_Tp(1) + __y / __k));
 169           __bet *= __term;
 170         }
 171
 172       return __bet;
 173     }
 174
 175
 176     /**
 177      *   @brief  Return the beta function \f$ B(x,y) \f$.
 178      *
 179      *   The beta function is defined by
 180      *   @f[
 181      *     B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}
 182      *   @f]
 183      *
 184      *   @param __x The first argument of the beta function.
 185      *   @param __y The second argument of the beta function.
 186      *   @return  The beta function.
 187      */
 188     template<typename _Tp>
 189     inline _Tp
 190     __beta(_Tp __x, _Tp __y)
 191     {
 192       if (__isnan(__x) || __isnan(__y))
 193         return std::numeric_limits<_Tp>::quiet_NaN();
 194       else
 195         return __beta_lgamma(__x, __y);
 196     }
 197
 198   } // namespace std::tr1::__detail
 199 }
 200 }
 201
 202 #endif // __GLIBCXX_TR1_BETA_FUNCTION_TCC