visus
/
toolchain-gcc-arm-embedded
Watch 1
Star 0
Fork 0
Code Issues 0 Pull Requests 0 Releases 0 Wiki Activity
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
2 Commits
1 Branch
Branch: main
Branches Tags
${ item.name }
Create branch ${ searchTerm }
from 'main'
${ noResults }
toolchain-gcc-arm-embedded/arm-none-eabi/include/c++/10.2.1/tr1/hypergeometric.tcc

786 lines
27KB
Raw Permalink Blame History 

  1. // Special functions -*- C++ -*-

  2. 

  3. // Copyright (C) 2006-2020 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. _GLIBCXX_BEGIN_NAMESPACE_VERSION

  47. 

  48. #if _GLIBCXX_USE_STD_SPEC_FUNCS

  49. # define _GLIBCXX_MATH_NS ::std

  50. #elif defined(_GLIBCXX_TR1_CMATH)

  51. namespace tr1

  52. {

  53. # define _GLIBCXX_MATH_NS ::std::tr1

  54. #else

  55. # error do not include this header directly, use <cmath> or <tr1/cmath>

  56. #endif

  57.   // [5.2] Special functions

  58. 

  59.   // Implementation-space details.

  60.   namespace __detail

  61.   {

  62.     /**

  63.      *   @brief This routine returns the confluent hypergeometric function

  64.      *          by series expansion.

  65.      * 

  66.      *   @f[

  67.      *     _1F_1(a;c;x) = \frac{\Gamma(c)}{\Gamma(a)}

  68.      *                      \sum_{n=0}^{\infty}

  69.      *                      \frac{\Gamma(a+n)}{\Gamma(c+n)}

  70.      *                      \frac{x^n}{n!}

  71.      *   @f]

  72.      * 

  73.      *   If a and b are integers and a < 0 and either b > 0 or b < a

  74.      *   then the series is a polynomial with a finite number of

  75.      *   terms.  If b is an integer and b <= 0 the confluent

  76.      *   hypergeometric function is undefined.

  77.      *

  78.      *   @param  __a  The "numerator" parameter.

  79.      *   @param  __c  The "denominator" parameter.

  80.      *   @param  __x  The argument of the confluent hypergeometric function.

  81.      *   @return  The confluent hypergeometric function.

  82.      */

  83.     template<typename _Tp>

  84.     _Tp

  85.     __conf_hyperg_series(_Tp __a, _Tp __c, _Tp __x)

  86.     {

  87.       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();

  88. 

  89.       _Tp __term = _Tp(1);

  90.       _Tp __Fac = _Tp(1);

  91.       const unsigned int __max_iter = 100000;

  92.       unsigned int __i;

  93.       for (__i = 0; __i < __max_iter; ++__i)

  94.         {

  95.           __term *= (__a + _Tp(__i)) * __x

  96.                   / ((__c + _Tp(__i)) * _Tp(1 + __i));

  97.           if (std::abs(__term) < __eps)

  98.             {

  99.               break;

  100.             }

  101.           __Fac += __term;

  102.         }

  103.       if (__i == __max_iter)

  104.         std::__throw_runtime_error(__N("Series failed to converge "

  105.                                        "in __conf_hyperg_series."));

  106. 

  107.       return __Fac;

  108.     }

  109. 

  110. 

  111.     /**

  112.      *  @brief  Return the hypogeometric function @f$ _2F_1(a,b;c;x) @f$

  113.      *          by an iterative procedure described in

  114.      *          Luke, Algorithms for the Computation of Mathematical Functions.

  115.      *

  116.      *  Like the case of the 2F1 rational approximations, these are 

  117.      *  probably guaranteed to converge for x < 0, barring gross    

  118.      *  numerical instability in the pre-asymptotic regime.         

  119.      */

  120.     template<typename _Tp>

  121.     _Tp

  122.     __conf_hyperg_luke(_Tp __a, _Tp __c, _Tp __xin)

  123.     {

  124.       const _Tp __big = std::pow(std::numeric_limits<_Tp>::max(), _Tp(0.16L));

  125.       const int __nmax = 20000;

  126.       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();

  127.       const _Tp __x  = -__xin;

  128.       const _Tp __x3 = __x * __x * __x;

  129.       const _Tp __t0 = __a / __c;

  130.       const _Tp __t1 = (__a + _Tp(1)) / (_Tp(2) * __c);

  131.       const _Tp __t2 = (__a + _Tp(2)) / (_Tp(2) * (__c + _Tp(1)));

  132.       _Tp __F = _Tp(1);

  133.       _Tp __prec;

  134. 

  135.       _Tp __Bnm3 = _Tp(1);

  136.       _Tp __Bnm2 = _Tp(1) + __t1 * __x;

  137.       _Tp __Bnm1 = _Tp(1) + __t2 * __x * (_Tp(1) + __t1 / _Tp(3) * __x);

  138. 

  139.       _Tp __Anm3 = _Tp(1);

  140.       _Tp __Anm2 = __Bnm2 - __t0 * __x;

  141.       _Tp __Anm1 = __Bnm1 - __t0 * (_Tp(1) + __t2 * __x) * __x

  142.                  + __t0 * __t1 * (__c / (__c + _Tp(1))) * __x * __x;

  143. 

  144.       int __n = 3;

  145.       while(1)

  146.         {

  147.           _Tp __npam1 = _Tp(__n - 1) + __a;

  148.           _Tp __npcm1 = _Tp(__n - 1) + __c;

  149.           _Tp __npam2 = _Tp(__n - 2) + __a;

  150.           _Tp __npcm2 = _Tp(__n - 2) + __c;

  151.           _Tp __tnm1  = _Tp(2 * __n - 1);

  152.           _Tp __tnm3  = _Tp(2 * __n - 3);

  153.           _Tp __tnm5  = _Tp(2 * __n - 5);

  154.           _Tp __F1 =  (_Tp(__n - 2) - __a) / (_Tp(2) * __tnm3 * __npcm1);

  155.           _Tp __F2 =  (_Tp(__n) + __a) * __npam1

  156.                    / (_Tp(4) * __tnm1 * __tnm3 * __npcm2 * __npcm1);

  157.           _Tp __F3 = -__npam2 * __npam1 * (_Tp(__n - 2) - __a)

  158.                    / (_Tp(8) * __tnm3 * __tnm3 * __tnm5

  159.                    * (_Tp(__n - 3) + __c) * __npcm2 * __npcm1);

  160.           _Tp __E  = -__npam1 * (_Tp(__n - 1) - __c)

  161.                    / (_Tp(2) * __tnm3 * __npcm2 * __npcm1);

  162. 

  163.           _Tp __An = (_Tp(1) + __F1 * __x) * __Anm1

  164.                    + (__E + __F2 * __x) * __x * __Anm2 + __F3 * __x3 * __Anm3;

  165.           _Tp __Bn = (_Tp(1) + __F1 * __x) * __Bnm1

  166.                    + (__E + __F2 * __x) * __x * __Bnm2 + __F3 * __x3 * __Bnm3;

  167.           _Tp __r = __An / __Bn;

  168. 

  169.           __prec = std::abs((__F - __r) / __F);

  170.           __F = __r;

  171. 

  172.           if (__prec < __eps || __n > __nmax)

  173.             break;

  174. 

  175.           if (std::abs(__An) > __big || std::abs(__Bn) > __big)

  176.             {

  177.               __An   /= __big;

  178.               __Bn   /= __big;

  179.               __Anm1 /= __big;

  180.               __Bnm1 /= __big;

  181.               __Anm2 /= __big;

  182.               __Bnm2 /= __big;

  183.               __Anm3 /= __big;

  184.               __Bnm3 /= __big;

  185.             }

  186.           else if (std::abs(__An) < _Tp(1) / __big

  187.                 || std::abs(__Bn) < _Tp(1) / __big)

  188.             {

  189.               __An   *= __big;

  190.               __Bn   *= __big;

  191.               __Anm1 *= __big;

  192.               __Bnm1 *= __big;

  193.               __Anm2 *= __big;

  194.               __Bnm2 *= __big;

  195.               __Anm3 *= __big;

  196.               __Bnm3 *= __big;

  197.             }

  198. 

  199.           ++__n;

  200.           __Bnm3 = __Bnm2;

  201.           __Bnm2 = __Bnm1;

  202.           __Bnm1 = __Bn;

  203.           __Anm3 = __Anm2;

  204.           __Anm2 = __Anm1;

  205.           __Anm1 = __An;

  206.         }

  207. 

  208.       if (__n >= __nmax)

  209.         std::__throw_runtime_error(__N("Iteration failed to converge "

  210.                                        "in __conf_hyperg_luke."));

  211. 

  212.       return __F;

  213.     }

  214. 

  215. 

  216.     /**

  217.      *   @brief  Return the confluent hypogeometric function

  218.      *           @f$ _1F_1(a;c;x) @f$.

  219.      * 

  220.      *   @todo  Handle b == nonpositive integer blowup - return NaN.

  221.      *

  222.      *   @param  __a  The @a numerator parameter.

  223.      *   @param  __c  The @a denominator parameter.

  224.      *   @param  __x  The argument of the confluent hypergeometric function.

  225.      *   @return  The confluent hypergeometric function.

  226.      */

  227.     template<typename _Tp>

  228.     _Tp

  229.     __conf_hyperg(_Tp __a, _Tp __c, _Tp __x)

  230.     {

  231. #if _GLIBCXX_USE_C99_MATH_TR1

  232.       const _Tp __c_nint = _GLIBCXX_MATH_NS::nearbyint(__c);

  233. #else

  234.       const _Tp __c_nint = static_cast<int>(__c + _Tp(0.5L));

  235. #endif

  236.       if (__isnan(__a) || __isnan(__c) || __isnan(__x))

  237.         return std::numeric_limits<_Tp>::quiet_NaN();

  238.       else if (__c_nint == __c && __c_nint <= 0)

  239.         return std::numeric_limits<_Tp>::infinity();

  240.       else if (__a == _Tp(0))

  241.         return _Tp(1);

  242.       else if (__c == __a)

  243.         return std::exp(__x);

  244.       else if (__x < _Tp(0))

  245.         return __conf_hyperg_luke(__a, __c, __x);

  246.       else

  247.         return __conf_hyperg_series(__a, __c, __x);

  248.     }

  249. 

  250. 

  251.     /**

  252.      *   @brief Return the hypogeometric function @f$ _2F_1(a,b;c;x) @f$

  253.      *   by series expansion.

  254.      * 

  255.      *   The hypogeometric function is defined by

  256.      *   @f[

  257.      *     _2F_1(a,b;c;x) = \frac{\Gamma(c)}{\Gamma(a)\Gamma(b)}

  258.      *                      \sum_{n=0}^{\infty}

  259.      *                      \frac{\Gamma(a+n)\Gamma(b+n)}{\Gamma(c+n)}

  260.      *                      \frac{x^n}{n!}

  261.      *   @f]

  262.      * 

  263.      *   This works and it's pretty fast.

  264.      *

  265.      *   @param  __a  The first @a numerator parameter.

  266.      *   @param  __a  The second @a numerator parameter.

  267.      *   @param  __c  The @a denominator parameter.

  268.      *   @param  __x  The argument of the confluent hypergeometric function.

  269.      *   @return  The confluent hypergeometric function.

  270.      */

  271.     template<typename _Tp>

  272.     _Tp

  273.     __hyperg_series(_Tp __a, _Tp __b, _Tp __c, _Tp __x)

  274.     {

  275.       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();

  276. 

  277.       _Tp __term = _Tp(1);

  278.       _Tp __Fabc = _Tp(1);

  279.       const unsigned int __max_iter = 100000;

  280.       unsigned int __i;

  281.       for (__i = 0; __i < __max_iter; ++__i)

  282.         {

  283.           __term *= (__a + _Tp(__i)) * (__b + _Tp(__i)) * __x

  284.                   / ((__c + _Tp(__i)) * _Tp(1 + __i));

  285.           if (std::abs(__term) < __eps)

  286.             {

  287.               break;

  288.             }

  289.           __Fabc += __term;

  290.         }

  291.       if (__i == __max_iter)

  292.         std::__throw_runtime_error(__N("Series failed to converge "

  293.                                        "in __hyperg_series."));

  294. 

  295.       return __Fabc;

  296.     }

  297. 

  298. 

  299.     /**

  300.      *   @brief  Return the hypogeometric function @f$ _2F_1(a,b;c;x) @f$

  301.      *           by an iterative procedure described in

  302.      *           Luke, Algorithms for the Computation of Mathematical Functions.

  303.      */

  304.     template<typename _Tp>

  305.     _Tp

  306.     __hyperg_luke(_Tp __a, _Tp __b, _Tp __c, _Tp __xin)

  307.     {

  308.       const _Tp __big = std::pow(std::numeric_limits<_Tp>::max(), _Tp(0.16L));

  309.       const int __nmax = 20000;

  310.       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();

  311.       const _Tp __x  = -__xin;

  312.       const _Tp __x3 = __x * __x * __x;

  313.       const _Tp __t0 = __a * __b / __c;

  314.       const _Tp __t1 = (__a + _Tp(1)) * (__b + _Tp(1)) / (_Tp(2) * __c);

  315.       const _Tp __t2 = (__a + _Tp(2)) * (__b + _Tp(2))

  316.                      / (_Tp(2) * (__c + _Tp(1)));

  317. 

  318.       _Tp __F = _Tp(1);

  319. 

  320.       _Tp __Bnm3 = _Tp(1);

  321.       _Tp __Bnm2 = _Tp(1) + __t1 * __x;

  322.       _Tp __Bnm1 = _Tp(1) + __t2 * __x * (_Tp(1) + __t1 / _Tp(3) * __x);

  323. 

  324.       _Tp __Anm3 = _Tp(1);

  325.       _Tp __Anm2 = __Bnm2 - __t0 * __x;

  326.       _Tp __Anm1 = __Bnm1 - __t0 * (_Tp(1) + __t2 * __x) * __x

  327.                  + __t0 * __t1 * (__c / (__c + _Tp(1))) * __x * __x;

  328. 

  329.       int __n = 3;

  330.       while (1)

  331.         {

  332.           const _Tp __npam1 = _Tp(__n - 1) + __a;

  333.           const _Tp __npbm1 = _Tp(__n - 1) + __b;

  334.           const _Tp __npcm1 = _Tp(__n - 1) + __c;

  335.           const _Tp __npam2 = _Tp(__n - 2) + __a;

  336.           const _Tp __npbm2 = _Tp(__n - 2) + __b;

  337.           const _Tp __npcm2 = _Tp(__n - 2) + __c;

  338.           const _Tp __tnm1  = _Tp(2 * __n - 1);

  339.           const _Tp __tnm3  = _Tp(2 * __n - 3);

  340.           const _Tp __tnm5  = _Tp(2 * __n - 5);

  341.           const _Tp __n2 = __n * __n;

  342.           const _Tp __F1 = (_Tp(3) * __n2 + (__a + __b - _Tp(6)) * __n

  343.                          + _Tp(2) - __a * __b - _Tp(2) * (__a + __b))

  344.                          / (_Tp(2) * __tnm3 * __npcm1);

  345.           const _Tp __F2 = -(_Tp(3) * __n2 - (__a + __b + _Tp(6)) * __n

  346.                          + _Tp(2) - __a * __b) * __npam1 * __npbm1

  347.                          / (_Tp(4) * __tnm1 * __tnm3 * __npcm2 * __npcm1);

  348.           const _Tp __F3 = (__npam2 * __npam1 * __npbm2 * __npbm1

  349.                          * (_Tp(__n - 2) - __a) * (_Tp(__n - 2) - __b))

  350.                          / (_Tp(8) * __tnm3 * __tnm3 * __tnm5

  351.                          * (_Tp(__n - 3) + __c) * __npcm2 * __npcm1);

  352.           const _Tp __E  = -__npam1 * __npbm1 * (_Tp(__n - 1) - __c)

  353.                          / (_Tp(2) * __tnm3 * __npcm2 * __npcm1);

  354. 

  355.           _Tp __An = (_Tp(1) + __F1 * __x) * __Anm1

  356.                    + (__E + __F2 * __x) * __x * __Anm2 + __F3 * __x3 * __Anm3;

  357.           _Tp __Bn = (_Tp(1) + __F1 * __x) * __Bnm1

  358.                    + (__E + __F2 * __x) * __x * __Bnm2 + __F3 * __x3 * __Bnm3;

  359.           const _Tp __r = __An / __Bn;

  360. 

  361.           const _Tp __prec = std::abs((__F - __r) / __F);

  362.           __F = __r;

  363. 

  364.           if (__prec < __eps || __n > __nmax)

  365.             break;

  366. 

  367.           if (std::abs(__An) > __big || std::abs(__Bn) > __big)

  368.             {

  369.               __An   /= __big;

  370.               __Bn   /= __big;

  371.               __Anm1 /= __big;

  372.               __Bnm1 /= __big;

  373.               __Anm2 /= __big;

  374.               __Bnm2 /= __big;

  375.               __Anm3 /= __big;

  376.               __Bnm3 /= __big;

  377.             }

  378.           else if (std::abs(__An) < _Tp(1) / __big

  379.                 || std::abs(__Bn) < _Tp(1) / __big)

  380.             {

  381.               __An   *= __big;

  382.               __Bn   *= __big;

  383.               __Anm1 *= __big;

  384.               __Bnm1 *= __big;

  385.               __Anm2 *= __big;

  386.               __Bnm2 *= __big;

  387.               __Anm3 *= __big;

  388.               __Bnm3 *= __big;

  389.             }

  390. 

  391.           ++__n;

  392.           __Bnm3 = __Bnm2;

  393.           __Bnm2 = __Bnm1;

  394.           __Bnm1 = __Bn;

  395.           __Anm3 = __Anm2;

  396.           __Anm2 = __Anm1;

  397.           __Anm1 = __An;

  398.         }

  399. 

  400.       if (__n >= __nmax)

  401.         std::__throw_runtime_error(__N("Iteration failed to converge "

  402.                                        "in __hyperg_luke."));

  403. 

  404.       return __F;

  405.     }

  406. 

  407. 

  408.     /**

  409.      *  @brief  Return the hypogeometric function @f$ _2F_1(a,b;c;x) @f$ 

  410.      *  by the reflection formulae in Abramowitz & Stegun formula

  411.      *  15.3.6 for d = c - a - b not integral and formula 15.3.11 for

  412.      *  d = c - a - b integral.  This assumes a, b, c != negative

  413.      *  integer.

  414.      *

  415.      *   The hypogeometric function is defined by

  416.      *   @f[

  417.      *     _2F_1(a,b;c;x) = \frac{\Gamma(c)}{\Gamma(a)\Gamma(b)}

  418.      *                      \sum_{n=0}^{\infty}

  419.      *                      \frac{\Gamma(a+n)\Gamma(b+n)}{\Gamma(c+n)}

  420.      *                      \frac{x^n}{n!}

  421.      *   @f]

  422.      *

  423.      *   The reflection formula for nonintegral @f$ d = c - a - b @f$ is:

  424.      *   @f[

  425.      *     _2F_1(a,b;c;x) = \frac{\Gamma(c)\Gamma(d)}{\Gamma(c-a)\Gamma(c-b)}

  426.      *                            _2F_1(a,b;1-d;1-x)

  427.      *                    + \frac{\Gamma(c)\Gamma(-d)}{\Gamma(a)\Gamma(b)}

  428.      *                            _2F_1(c-a,c-b;1+d;1-x)

  429.      *   @f]

  430.      *

  431.      *   The reflection formula for integral @f$ m = c - a - b @f$ is:

  432.      *   @f[

  433.      *     _2F_1(a,b;a+b+m;x) = \frac{\Gamma(m)\Gamma(a+b+m)}{\Gamma(a+m)\Gamma(b+m)}

  434.      *                        \sum_{k=0}^{m-1} \frac{(m+a)_k(m+b)_k}{k!(1-m)_k}

  435.      *                      - 

  436.      *   @f]

  437.      */

  438.     template<typename _Tp>

  439.     _Tp

  440.     __hyperg_reflect(_Tp __a, _Tp __b, _Tp __c, _Tp __x)

  441.     {

  442.       const _Tp __d = __c - __a - __b;

  443.       const int __intd  = std::floor(__d + _Tp(0.5L));

  444.       const _Tp __eps = std::numeric_limits<_Tp>::epsilon();

  445.       const _Tp __toler = _Tp(1000) * __eps;

  446.       const _Tp __log_max = std::log(std::numeric_limits<_Tp>::max());

  447.       const bool __d_integer = (std::abs(__d - __intd) < __toler);

  448. 

  449.       if (__d_integer)

  450.         {

  451.           const _Tp __ln_omx = std::log(_Tp(1) - __x);

  452.           const _Tp __ad = std::abs(__d);

  453.           _Tp __F1, __F2;

  454. 

  455.           _Tp __d1, __d2;

  456.           if (__d >= _Tp(0))

  457.             {

  458.               __d1 = __d;

  459.               __d2 = _Tp(0);

  460.             }

  461.           else

  462.             {

  463.               __d1 = _Tp(0);

  464.               __d2 = __d;

  465.             }

  466. 

  467.           const _Tp __lng_c = __log_gamma(__c);

  468. 

  469.           //  Evaluate F1.

  470.           if (__ad < __eps)

  471.             {

  472.               //  d = c - a - b = 0.

  473.               __F1 = _Tp(0);

  474.             }

  475.           else

  476.             {

  477. 

  478.               bool __ok_d1 = true;

  479.               _Tp __lng_ad, __lng_ad1, __lng_bd1;

  480.               __try

  481.                 {

  482.                   __lng_ad = __log_gamma(__ad);

  483.                   __lng_ad1 = __log_gamma(__a + __d1);

  484.                   __lng_bd1 = __log_gamma(__b + __d1);

  485.                 }

  486.               __catch(...)

  487.                 {

  488.                   __ok_d1 = false;

  489.                 }

  490. 

  491.               if (__ok_d1)

  492.                 {

  493.                   /* Gamma functions in the denominator are ok.

  494.                    * Proceed with evaluation.

  495.                    */

  496.                   _Tp __sum1 = _Tp(1);

  497.                   _Tp __term = _Tp(1);

  498.                   _Tp __ln_pre1 = __lng_ad + __lng_c + __d2 * __ln_omx

  499.                                 - __lng_ad1 - __lng_bd1;

  500. 

  501.                   /* Do F1 sum.

  502.                    */

  503.                   for (int __i = 1; __i < __ad; ++__i)

  504.                     {

  505.                       const int __j = __i - 1;

  506.                       __term *= (__a + __d2 + __j) * (__b + __d2 + __j)

  507.                               / (_Tp(1) + __d2 + __j) / __i * (_Tp(1) - __x);

  508.                       __sum1 += __term;

  509.                     }

  510. 

  511.                   if (__ln_pre1 > __log_max)

  512.                     std::__throw_runtime_error(__N("Overflow of gamma functions"

  513.                                                    " in __hyperg_luke."));

  514.                   else

  515.                     __F1 = std::exp(__ln_pre1) * __sum1;

  516.                 }

  517.               else

  518.                 {

  519.                   //  Gamma functions in the denominator were not ok.

  520.                   //  So the F1 term is zero.

  521.                   __F1 = _Tp(0);

  522.                 }

  523.             } // end F1 evaluation

  524. 

  525.           // Evaluate F2.

  526.           bool __ok_d2 = true;

  527.           _Tp __lng_ad2, __lng_bd2;

  528.           __try

  529.             {

  530.               __lng_ad2 = __log_gamma(__a + __d2);

  531.               __lng_bd2 = __log_gamma(__b + __d2);

  532.             }

  533.           __catch(...)

  534.             {

  535.               __ok_d2 = false;

  536.             }

  537. 

  538.           if (__ok_d2)

  539.             {

  540.               //  Gamma functions in the denominator are ok.

  541.               //  Proceed with evaluation.

  542.               const int __maxiter = 2000;

  543.               const _Tp __psi_1 = -__numeric_constants<_Tp>::__gamma_e();

  544.               const _Tp __psi_1pd = __psi(_Tp(1) + __ad);

  545.               const _Tp __psi_apd1 = __psi(__a + __d1);

  546.               const _Tp __psi_bpd1 = __psi(__b + __d1);

  547. 

  548.               _Tp __psi_term = __psi_1 + __psi_1pd - __psi_apd1

  549.                              - __psi_bpd1 - __ln_omx;

  550.               _Tp __fact = _Tp(1);

  551.               _Tp __sum2 = __psi_term;

  552.               _Tp __ln_pre2 = __lng_c + __d1 * __ln_omx

  553.                             - __lng_ad2 - __lng_bd2;

  554. 

  555.               // Do F2 sum.

  556.               int __j;

  557.               for (__j = 1; __j < __maxiter; ++__j)

  558.                 {

  559.                   //  Values for psi functions use recurrence;

  560.                   //  Abramowitz & Stegun 6.3.5

  561.                   const _Tp __term1 = _Tp(1) / _Tp(__j)

  562.                                     + _Tp(1) / (__ad + __j);

  563.                   const _Tp __term2 = _Tp(1) / (__a + __d1 + _Tp(__j - 1))

  564.                                     + _Tp(1) / (__b + __d1 + _Tp(__j - 1));

  565.                   __psi_term += __term1 - __term2;

  566.                   __fact *= (__a + __d1 + _Tp(__j - 1))

  567.                           * (__b + __d1 + _Tp(__j - 1))

  568.                           / ((__ad + __j) * __j) * (_Tp(1) - __x);

  569.                   const _Tp __delta = __fact * __psi_term;

  570.                   __sum2 += __delta;

  571.                   if (std::abs(__delta) < __eps * std::abs(__sum2))

  572.                     break;

  573.                 }

  574.               if (__j == __maxiter)

  575.                 std::__throw_runtime_error(__N("Sum F2 failed to converge "

  576.                                                "in __hyperg_reflect"));

  577. 

  578.               if (__sum2 == _Tp(0))

  579.                 __F2 = _Tp(0);

  580.               else

  581.                 __F2 = std::exp(__ln_pre2) * __sum2;

  582.             }

  583.           else

  584.             {

  585.               // Gamma functions in the denominator not ok.

  586.               // So the F2 term is zero.

  587.               __F2 = _Tp(0);

  588.             } // end F2 evaluation

  589. 

  590.           const _Tp __sgn_2 = (__intd % 2 == 1 ? -_Tp(1) : _Tp(1));

  591.           const _Tp __F = __F1 + __sgn_2 * __F2;

  592. 

  593.           return __F;

  594.         }

  595.       else

  596.         {

  597.           //  d = c - a - b not an integer.

  598. 

  599.           //  These gamma functions appear in the denominator, so we

  600.           //  catch their harmless domain errors and set the terms to zero.

  601.           bool __ok1 = true;

  602.           _Tp __sgn_g1ca = _Tp(0), __ln_g1ca = _Tp(0);

  603.           _Tp __sgn_g1cb = _Tp(0), __ln_g1cb = _Tp(0);

  604.           __try

  605.             {

  606.               __sgn_g1ca = __log_gamma_sign(__c - __a);

  607.               __ln_g1ca = __log_gamma(__c - __a);

  608.               __sgn_g1cb = __log_gamma_sign(__c - __b);

  609.               __ln_g1cb = __log_gamma(__c - __b);

  610.             }

  611.           __catch(...)

  612.             {

  613.               __ok1 = false;

  614.             }

  615. 

  616.           bool __ok2 = true;

  617.           _Tp __sgn_g2a = _Tp(0), __ln_g2a = _Tp(0);

  618.           _Tp __sgn_g2b = _Tp(0), __ln_g2b = _Tp(0);

  619.           __try

  620.             {

  621.               __sgn_g2a = __log_gamma_sign(__a);

  622.               __ln_g2a = __log_gamma(__a);

  623.               __sgn_g2b = __log_gamma_sign(__b);

  624.               __ln_g2b = __log_gamma(__b);

  625.             }

  626.           __catch(...)

  627.             {

  628.               __ok2 = false;

  629.             }

  630. 

  631.           const _Tp __sgn_gc = __log_gamma_sign(__c);

  632.           const _Tp __ln_gc = __log_gamma(__c);

  633.           const _Tp __sgn_gd = __log_gamma_sign(__d);

  634.           const _Tp __ln_gd = __log_gamma(__d);

  635.           const _Tp __sgn_gmd = __log_gamma_sign(-__d);

  636.           const _Tp __ln_gmd = __log_gamma(-__d);

  637. 

  638.           const _Tp __sgn1 = __sgn_gc * __sgn_gd  * __sgn_g1ca * __sgn_g1cb;

  639.           const _Tp __sgn2 = __sgn_gc * __sgn_gmd * __sgn_g2a  * __sgn_g2b;

  640. 

  641.           _Tp __pre1, __pre2;

  642.           if (__ok1 && __ok2)

  643.             {

  644.               _Tp __ln_pre1 = __ln_gc + __ln_gd  - __ln_g1ca - __ln_g1cb;

  645.               _Tp __ln_pre2 = __ln_gc + __ln_gmd - __ln_g2a  - __ln_g2b

  646.                             + __d * std::log(_Tp(1) - __x);

  647.               if (__ln_pre1 < __log_max && __ln_pre2 < __log_max)

  648.                 {

  649.                   __pre1 = std::exp(__ln_pre1);

  650.                   __pre2 = std::exp(__ln_pre2);

  651.                   __pre1 *= __sgn1;

  652.                   __pre2 *= __sgn2;

  653.                 }

  654.               else

  655.                 {

  656.                   std::__throw_runtime_error(__N("Overflow of gamma functions "

  657.                                                  "in __hyperg_reflect"));

  658.                 }

  659.             }

  660.           else if (__ok1 && !__ok2)

  661.             {

  662.               _Tp __ln_pre1 = __ln_gc + __ln_gd - __ln_g1ca - __ln_g1cb;

  663.               if (__ln_pre1 < __log_max)

  664.                 {

  665.                   __pre1 = std::exp(__ln_pre1);

  666.                   __pre1 *= __sgn1;

  667.                   __pre2 = _Tp(0);

  668.                 }

  669.               else

  670.                 {

  671.                   std::__throw_runtime_error(__N("Overflow of gamma functions "

  672.                                                  "in __hyperg_reflect"));

  673.                 }

  674.             }

  675.           else if (!__ok1 && __ok2)

  676.             {

  677.               _Tp __ln_pre2 = __ln_gc + __ln_gmd - __ln_g2a - __ln_g2b

  678.                             + __d * std::log(_Tp(1) - __x);

  679.               if (__ln_pre2 < __log_max)

  680.                 {

  681.                   __pre1 = _Tp(0);

  682.                   __pre2 = std::exp(__ln_pre2);

  683.                   __pre2 *= __sgn2;

  684.                 }

  685.               else

  686.                 {

  687.                   std::__throw_runtime_error(__N("Overflow of gamma functions "

  688.                                                  "in __hyperg_reflect"));

  689.                 }

  690.             }

  691.           else

  692.             {

  693.               __pre1 = _Tp(0);

  694.               __pre2 = _Tp(0);

  695.               std::__throw_runtime_error(__N("Underflow of gamma functions "

  696.                                              "in __hyperg_reflect"));

  697.             }

  698. 

  699.           const _Tp __F1 = __hyperg_series(__a, __b, _Tp(1) - __d,

  700.                                            _Tp(1) - __x);

  701.           const _Tp __F2 = __hyperg_series(__c - __a, __c - __b, _Tp(1) + __d,

  702.                                            _Tp(1) - __x);

  703. 

  704.           const _Tp __F = __pre1 * __F1 + __pre2 * __F2;

  705. 

  706.           return __F;

  707.         }

  708.     }

  709. 

  710. 

  711.     /**

  712.      *   @brief Return the hypogeometric function @f$ _2F_1(a,b;c;x) @f$.

  713.      *

  714.      *   The hypogeometric function is defined by

  715.      *   @f[

  716.      *     _2F_1(a,b;c;x) = \frac{\Gamma(c)}{\Gamma(a)\Gamma(b)}

  717.      *                      \sum_{n=0}^{\infty}

  718.      *                      \frac{\Gamma(a+n)\Gamma(b+n)}{\Gamma(c+n)}

  719.      *                      \frac{x^n}{n!}

  720.      *   @f]

  721.      *

  722.      *   @param  __a  The first @a numerator parameter.

  723.      *   @param  __a  The second @a numerator parameter.

  724.      *   @param  __c  The @a denominator parameter.

  725.      *   @param  __x  The argument of the confluent hypergeometric function.

  726.      *   @return  The confluent hypergeometric function.

  727.      */

  728.     template<typename _Tp>

  729.     _Tp

  730.     __hyperg(_Tp __a, _Tp __b, _Tp __c, _Tp __x)

  731.     {

  732. #if _GLIBCXX_USE_C99_MATH_TR1

  733.       const _Tp __a_nint = _GLIBCXX_MATH_NS::nearbyint(__a);

  734.       const _Tp __b_nint = _GLIBCXX_MATH_NS::nearbyint(__b);

  735.       const _Tp __c_nint = _GLIBCXX_MATH_NS::nearbyint(__c);

  736. #else

  737.       const _Tp __a_nint = static_cast<int>(__a + _Tp(0.5L));

  738.       const _Tp __b_nint = static_cast<int>(__b + _Tp(0.5L));

  739.       const _Tp __c_nint = static_cast<int>(__c + _Tp(0.5L));

  740. #endif

  741.       const _Tp __toler = _Tp(1000) * std::numeric_limits<_Tp>::epsilon();

  742.       if (std::abs(__x) >= _Tp(1))

  743.         std::__throw_domain_error(__N("Argument outside unit circle "

  744.                                       "in __hyperg."));

  745.       else if (__isnan(__a) || __isnan(__b)

  746.             || __isnan(__c) || __isnan(__x))

  747.         return std::numeric_limits<_Tp>::quiet_NaN();

  748.       else if (__c_nint == __c && __c_nint <= _Tp(0))

  749.         return std::numeric_limits<_Tp>::infinity();

  750.       else if (std::abs(__c - __b) < __toler || std::abs(__c - __a) < __toler)

  751.         return std::pow(_Tp(1) - __x, __c - __a - __b);

  752.       else if (__a >= _Tp(0) && __b >= _Tp(0) && __c >= _Tp(0)

  753.             && __x >= _Tp(0) && __x < _Tp(0.995L))

  754.         return __hyperg_series(__a, __b, __c, __x);

  755.       else if (std::abs(__a) < _Tp(10) && std::abs(__b) < _Tp(10))

  756.         {

  757.           //  For integer a and b the hypergeometric function is a

  758.           //  finite polynomial.

  759.           if (__a < _Tp(0)  &&  std::abs(__a - __a_nint) < __toler)

  760.             return __hyperg_series(__a_nint, __b, __c, __x);

  761.           else if (__b < _Tp(0)  &&  std::abs(__b - __b_nint) < __toler)

  762.             return __hyperg_series(__a, __b_nint, __c, __x);

  763.           else if (__x < -_Tp(0.25L))

  764.             return __hyperg_luke(__a, __b, __c, __x);

  765.           else if (__x < _Tp(0.5L))

  766.             return __hyperg_series(__a, __b, __c, __x);

  767.           else

  768.             if (std::abs(__c) > _Tp(10))

  769.               return __hyperg_series(__a, __b, __c, __x);

  770.             else

  771.               return __hyperg_reflect(__a, __b, __c, __x);

  772.         }

  773.       else

  774.         return __hyperg_luke(__a, __b, __c, __x);

  775.     }

  776.   } // namespace __detail

  777. #undef _GLIBCXX_MATH_NS

  778. #if ! _GLIBCXX_USE_STD_SPEC_FUNCS && defined(_GLIBCXX_TR1_CMATH)

  779. } // namespace tr1

  780. #endif

  781. 

  782. _GLIBCXX_END_NAMESPACE_VERSION

  783. }

  784. 

  785. #endif // _GLIBCXX_TR1_HYPERGEOMETRIC_TCC