Subversion Repositories Kolibri OS

Rev

Blame | Last modification | View Log | RSS feed

  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/hypergeometric.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:
  35. //   (1) Handbook of Mathematical Functions,
  36. //       ed. Milton Abramowitz and Irene A. Stegun,
  37. //       Dover Publications,
  38. //       Section 6, pp. 555-566
  39. //   (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl
  40.  
  41. #ifndef _GLIBCXX_TR1_HYPERGEOMETRIC_TCC
  42. #define _GLIBCXX_TR1_HYPERGEOMETRIC_TCC 1
  43.  
  44. namespace std _GLIBCXX_VISIBILITY(default)
  45. {
  46. namespace tr1
  47. {
  48.   // [5.2] Special functions
  49.  
  50.   // Implementation-space details.
  51.   namespace __detail
  52.   {
  53.   _GLIBCXX_BEGIN_NAMESPACE_VERSION
  54.  
  55.     /**
  56.      *   @brief This routine returns the confluent hypergeometric function
  57.      *          by series expansion.
  58.      *
  59.      *   @f[
  60.      *     _1F_1(a;c;x) = \frac{\Gamma(c)}{\Gamma(a)}
  61.      *                      \sum_{n=0}^{\infty}
  62.      *                      \frac{\Gamma(a+n)}{\Gamma(c+n)}
  63.      *                      \frac{x^n}{n!}
  64.      *   @f]
  65.      *
  66.      *   If a and b are integers and a < 0 and either b > 0 or b < a
  67.      *   then the series is a polynomial with a finite number of
  68.      *   terms.  If b is an integer and b <= 0 the confluent
  69.      *   hypergeometric function is undefined.
  70.      *
  71.      *   @param  __a  The "numerator" parameter.
  72.      *   @param  __c  The "denominator" parameter.
  73.      *   @param  __x  The argument of the confluent hypergeometric function.
  74.      *   @return  The confluent hypergeometric function.
  75.      */
  76.     template<typename _Tp>
  77.     _Tp
  78.     __conf_hyperg_series(_Tp __a, _Tp __c, _Tp __x)
  79.     {
  80.       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
  81.  
  82.       _Tp __term = _Tp(1);
  83.       _Tp __Fac = _Tp(1);
  84.       const unsigned int __max_iter = 100000;
  85.       unsigned int __i;
  86.       for (__i = 0; __i < __max_iter; ++__i)
  87.         {
  88.           __term *= (__a + _Tp(__i)) * __x
  89.                   / ((__c + _Tp(__i)) * _Tp(1 + __i));
  90.           if (std::abs(__term) < __eps)
  91.             {
  92.               break;
  93.             }
  94.           __Fac += __term;
  95.         }
  96.       if (__i == __max_iter)
  97.         std::__throw_runtime_error(__N("Series failed to converge "
  98.                                        "in __conf_hyperg_series."));
  99.  
  100.       return __Fac;
  101.     }
  102.  
  103.  
  104.     /**
  105.      *  @brief  Return the hypogeometric function @f$ _2F_1(a,b;c;x) @f$
  106.      *          by an iterative procedure described in
  107.      *          Luke, Algorithms for the Computation of Mathematical Functions.
  108.      *
  109.      *  Like the case of the 2F1 rational approximations, these are
  110.      *  probably guaranteed to converge for x < 0, barring gross    
  111.      *  numerical instability in the pre-asymptotic regime.        
  112.      */
  113.     template<typename _Tp>
  114.     _Tp
  115.     __conf_hyperg_luke(_Tp __a, _Tp __c, _Tp __xin)
  116.     {
  117.       const _Tp __big = std::pow(std::numeric_limits<_Tp>::max(), _Tp(0.16L));
  118.       const int __nmax = 20000;
  119.       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
  120.       const _Tp __x  = -__xin;
  121.       const _Tp __x3 = __x * __x * __x;
  122.       const _Tp __t0 = __a / __c;
  123.       const _Tp __t1 = (__a + _Tp(1)) / (_Tp(2) * __c);
  124.       const _Tp __t2 = (__a + _Tp(2)) / (_Tp(2) * (__c + _Tp(1)));
  125.       _Tp __F = _Tp(1);
  126.       _Tp __prec;
  127.  
  128.       _Tp __Bnm3 = _Tp(1);
  129.       _Tp __Bnm2 = _Tp(1) + __t1 * __x;
  130.       _Tp __Bnm1 = _Tp(1) + __t2 * __x * (_Tp(1) + __t1 / _Tp(3) * __x);
  131.  
  132.       _Tp __Anm3 = _Tp(1);
  133.       _Tp __Anm2 = __Bnm2 - __t0 * __x;
  134.       _Tp __Anm1 = __Bnm1 - __t0 * (_Tp(1) + __t2 * __x) * __x
  135.                  + __t0 * __t1 * (__c / (__c + _Tp(1))) * __x * __x;
  136.  
  137.       int __n = 3;
  138.       while(1)
  139.         {
  140.           _Tp __npam1 = _Tp(__n - 1) + __a;
  141.           _Tp __npcm1 = _Tp(__n - 1) + __c;
  142.           _Tp __npam2 = _Tp(__n - 2) + __a;
  143.           _Tp __npcm2 = _Tp(__n - 2) + __c;
  144.           _Tp __tnm1  = _Tp(2 * __n - 1);
  145.           _Tp __tnm3  = _Tp(2 * __n - 3);
  146.           _Tp __tnm5  = _Tp(2 * __n - 5);
  147.           _Tp __F1 =  (_Tp(__n - 2) - __a) / (_Tp(2) * __tnm3 * __npcm1);
  148.           _Tp __F2 =  (_Tp(__n) + __a) * __npam1
  149.                    / (_Tp(4) * __tnm1 * __tnm3 * __npcm2 * __npcm1);
  150.           _Tp __F3 = -__npam2 * __npam1 * (_Tp(__n - 2) - __a)
  151.                    / (_Tp(8) * __tnm3 * __tnm3 * __tnm5
  152.                    * (_Tp(__n - 3) + __c) * __npcm2 * __npcm1);
  153.           _Tp __E  = -__npam1 * (_Tp(__n - 1) - __c)
  154.                    / (_Tp(2) * __tnm3 * __npcm2 * __npcm1);
  155.  
  156.           _Tp __An = (_Tp(1) + __F1 * __x) * __Anm1
  157.                    + (__E + __F2 * __x) * __x * __Anm2 + __F3 * __x3 * __Anm3;
  158.           _Tp __Bn = (_Tp(1) + __F1 * __x) * __Bnm1
  159.                    + (__E + __F2 * __x) * __x * __Bnm2 + __F3 * __x3 * __Bnm3;
  160.           _Tp __r = __An / __Bn;
  161.  
  162.           __prec = std::abs((__F - __r) / __F);
  163.           __F = __r;
  164.  
  165.           if (__prec < __eps || __n > __nmax)
  166.             break;
  167.  
  168.           if (std::abs(__An) > __big || std::abs(__Bn) > __big)
  169.             {
  170.               __An   /= __big;
  171.               __Bn   /= __big;
  172.               __Anm1 /= __big;
  173.               __Bnm1 /= __big;
  174.               __Anm2 /= __big;
  175.               __Bnm2 /= __big;
  176.               __Anm3 /= __big;
  177.               __Bnm3 /= __big;
  178.             }
  179.           else if (std::abs(__An) < _Tp(1) / __big
  180.                 || std::abs(__Bn) < _Tp(1) / __big)
  181.             {
  182.               __An   *= __big;
  183.               __Bn   *= __big;
  184.               __Anm1 *= __big;
  185.               __Bnm1 *= __big;
  186.               __Anm2 *= __big;
  187.               __Bnm2 *= __big;
  188.               __Anm3 *= __big;
  189.               __Bnm3 *= __big;
  190.             }
  191.  
  192.           ++__n;
  193.           __Bnm3 = __Bnm2;
  194.           __Bnm2 = __Bnm1;
  195.           __Bnm1 = __Bn;
  196.           __Anm3 = __Anm2;
  197.           __Anm2 = __Anm1;
  198.           __Anm1 = __An;
  199.         }
  200.  
  201.       if (__n >= __nmax)
  202.         std::__throw_runtime_error(__N("Iteration failed to converge "
  203.                                        "in __conf_hyperg_luke."));
  204.  
  205.       return __F;
  206.     }
  207.  
  208.  
  209.     /**
  210.      *   @brief  Return the confluent hypogeometric function
  211.      *           @f$ _1F_1(a;c;x) @f$.
  212.      *
  213.      *   @todo  Handle b == nonpositive integer blowup - return NaN.
  214.      *
  215.      *   @param  __a  The @a numerator parameter.
  216.      *   @param  __c  The @a denominator parameter.
  217.      *   @param  __x  The argument of the confluent hypergeometric function.
  218.      *   @return  The confluent hypergeometric function.
  219.      */
  220.     template<typename _Tp>
  221.     _Tp
  222.     __conf_hyperg(_Tp __a, _Tp __c, _Tp __x)
  223.     {
  224. #if _GLIBCXX_USE_C99_MATH_TR1
  225.       const _Tp __c_nint = std::tr1::nearbyint(__c);
  226. #else
  227.       const _Tp __c_nint = static_cast<int>(__c + _Tp(0.5L));
  228. #endif
  229.       if (__isnan(__a) || __isnan(__c) || __isnan(__x))
  230.         return std::numeric_limits<_Tp>::quiet_NaN();
  231.       else if (__c_nint == __c && __c_nint <= 0)
  232.         return std::numeric_limits<_Tp>::infinity();
  233.       else if (__a == _Tp(0))
  234.         return _Tp(1);
  235.       else if (__c == __a)
  236.         return std::exp(__x);
  237.       else if (__x < _Tp(0))
  238.         return __conf_hyperg_luke(__a, __c, __x);
  239.       else
  240.         return __conf_hyperg_series(__a, __c, __x);
  241.     }
  242.  
  243.  
  244.     /**
  245.      *   @brief Return the hypogeometric function @f$ _2F_1(a,b;c;x) @f$
  246.      *   by series expansion.
  247.      *
  248.      *   The hypogeometric function is defined by
  249.      *   @f[
  250.      *     _2F_1(a,b;c;x) = \frac{\Gamma(c)}{\Gamma(a)\Gamma(b)}
  251.      *                      \sum_{n=0}^{\infty}
  252.      *                      \frac{\Gamma(a+n)\Gamma(b+n)}{\Gamma(c+n)}
  253.      *                      \frac{x^n}{n!}
  254.      *   @f]
  255.      *
  256.      *   This works and it's pretty fast.
  257.      *
  258.      *   @param  __a  The first @a numerator parameter.
  259.      *   @param  __a  The second @a numerator parameter.
  260.      *   @param  __c  The @a denominator parameter.
  261.      *   @param  __x  The argument of the confluent hypergeometric function.
  262.      *   @return  The confluent hypergeometric function.
  263.      */
  264.     template<typename _Tp>
  265.     _Tp
  266.     __hyperg_series(_Tp __a, _Tp __b, _Tp __c, _Tp __x)
  267.     {
  268.       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
  269.  
  270.       _Tp __term = _Tp(1);
  271.       _Tp __Fabc = _Tp(1);
  272.       const unsigned int __max_iter = 100000;
  273.       unsigned int __i;
  274.       for (__i = 0; __i < __max_iter; ++__i)
  275.         {
  276.           __term *= (__a + _Tp(__i)) * (__b + _Tp(__i)) * __x
  277.                   / ((__c + _Tp(__i)) * _Tp(1 + __i));
  278.           if (std::abs(__term) < __eps)
  279.             {
  280.               break;
  281.             }
  282.           __Fabc += __term;
  283.         }
  284.       if (__i == __max_iter)
  285.         std::__throw_runtime_error(__N("Series failed to converge "
  286.                                        "in __hyperg_series."));
  287.  
  288.       return __Fabc;
  289.     }
  290.  
  291.  
  292.     /**
  293.      *   @brief  Return the hypogeometric function @f$ _2F_1(a,b;c;x) @f$
  294.      *           by an iterative procedure described in
  295.      *           Luke, Algorithms for the Computation of Mathematical Functions.
  296.      */
  297.     template<typename _Tp>
  298.     _Tp
  299.     __hyperg_luke(_Tp __a, _Tp __b, _Tp __c, _Tp __xin)
  300.     {
  301.       const _Tp __big = std::pow(std::numeric_limits<_Tp>::max(), _Tp(0.16L));
  302.       const int __nmax = 20000;
  303.       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
  304.       const _Tp __x  = -__xin;
  305.       const _Tp __x3 = __x * __x * __x;
  306.       const _Tp __t0 = __a * __b / __c;
  307.       const _Tp __t1 = (__a + _Tp(1)) * (__b + _Tp(1)) / (_Tp(2) * __c);
  308.       const _Tp __t2 = (__a + _Tp(2)) * (__b + _Tp(2))
  309.                      / (_Tp(2) * (__c + _Tp(1)));
  310.  
  311.       _Tp __F = _Tp(1);
  312.  
  313.       _Tp __Bnm3 = _Tp(1);
  314.       _Tp __Bnm2 = _Tp(1) + __t1 * __x;
  315.       _Tp __Bnm1 = _Tp(1) + __t2 * __x * (_Tp(1) + __t1 / _Tp(3) * __x);
  316.  
  317.       _Tp __Anm3 = _Tp(1);
  318.       _Tp __Anm2 = __Bnm2 - __t0 * __x;
  319.       _Tp __Anm1 = __Bnm1 - __t0 * (_Tp(1) + __t2 * __x) * __x
  320.                  + __t0 * __t1 * (__c / (__c + _Tp(1))) * __x * __x;
  321.  
  322.       int __n = 3;
  323.       while (1)
  324.         {
  325.           const _Tp __npam1 = _Tp(__n - 1) + __a;
  326.           const _Tp __npbm1 = _Tp(__n - 1) + __b;
  327.           const _Tp __npcm1 = _Tp(__n - 1) + __c;
  328.           const _Tp __npam2 = _Tp(__n - 2) + __a;
  329.           const _Tp __npbm2 = _Tp(__n - 2) + __b;
  330.           const _Tp __npcm2 = _Tp(__n - 2) + __c;
  331.           const _Tp __tnm1  = _Tp(2 * __n - 1);
  332.           const _Tp __tnm3  = _Tp(2 * __n - 3);
  333.           const _Tp __tnm5  = _Tp(2 * __n - 5);
  334.           const _Tp __n2 = __n * __n;
  335.           const _Tp __F1 = (_Tp(3) * __n2 + (__a + __b - _Tp(6)) * __n
  336.                          + _Tp(2) - __a * __b - _Tp(2) * (__a + __b))
  337.                          / (_Tp(2) * __tnm3 * __npcm1);
  338.           const _Tp __F2 = -(_Tp(3) * __n2 - (__a + __b + _Tp(6)) * __n
  339.                          + _Tp(2) - __a * __b) * __npam1 * __npbm1
  340.                          / (_Tp(4) * __tnm1 * __tnm3 * __npcm2 * __npcm1);
  341.           const _Tp __F3 = (__npam2 * __npam1 * __npbm2 * __npbm1
  342.                          * (_Tp(__n - 2) - __a) * (_Tp(__n - 2) - __b))
  343.                          / (_Tp(8) * __tnm3 * __tnm3 * __tnm5
  344.                          * (_Tp(__n - 3) + __c) * __npcm2 * __npcm1);
  345.           const _Tp __E  = -__npam1 * __npbm1 * (_Tp(__n - 1) - __c)
  346.                          / (_Tp(2) * __tnm3 * __npcm2 * __npcm1);
  347.  
  348.           _Tp __An = (_Tp(1) + __F1 * __x) * __Anm1
  349.                    + (__E + __F2 * __x) * __x * __Anm2 + __F3 * __x3 * __Anm3;
  350.           _Tp __Bn = (_Tp(1) + __F1 * __x) * __Bnm1
  351.                    + (__E + __F2 * __x) * __x * __Bnm2 + __F3 * __x3 * __Bnm3;
  352.           const _Tp __r = __An / __Bn;
  353.  
  354.           const _Tp __prec = std::abs((__F - __r) / __F);
  355.           __F = __r;
  356.  
  357.           if (__prec < __eps || __n > __nmax)
  358.             break;
  359.  
  360.           if (std::abs(__An) > __big || std::abs(__Bn) > __big)
  361.             {
  362.               __An   /= __big;
  363.               __Bn   /= __big;
  364.               __Anm1 /= __big;
  365.               __Bnm1 /= __big;
  366.               __Anm2 /= __big;
  367.               __Bnm2 /= __big;
  368.               __Anm3 /= __big;
  369.               __Bnm3 /= __big;
  370.             }
  371.           else if (std::abs(__An) < _Tp(1) / __big
  372.                 || std::abs(__Bn) < _Tp(1) / __big)
  373.             {
  374.               __An   *= __big;
  375.               __Bn   *= __big;
  376.               __Anm1 *= __big;
  377.               __Bnm1 *= __big;
  378.               __Anm2 *= __big;
  379.               __Bnm2 *= __big;
  380.               __Anm3 *= __big;
  381.               __Bnm3 *= __big;
  382.             }
  383.  
  384.           ++__n;
  385.           __Bnm3 = __Bnm2;
  386.           __Bnm2 = __Bnm1;
  387.           __Bnm1 = __Bn;
  388.           __Anm3 = __Anm2;
  389.           __Anm2 = __Anm1;
  390.           __Anm1 = __An;
  391.         }
  392.  
  393.       if (__n >= __nmax)
  394.         std::__throw_runtime_error(__N("Iteration failed to converge "
  395.                                        "in __hyperg_luke."));
  396.  
  397.       return __F;
  398.     }
  399.  
  400.  
  401.     /**
  402.      *  @brief  Return the hypogeometric function @f$ _2F_1(a,b;c;x) @f$
  403.      *  by the reflection formulae in Abramowitz & Stegun formula
  404.      *  15.3.6 for d = c - a - b not integral and formula 15.3.11 for
  405.      *  d = c - a - b integral.  This assumes a, b, c != negative
  406.      *  integer.
  407.      *
  408.      *   The hypogeometric function is defined by
  409.      *   @f[
  410.      *     _2F_1(a,b;c;x) = \frac{\Gamma(c)}{\Gamma(a)\Gamma(b)}
  411.      *                      \sum_{n=0}^{\infty}
  412.      *                      \frac{\Gamma(a+n)\Gamma(b+n)}{\Gamma(c+n)}
  413.      *                      \frac{x^n}{n!}
  414.      *   @f]
  415.      *
  416.      *   The reflection formula for nonintegral @f$ d = c - a - b @f$ is:
  417.      *   @f[
  418.      *     _2F_1(a,b;c;x) = \frac{\Gamma(c)\Gamma(d)}{\Gamma(c-a)\Gamma(c-b)}
  419.      *                            _2F_1(a,b;1-d;1-x)
  420.      *                    + \frac{\Gamma(c)\Gamma(-d)}{\Gamma(a)\Gamma(b)}
  421.      *                            _2F_1(c-a,c-b;1+d;1-x)
  422.      *   @f]
  423.      *
  424.      *   The reflection formula for integral @f$ m = c - a - b @f$ is:
  425.      *   @f[
  426.      *     _2F_1(a,b;a+b+m;x) = \frac{\Gamma(m)\Gamma(a+b+m)}{\Gamma(a+m)\Gamma(b+m)}
  427.      *                        \sum_{k=0}^{m-1} \frac{(m+a)_k(m+b)_k}{k!(1-m)_k}
  428.      *                      -
  429.      *   @f]
  430.      */
  431.     template<typename _Tp>
  432.     _Tp
  433.     __hyperg_reflect(_Tp __a, _Tp __b, _Tp __c, _Tp __x)
  434.     {
  435.       const _Tp __d = __c - __a - __b;
  436.       const int __intd  = std::floor(__d + _Tp(0.5L));
  437.       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();
  438.       const _Tp __toler = _Tp(1000) * __eps;
  439.       const _Tp __log_max = std::log(std::numeric_limits<_Tp>::max());
  440.       const bool __d_integer = (std::abs(__d - __intd) < __toler);
  441.  
  442.       if (__d_integer)
  443.         {
  444.           const _Tp __ln_omx = std::log(_Tp(1) - __x);
  445.           const _Tp __ad = std::abs(__d);
  446.           _Tp __F1, __F2;
  447.  
  448.           _Tp __d1, __d2;
  449.           if (__d >= _Tp(0))
  450.             {
  451.               __d1 = __d;
  452.               __d2 = _Tp(0);
  453.             }
  454.           else
  455.             {
  456.               __d1 = _Tp(0);
  457.               __d2 = __d;
  458.             }
  459.  
  460.           const _Tp __lng_c = __log_gamma(__c);
  461.  
  462.           //  Evaluate F1.
  463.           if (__ad < __eps)
  464.             {
  465.               //  d = c - a - b = 0.
  466.               __F1 = _Tp(0);
  467.             }
  468.           else
  469.             {
  470.  
  471.               bool __ok_d1 = true;
  472.               _Tp __lng_ad, __lng_ad1, __lng_bd1;
  473.               __try
  474.                 {
  475.                   __lng_ad = __log_gamma(__ad);
  476.                   __lng_ad1 = __log_gamma(__a + __d1);
  477.                   __lng_bd1 = __log_gamma(__b + __d1);
  478.                 }
  479.               __catch(...)
  480.                 {
  481.                   __ok_d1 = false;
  482.                 }
  483.  
  484.               if (__ok_d1)
  485.                 {
  486.                   /* Gamma functions in the denominator are ok.
  487.                    * Proceed with evaluation.
  488.                    */
  489.                   _Tp __sum1 = _Tp(1);
  490.                   _Tp __term = _Tp(1);
  491.                   _Tp __ln_pre1 = __lng_ad + __lng_c + __d2 * __ln_omx
  492.                                 - __lng_ad1 - __lng_bd1;
  493.  
  494.                   /* Do F1 sum.
  495.                    */
  496.                   for (int __i = 1; __i < __ad; ++__i)
  497.                     {
  498.                       const int __j = __i - 1;
  499.                       __term *= (__a + __d2 + __j) * (__b + __d2 + __j)
  500.                               / (_Tp(1) + __d2 + __j) / __i * (_Tp(1) - __x);
  501.                       __sum1 += __term;
  502.                     }
  503.  
  504.                   if (__ln_pre1 > __log_max)
  505.                     std::__throw_runtime_error(__N("Overflow of gamma functions"
  506.                                                    " in __hyperg_luke."));
  507.                   else
  508.                     __F1 = std::exp(__ln_pre1) * __sum1;
  509.                 }
  510.               else
  511.                 {
  512.                   //  Gamma functions in the denominator were not ok.
  513.                   //  So the F1 term is zero.
  514.                   __F1 = _Tp(0);
  515.                 }
  516.             } // end F1 evaluation
  517.  
  518.           // Evaluate F2.
  519.           bool __ok_d2 = true;
  520.           _Tp __lng_ad2, __lng_bd2;
  521.           __try
  522.             {
  523.               __lng_ad2 = __log_gamma(__a + __d2);
  524.               __lng_bd2 = __log_gamma(__b + __d2);
  525.             }
  526.           __catch(...)
  527.             {
  528.               __ok_d2 = false;
  529.             }
  530.  
  531.           if (__ok_d2)
  532.             {
  533.               //  Gamma functions in the denominator are ok.
  534.               //  Proceed with evaluation.
  535.               const int __maxiter = 2000;
  536.               const _Tp __psi_1 = -__numeric_constants<_Tp>::__gamma_e();
  537.               const _Tp __psi_1pd = __psi(_Tp(1) + __ad);
  538.               const _Tp __psi_apd1 = __psi(__a + __d1);
  539.               const _Tp __psi_bpd1 = __psi(__b + __d1);
  540.  
  541.               _Tp __psi_term = __psi_1 + __psi_1pd - __psi_apd1
  542.                              - __psi_bpd1 - __ln_omx;
  543.               _Tp __fact = _Tp(1);
  544.               _Tp __sum2 = __psi_term;
  545.               _Tp __ln_pre2 = __lng_c + __d1 * __ln_omx
  546.                             - __lng_ad2 - __lng_bd2;
  547.  
  548.               // Do F2 sum.
  549.               int __j;
  550.               for (__j = 1; __j < __maxiter; ++__j)
  551.                 {
  552.                   //  Values for psi functions use recurrence;
  553.                   //  Abramowitz & Stegun 6.3.5
  554.                   const _Tp __term1 = _Tp(1) / _Tp(__j)
  555.                                     + _Tp(1) / (__ad + __j);
  556.                   const _Tp __term2 = _Tp(1) / (__a + __d1 + _Tp(__j - 1))
  557.                                     + _Tp(1) / (__b + __d1 + _Tp(__j - 1));
  558.                   __psi_term += __term1 - __term2;
  559.                   __fact *= (__a + __d1 + _Tp(__j - 1))
  560.                           * (__b + __d1 + _Tp(__j - 1))
  561.                           / ((__ad + __j) * __j) * (_Tp(1) - __x);
  562.                   const _Tp __delta = __fact * __psi_term;
  563.                   __sum2 += __delta;
  564.                   if (std::abs(__delta) < __eps * std::abs(__sum2))
  565.                     break;
  566.                 }
  567.               if (__j == __maxiter)
  568.                 std::__throw_runtime_error(__N("Sum F2 failed to converge "
  569.                                                "in __hyperg_reflect"));
  570.  
  571.               if (__sum2 == _Tp(0))
  572.                 __F2 = _Tp(0);
  573.               else
  574.                 __F2 = std::exp(__ln_pre2) * __sum2;
  575.             }
  576.           else
  577.             {
  578.               // Gamma functions in the denominator not ok.
  579.               // So the F2 term is zero.
  580.               __F2 = _Tp(0);
  581.             } // end F2 evaluation
  582.  
  583.           const _Tp __sgn_2 = (__intd % 2 == 1 ? -_Tp(1) : _Tp(1));
  584.           const _Tp __F = __F1 + __sgn_2 * __F2;
  585.  
  586.           return __F;
  587.         }
  588.       else
  589.         {
  590.           //  d = c - a - b not an integer.
  591.  
  592.           //  These gamma functions appear in the denominator, so we
  593.           //  catch their harmless domain errors and set the terms to zero.
  594.           bool __ok1 = true;
  595.           _Tp __sgn_g1ca = _Tp(0), __ln_g1ca = _Tp(0);
  596.           _Tp __sgn_g1cb = _Tp(0), __ln_g1cb = _Tp(0);
  597.           __try
  598.             {
  599.               __sgn_g1ca = __log_gamma_sign(__c - __a);
  600.               __ln_g1ca = __log_gamma(__c - __a);
  601.               __sgn_g1cb = __log_gamma_sign(__c - __b);
  602.               __ln_g1cb = __log_gamma(__c - __b);
  603.             }
  604.           __catch(...)
  605.             {
  606.               __ok1 = false;
  607.             }
  608.  
  609.           bool __ok2 = true;
  610.           _Tp __sgn_g2a = _Tp(0), __ln_g2a = _Tp(0);
  611.           _Tp __sgn_g2b = _Tp(0), __ln_g2b = _Tp(0);
  612.           __try
  613.             {
  614.               __sgn_g2a = __log_gamma_sign(__a);
  615.               __ln_g2a = __log_gamma(__a);
  616.               __sgn_g2b = __log_gamma_sign(__b);
  617.               __ln_g2b = __log_gamma(__b);
  618.             }
  619.           __catch(...)
  620.             {
  621.               __ok2 = false;
  622.             }
  623.  
  624.           const _Tp __sgn_gc = __log_gamma_sign(__c);
  625.           const _Tp __ln_gc = __log_gamma(__c);
  626.           const _Tp __sgn_gd = __log_gamma_sign(__d);
  627.           const _Tp __ln_gd = __log_gamma(__d);
  628.           const _Tp __sgn_gmd = __log_gamma_sign(-__d);
  629.           const _Tp __ln_gmd = __log_gamma(-__d);
  630.  
  631.           const _Tp __sgn1 = __sgn_gc * __sgn_gd  * __sgn_g1ca * __sgn_g1cb;
  632.           const _Tp __sgn2 = __sgn_gc * __sgn_gmd * __sgn_g2a  * __sgn_g2b;
  633.  
  634.           _Tp __pre1, __pre2;
  635.           if (__ok1 && __ok2)
  636.             {
  637.               _Tp __ln_pre1 = __ln_gc + __ln_gd  - __ln_g1ca - __ln_g1cb;
  638.               _Tp __ln_pre2 = __ln_gc + __ln_gmd - __ln_g2a  - __ln_g2b
  639.                             + __d * std::log(_Tp(1) - __x);
  640.               if (__ln_pre1 < __log_max && __ln_pre2 < __log_max)
  641.                 {
  642.                   __pre1 = std::exp(__ln_pre1);
  643.                   __pre2 = std::exp(__ln_pre2);
  644.                   __pre1 *= __sgn1;
  645.                   __pre2 *= __sgn2;
  646.                 }
  647.               else
  648.                 {
  649.                   std::__throw_runtime_error(__N("Overflow of gamma functions "
  650.                                                  "in __hyperg_reflect"));
  651.                 }
  652.             }
  653.           else if (__ok1 && !__ok2)
  654.             {
  655.               _Tp __ln_pre1 = __ln_gc + __ln_gd - __ln_g1ca - __ln_g1cb;
  656.               if (__ln_pre1 < __log_max)
  657.                 {
  658.                   __pre1 = std::exp(__ln_pre1);
  659.                   __pre1 *= __sgn1;
  660.                   __pre2 = _Tp(0);
  661.                 }
  662.               else
  663.                 {
  664.                   std::__throw_runtime_error(__N("Overflow of gamma functions "
  665.                                                  "in __hyperg_reflect"));
  666.                 }
  667.             }
  668.           else if (!__ok1 && __ok2)
  669.             {
  670.               _Tp __ln_pre2 = __ln_gc + __ln_gmd - __ln_g2a - __ln_g2b
  671.                             + __d * std::log(_Tp(1) - __x);
  672.               if (__ln_pre2 < __log_max)
  673.                 {
  674.                   __pre1 = _Tp(0);
  675.                   __pre2 = std::exp(__ln_pre2);
  676.                   __pre2 *= __sgn2;
  677.                 }
  678.               else
  679.                 {
  680.                   std::__throw_runtime_error(__N("Overflow of gamma functions "
  681.                                                  "in __hyperg_reflect"));
  682.                 }
  683.             }
  684.           else
  685.             {
  686.               __pre1 = _Tp(0);
  687.               __pre2 = _Tp(0);
  688.               std::__throw_runtime_error(__N("Underflow of gamma functions "
  689.                                              "in __hyperg_reflect"));
  690.             }
  691.  
  692.           const _Tp __F1 = __hyperg_series(__a, __b, _Tp(1) - __d,
  693.                                            _Tp(1) - __x);
  694.           const _Tp __F2 = __hyperg_series(__c - __a, __c - __b, _Tp(1) + __d,
  695.                                            _Tp(1) - __x);
  696.  
  697.           const _Tp __F = __pre1 * __F1 + __pre2 * __F2;
  698.  
  699.           return __F;
  700.         }
  701.     }
  702.  
  703.  
  704.     /**
  705.      *   @brief Return the hypogeometric function @f$ _2F_1(a,b;c;x) @f$.
  706.      *
  707.      *   The hypogeometric function is defined by
  708.      *   @f[
  709.      *     _2F_1(a,b;c;x) = \frac{\Gamma(c)}{\Gamma(a)\Gamma(b)}
  710.      *                      \sum_{n=0}^{\infty}
  711.      *                      \frac{\Gamma(a+n)\Gamma(b+n)}{\Gamma(c+n)}
  712.      *                      \frac{x^n}{n!}
  713.      *   @f]
  714.      *
  715.      *   @param  __a  The first @a numerator parameter.
  716.      *   @param  __a  The second @a numerator parameter.
  717.      *   @param  __c  The @a denominator parameter.
  718.      *   @param  __x  The argument of the confluent hypergeometric function.
  719.      *   @return  The confluent hypergeometric function.
  720.      */
  721.     template<typename _Tp>
  722.     _Tp
  723.     __hyperg(_Tp __a, _Tp __b, _Tp __c, _Tp __x)
  724.     {
  725. #if _GLIBCXX_USE_C99_MATH_TR1
  726.       const _Tp __a_nint = std::tr1::nearbyint(__a);
  727.       const _Tp __b_nint = std::tr1::nearbyint(__b);
  728.       const _Tp __c_nint = std::tr1::nearbyint(__c);
  729. #else
  730.       const _Tp __a_nint = static_cast<int>(__a + _Tp(0.5L));
  731.       const _Tp __b_nint = static_cast<int>(__b + _Tp(0.5L));
  732.       const _Tp __c_nint = static_cast<int>(__c + _Tp(0.5L));
  733. #endif
  734.       const _Tp __toler = _Tp(1000) * std::numeric_limits<_Tp>::epsilon();
  735.       if (std::abs(__x) >= _Tp(1))
  736.         std::__throw_domain_error(__N("Argument outside unit circle "
  737.                                       "in __hyperg."));
  738.       else if (__isnan(__a) || __isnan(__b)
  739.             || __isnan(__c) || __isnan(__x))
  740.         return std::numeric_limits<_Tp>::quiet_NaN();
  741.       else if (__c_nint == __c && __c_nint <= _Tp(0))
  742.         return std::numeric_limits<_Tp>::infinity();
  743.       else if (std::abs(__c - __b) < __toler || std::abs(__c - __a) < __toler)
  744.         return std::pow(_Tp(1) - __x, __c - __a - __b);
  745.       else if (__a >= _Tp(0) && __b >= _Tp(0) && __c >= _Tp(0)
  746.             && __x >= _Tp(0) && __x < _Tp(0.995L))
  747.         return __hyperg_series(__a, __b, __c, __x);
  748.       else if (std::abs(__a) < _Tp(10) && std::abs(__b) < _Tp(10))
  749.         {
  750.           //  For integer a and b the hypergeometric function is a
  751.           //  finite polynomial.
  752.           if (__a < _Tp(0)  &&  std::abs(__a - __a_nint) < __toler)
  753.             return __hyperg_series(__a_nint, __b, __c, __x);
  754.           else if (__b < _Tp(0)  &&  std::abs(__b - __b_nint) < __toler)
  755.             return __hyperg_series(__a, __b_nint, __c, __x);
  756.           else if (__x < -_Tp(0.25L))
  757.             return __hyperg_luke(__a, __b, __c, __x);
  758.           else if (__x < _Tp(0.5L))
  759.             return __hyperg_series(__a, __b, __c, __x);
  760.           else
  761.             if (std::abs(__c) > _Tp(10))
  762.               return __hyperg_series(__a, __b, __c, __x);
  763.             else
  764.               return __hyperg_reflect(__a, __b, __c, __x);
  765.         }
  766.       else
  767.         return __hyperg_luke(__a, __b, __c, __x);
  768.     }
  769.  
  770.   _GLIBCXX_END_NAMESPACE_VERSION
  771.   } // namespace std::tr1::__detail
  772. }
  773. }
  774.  
  775. #endif // _GLIBCXX_TR1_HYPERGEOMETRIC_TCC
  776.