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

   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  *  Do not attempt to use it directly. @headername{tr1/cmath}
  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
  50 namespace std
  51 {
  52 namespace 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