Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | RSS feed

  1. // Special functions -*- C++ -*-
  2.  
  3. // Copyright (C) 2006-2013 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
  198.