163d1a8abSmrg // Special functions -*- C++ -*- 263d1a8abSmrg 3*ec02198aSmrg // Copyright (C) 2006-2020 Free Software Foundation, Inc. 463d1a8abSmrg // 563d1a8abSmrg // This file is part of the GNU ISO C++ Library. This library is free 663d1a8abSmrg // software; you can redistribute it and/or modify it under the 763d1a8abSmrg // terms of the GNU General Public License as published by the 863d1a8abSmrg // Free Software Foundation; either version 3, or (at your option) 963d1a8abSmrg // any later version. 1063d1a8abSmrg // 1163d1a8abSmrg // This library is distributed in the hope that it will be useful, 1263d1a8abSmrg // but WITHOUT ANY WARRANTY; without even the implied warranty of 1363d1a8abSmrg // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 1463d1a8abSmrg // GNU General Public License for more details. 1563d1a8abSmrg // 1663d1a8abSmrg // Under Section 7 of GPL version 3, you are granted additional 1763d1a8abSmrg // permissions described in the GCC Runtime Library Exception, version 1863d1a8abSmrg // 3.1, as published by the Free Software Foundation. 1963d1a8abSmrg 2063d1a8abSmrg // You should have received a copy of the GNU General Public License and 2163d1a8abSmrg // a copy of the GCC Runtime Library Exception along with this program; 2263d1a8abSmrg // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 2363d1a8abSmrg // <http://www.gnu.org/licenses/>. 2463d1a8abSmrg 2563d1a8abSmrg /** @file tr1/hypergeometric.tcc 2663d1a8abSmrg * This is an internal header file, included by other library headers. 2763d1a8abSmrg * Do not attempt to use it directly. @headername{tr1/cmath} 2863d1a8abSmrg */ 2963d1a8abSmrg 3063d1a8abSmrg // 3163d1a8abSmrg // ISO C++ 14882 TR1: 5.2 Special functions 3263d1a8abSmrg // 3363d1a8abSmrg 3463d1a8abSmrg // Written by Edward Smith-Rowland based: 3563d1a8abSmrg // (1) Handbook of Mathematical Functions, 3663d1a8abSmrg // ed. Milton Abramowitz and Irene A. Stegun, 3763d1a8abSmrg // Dover Publications, 3863d1a8abSmrg // Section 6, pp. 555-566 3963d1a8abSmrg // (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl 4063d1a8abSmrg 4163d1a8abSmrg #ifndef _GLIBCXX_TR1_HYPERGEOMETRIC_TCC 4263d1a8abSmrg #define _GLIBCXX_TR1_HYPERGEOMETRIC_TCC 1 4363d1a8abSmrg 4463d1a8abSmrg namespace std _GLIBCXX_VISIBILITY(default) 4563d1a8abSmrg { 46c7a68eb7Smrg _GLIBCXX_BEGIN_NAMESPACE_VERSION 47c7a68eb7Smrg 4863d1a8abSmrg #if _GLIBCXX_USE_STD_SPEC_FUNCS 4963d1a8abSmrg # define _GLIBCXX_MATH_NS ::std 5063d1a8abSmrg #elif defined(_GLIBCXX_TR1_CMATH) 5163d1a8abSmrg namespace tr1 5263d1a8abSmrg { 5363d1a8abSmrg # define _GLIBCXX_MATH_NS ::std::tr1 5463d1a8abSmrg #else 5563d1a8abSmrg # error do not include this header directly, use <cmath> or <tr1/cmath> 5663d1a8abSmrg #endif 5763d1a8abSmrg // [5.2] Special functions 5863d1a8abSmrg 5963d1a8abSmrg // Implementation-space details. 6063d1a8abSmrg namespace __detail 6163d1a8abSmrg { 6263d1a8abSmrg /** 6363d1a8abSmrg * @brief This routine returns the confluent hypergeometric function 6463d1a8abSmrg * by series expansion. 6563d1a8abSmrg * 6663d1a8abSmrg * @f[ 6763d1a8abSmrg * _1F_1(a;c;x) = \frac{\Gamma(c)}{\Gamma(a)} 6863d1a8abSmrg * \sum_{n=0}^{\infty} 6963d1a8abSmrg * \frac{\Gamma(a+n)}{\Gamma(c+n)} 7063d1a8abSmrg * \frac{x^n}{n!} 7163d1a8abSmrg * @f] 7263d1a8abSmrg * 7363d1a8abSmrg * If a and b are integers and a < 0 and either b > 0 or b < a 7463d1a8abSmrg * then the series is a polynomial with a finite number of 7563d1a8abSmrg * terms. If b is an integer and b <= 0 the confluent 7663d1a8abSmrg * hypergeometric function is undefined. 7763d1a8abSmrg * 7863d1a8abSmrg * @param __a The "numerator" parameter. 7963d1a8abSmrg * @param __c The "denominator" parameter. 8063d1a8abSmrg * @param __x The argument of the confluent hypergeometric function. 8163d1a8abSmrg * @return The confluent hypergeometric function. 8263d1a8abSmrg */ 8363d1a8abSmrg template<typename _Tp> 8463d1a8abSmrg _Tp __conf_hyperg_series(_Tp __a,_Tp __c,_Tp __x)8563d1a8abSmrg __conf_hyperg_series(_Tp __a, _Tp __c, _Tp __x) 8663d1a8abSmrg { 8763d1a8abSmrg const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 8863d1a8abSmrg 8963d1a8abSmrg _Tp __term = _Tp(1); 9063d1a8abSmrg _Tp __Fac = _Tp(1); 9163d1a8abSmrg const unsigned int __max_iter = 100000; 9263d1a8abSmrg unsigned int __i; 9363d1a8abSmrg for (__i = 0; __i < __max_iter; ++__i) 9463d1a8abSmrg { 9563d1a8abSmrg __term *= (__a + _Tp(__i)) * __x 9663d1a8abSmrg / ((__c + _Tp(__i)) * _Tp(1 + __i)); 9763d1a8abSmrg if (std::abs(__term) < __eps) 9863d1a8abSmrg { 9963d1a8abSmrg break; 10063d1a8abSmrg } 10163d1a8abSmrg __Fac += __term; 10263d1a8abSmrg } 10363d1a8abSmrg if (__i == __max_iter) 10463d1a8abSmrg std::__throw_runtime_error(__N("Series failed to converge " 10563d1a8abSmrg "in __conf_hyperg_series.")); 10663d1a8abSmrg 10763d1a8abSmrg return __Fac; 10863d1a8abSmrg } 10963d1a8abSmrg 11063d1a8abSmrg 11163d1a8abSmrg /** 11263d1a8abSmrg * @brief Return the hypogeometric function @f$ _2F_1(a,b;c;x) @f$ 11363d1a8abSmrg * by an iterative procedure described in 11463d1a8abSmrg * Luke, Algorithms for the Computation of Mathematical Functions. 11563d1a8abSmrg * 11663d1a8abSmrg * Like the case of the 2F1 rational approximations, these are 11763d1a8abSmrg * probably guaranteed to converge for x < 0, barring gross 11863d1a8abSmrg * numerical instability in the pre-asymptotic regime. 11963d1a8abSmrg */ 12063d1a8abSmrg template<typename _Tp> 12163d1a8abSmrg _Tp __conf_hyperg_luke(_Tp __a,_Tp __c,_Tp __xin)12263d1a8abSmrg __conf_hyperg_luke(_Tp __a, _Tp __c, _Tp __xin) 12363d1a8abSmrg { 12463d1a8abSmrg const _Tp __big = std::pow(std::numeric_limits<_Tp>::max(), _Tp(0.16L)); 12563d1a8abSmrg const int __nmax = 20000; 12663d1a8abSmrg const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 12763d1a8abSmrg const _Tp __x = -__xin; 12863d1a8abSmrg const _Tp __x3 = __x * __x * __x; 12963d1a8abSmrg const _Tp __t0 = __a / __c; 13063d1a8abSmrg const _Tp __t1 = (__a + _Tp(1)) / (_Tp(2) * __c); 13163d1a8abSmrg const _Tp __t2 = (__a + _Tp(2)) / (_Tp(2) * (__c + _Tp(1))); 13263d1a8abSmrg _Tp __F = _Tp(1); 13363d1a8abSmrg _Tp __prec; 13463d1a8abSmrg 13563d1a8abSmrg _Tp __Bnm3 = _Tp(1); 13663d1a8abSmrg _Tp __Bnm2 = _Tp(1) + __t1 * __x; 13763d1a8abSmrg _Tp __Bnm1 = _Tp(1) + __t2 * __x * (_Tp(1) + __t1 / _Tp(3) * __x); 13863d1a8abSmrg 13963d1a8abSmrg _Tp __Anm3 = _Tp(1); 14063d1a8abSmrg _Tp __Anm2 = __Bnm2 - __t0 * __x; 14163d1a8abSmrg _Tp __Anm1 = __Bnm1 - __t0 * (_Tp(1) + __t2 * __x) * __x 14263d1a8abSmrg + __t0 * __t1 * (__c / (__c + _Tp(1))) * __x * __x; 14363d1a8abSmrg 14463d1a8abSmrg int __n = 3; 14563d1a8abSmrg while(1) 14663d1a8abSmrg { 14763d1a8abSmrg _Tp __npam1 = _Tp(__n - 1) + __a; 14863d1a8abSmrg _Tp __npcm1 = _Tp(__n - 1) + __c; 14963d1a8abSmrg _Tp __npam2 = _Tp(__n - 2) + __a; 15063d1a8abSmrg _Tp __npcm2 = _Tp(__n - 2) + __c; 15163d1a8abSmrg _Tp __tnm1 = _Tp(2 * __n - 1); 15263d1a8abSmrg _Tp __tnm3 = _Tp(2 * __n - 3); 15363d1a8abSmrg _Tp __tnm5 = _Tp(2 * __n - 5); 15463d1a8abSmrg _Tp __F1 = (_Tp(__n - 2) - __a) / (_Tp(2) * __tnm3 * __npcm1); 15563d1a8abSmrg _Tp __F2 = (_Tp(__n) + __a) * __npam1 15663d1a8abSmrg / (_Tp(4) * __tnm1 * __tnm3 * __npcm2 * __npcm1); 15763d1a8abSmrg _Tp __F3 = -__npam2 * __npam1 * (_Tp(__n - 2) - __a) 15863d1a8abSmrg / (_Tp(8) * __tnm3 * __tnm3 * __tnm5 15963d1a8abSmrg * (_Tp(__n - 3) + __c) * __npcm2 * __npcm1); 16063d1a8abSmrg _Tp __E = -__npam1 * (_Tp(__n - 1) - __c) 16163d1a8abSmrg / (_Tp(2) * __tnm3 * __npcm2 * __npcm1); 16263d1a8abSmrg 16363d1a8abSmrg _Tp __An = (_Tp(1) + __F1 * __x) * __Anm1 16463d1a8abSmrg + (__E + __F2 * __x) * __x * __Anm2 + __F3 * __x3 * __Anm3; 16563d1a8abSmrg _Tp __Bn = (_Tp(1) + __F1 * __x) * __Bnm1 16663d1a8abSmrg + (__E + __F2 * __x) * __x * __Bnm2 + __F3 * __x3 * __Bnm3; 16763d1a8abSmrg _Tp __r = __An / __Bn; 16863d1a8abSmrg 16963d1a8abSmrg __prec = std::abs((__F - __r) / __F); 17063d1a8abSmrg __F = __r; 17163d1a8abSmrg 17263d1a8abSmrg if (__prec < __eps || __n > __nmax) 17363d1a8abSmrg break; 17463d1a8abSmrg 17563d1a8abSmrg if (std::abs(__An) > __big || std::abs(__Bn) > __big) 17663d1a8abSmrg { 17763d1a8abSmrg __An /= __big; 17863d1a8abSmrg __Bn /= __big; 17963d1a8abSmrg __Anm1 /= __big; 18063d1a8abSmrg __Bnm1 /= __big; 18163d1a8abSmrg __Anm2 /= __big; 18263d1a8abSmrg __Bnm2 /= __big; 18363d1a8abSmrg __Anm3 /= __big; 18463d1a8abSmrg __Bnm3 /= __big; 18563d1a8abSmrg } 18663d1a8abSmrg else if (std::abs(__An) < _Tp(1) / __big 18763d1a8abSmrg || std::abs(__Bn) < _Tp(1) / __big) 18863d1a8abSmrg { 18963d1a8abSmrg __An *= __big; 19063d1a8abSmrg __Bn *= __big; 19163d1a8abSmrg __Anm1 *= __big; 19263d1a8abSmrg __Bnm1 *= __big; 19363d1a8abSmrg __Anm2 *= __big; 19463d1a8abSmrg __Bnm2 *= __big; 19563d1a8abSmrg __Anm3 *= __big; 19663d1a8abSmrg __Bnm3 *= __big; 19763d1a8abSmrg } 19863d1a8abSmrg 19963d1a8abSmrg ++__n; 20063d1a8abSmrg __Bnm3 = __Bnm2; 20163d1a8abSmrg __Bnm2 = __Bnm1; 20263d1a8abSmrg __Bnm1 = __Bn; 20363d1a8abSmrg __Anm3 = __Anm2; 20463d1a8abSmrg __Anm2 = __Anm1; 20563d1a8abSmrg __Anm1 = __An; 20663d1a8abSmrg } 20763d1a8abSmrg 20863d1a8abSmrg if (__n >= __nmax) 20963d1a8abSmrg std::__throw_runtime_error(__N("Iteration failed to converge " 21063d1a8abSmrg "in __conf_hyperg_luke.")); 21163d1a8abSmrg 21263d1a8abSmrg return __F; 21363d1a8abSmrg } 21463d1a8abSmrg 21563d1a8abSmrg 21663d1a8abSmrg /** 21763d1a8abSmrg * @brief Return the confluent hypogeometric function 21863d1a8abSmrg * @f$ _1F_1(a;c;x) @f$. 21963d1a8abSmrg * 22063d1a8abSmrg * @todo Handle b == nonpositive integer blowup - return NaN. 22163d1a8abSmrg * 22263d1a8abSmrg * @param __a The @a numerator parameter. 22363d1a8abSmrg * @param __c The @a denominator parameter. 22463d1a8abSmrg * @param __x The argument of the confluent hypergeometric function. 22563d1a8abSmrg * @return The confluent hypergeometric function. 22663d1a8abSmrg */ 22763d1a8abSmrg template<typename _Tp> 22863d1a8abSmrg _Tp __conf_hyperg(_Tp __a,_Tp __c,_Tp __x)22963d1a8abSmrg __conf_hyperg(_Tp __a, _Tp __c, _Tp __x) 23063d1a8abSmrg { 23163d1a8abSmrg #if _GLIBCXX_USE_C99_MATH_TR1 23263d1a8abSmrg const _Tp __c_nint = _GLIBCXX_MATH_NS::nearbyint(__c); 23363d1a8abSmrg #else 23463d1a8abSmrg const _Tp __c_nint = static_cast<int>(__c + _Tp(0.5L)); 23563d1a8abSmrg #endif 23663d1a8abSmrg if (__isnan(__a) || __isnan(__c) || __isnan(__x)) 23763d1a8abSmrg return std::numeric_limits<_Tp>::quiet_NaN(); 23863d1a8abSmrg else if (__c_nint == __c && __c_nint <= 0) 23963d1a8abSmrg return std::numeric_limits<_Tp>::infinity(); 24063d1a8abSmrg else if (__a == _Tp(0)) 24163d1a8abSmrg return _Tp(1); 24263d1a8abSmrg else if (__c == __a) 24363d1a8abSmrg return std::exp(__x); 24463d1a8abSmrg else if (__x < _Tp(0)) 24563d1a8abSmrg return __conf_hyperg_luke(__a, __c, __x); 24663d1a8abSmrg else 24763d1a8abSmrg return __conf_hyperg_series(__a, __c, __x); 24863d1a8abSmrg } 24963d1a8abSmrg 25063d1a8abSmrg 25163d1a8abSmrg /** 25263d1a8abSmrg * @brief Return the hypogeometric function @f$ _2F_1(a,b;c;x) @f$ 25363d1a8abSmrg * by series expansion. 25463d1a8abSmrg * 25563d1a8abSmrg * The hypogeometric function is defined by 25663d1a8abSmrg * @f[ 25763d1a8abSmrg * _2F_1(a,b;c;x) = \frac{\Gamma(c)}{\Gamma(a)\Gamma(b)} 25863d1a8abSmrg * \sum_{n=0}^{\infty} 25963d1a8abSmrg * \frac{\Gamma(a+n)\Gamma(b+n)}{\Gamma(c+n)} 26063d1a8abSmrg * \frac{x^n}{n!} 26163d1a8abSmrg * @f] 26263d1a8abSmrg * 26363d1a8abSmrg * This works and it's pretty fast. 26463d1a8abSmrg * 26563d1a8abSmrg * @param __a The first @a numerator parameter. 26663d1a8abSmrg * @param __a The second @a numerator parameter. 26763d1a8abSmrg * @param __c The @a denominator parameter. 26863d1a8abSmrg * @param __x The argument of the confluent hypergeometric function. 26963d1a8abSmrg * @return The confluent hypergeometric function. 27063d1a8abSmrg */ 27163d1a8abSmrg template<typename _Tp> 27263d1a8abSmrg _Tp __hyperg_series(_Tp __a,_Tp __b,_Tp __c,_Tp __x)27363d1a8abSmrg __hyperg_series(_Tp __a, _Tp __b, _Tp __c, _Tp __x) 27463d1a8abSmrg { 27563d1a8abSmrg const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 27663d1a8abSmrg 27763d1a8abSmrg _Tp __term = _Tp(1); 27863d1a8abSmrg _Tp __Fabc = _Tp(1); 27963d1a8abSmrg const unsigned int __max_iter = 100000; 28063d1a8abSmrg unsigned int __i; 28163d1a8abSmrg for (__i = 0; __i < __max_iter; ++__i) 28263d1a8abSmrg { 28363d1a8abSmrg __term *= (__a + _Tp(__i)) * (__b + _Tp(__i)) * __x 28463d1a8abSmrg / ((__c + _Tp(__i)) * _Tp(1 + __i)); 28563d1a8abSmrg if (std::abs(__term) < __eps) 28663d1a8abSmrg { 28763d1a8abSmrg break; 28863d1a8abSmrg } 28963d1a8abSmrg __Fabc += __term; 29063d1a8abSmrg } 29163d1a8abSmrg if (__i == __max_iter) 29263d1a8abSmrg std::__throw_runtime_error(__N("Series failed to converge " 29363d1a8abSmrg "in __hyperg_series.")); 29463d1a8abSmrg 29563d1a8abSmrg return __Fabc; 29663d1a8abSmrg } 29763d1a8abSmrg 29863d1a8abSmrg 29963d1a8abSmrg /** 30063d1a8abSmrg * @brief Return the hypogeometric function @f$ _2F_1(a,b;c;x) @f$ 30163d1a8abSmrg * by an iterative procedure described in 30263d1a8abSmrg * Luke, Algorithms for the Computation of Mathematical Functions. 30363d1a8abSmrg */ 30463d1a8abSmrg template<typename _Tp> 30563d1a8abSmrg _Tp __hyperg_luke(_Tp __a,_Tp __b,_Tp __c,_Tp __xin)30663d1a8abSmrg __hyperg_luke(_Tp __a, _Tp __b, _Tp __c, _Tp __xin) 30763d1a8abSmrg { 30863d1a8abSmrg const _Tp __big = std::pow(std::numeric_limits<_Tp>::max(), _Tp(0.16L)); 30963d1a8abSmrg const int __nmax = 20000; 31063d1a8abSmrg const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 31163d1a8abSmrg const _Tp __x = -__xin; 31263d1a8abSmrg const _Tp __x3 = __x * __x * __x; 31363d1a8abSmrg const _Tp __t0 = __a * __b / __c; 31463d1a8abSmrg const _Tp __t1 = (__a + _Tp(1)) * (__b + _Tp(1)) / (_Tp(2) * __c); 31563d1a8abSmrg const _Tp __t2 = (__a + _Tp(2)) * (__b + _Tp(2)) 31663d1a8abSmrg / (_Tp(2) * (__c + _Tp(1))); 31763d1a8abSmrg 31863d1a8abSmrg _Tp __F = _Tp(1); 31963d1a8abSmrg 32063d1a8abSmrg _Tp __Bnm3 = _Tp(1); 32163d1a8abSmrg _Tp __Bnm2 = _Tp(1) + __t1 * __x; 32263d1a8abSmrg _Tp __Bnm1 = _Tp(1) + __t2 * __x * (_Tp(1) + __t1 / _Tp(3) * __x); 32363d1a8abSmrg 32463d1a8abSmrg _Tp __Anm3 = _Tp(1); 32563d1a8abSmrg _Tp __Anm2 = __Bnm2 - __t0 * __x; 32663d1a8abSmrg _Tp __Anm1 = __Bnm1 - __t0 * (_Tp(1) + __t2 * __x) * __x 32763d1a8abSmrg + __t0 * __t1 * (__c / (__c + _Tp(1))) * __x * __x; 32863d1a8abSmrg 32963d1a8abSmrg int __n = 3; 33063d1a8abSmrg while (1) 33163d1a8abSmrg { 33263d1a8abSmrg const _Tp __npam1 = _Tp(__n - 1) + __a; 33363d1a8abSmrg const _Tp __npbm1 = _Tp(__n - 1) + __b; 33463d1a8abSmrg const _Tp __npcm1 = _Tp(__n - 1) + __c; 33563d1a8abSmrg const _Tp __npam2 = _Tp(__n - 2) + __a; 33663d1a8abSmrg const _Tp __npbm2 = _Tp(__n - 2) + __b; 33763d1a8abSmrg const _Tp __npcm2 = _Tp(__n - 2) + __c; 33863d1a8abSmrg const _Tp __tnm1 = _Tp(2 * __n - 1); 33963d1a8abSmrg const _Tp __tnm3 = _Tp(2 * __n - 3); 34063d1a8abSmrg const _Tp __tnm5 = _Tp(2 * __n - 5); 34163d1a8abSmrg const _Tp __n2 = __n * __n; 34263d1a8abSmrg const _Tp __F1 = (_Tp(3) * __n2 + (__a + __b - _Tp(6)) * __n 34363d1a8abSmrg + _Tp(2) - __a * __b - _Tp(2) * (__a + __b)) 34463d1a8abSmrg / (_Tp(2) * __tnm3 * __npcm1); 34563d1a8abSmrg const _Tp __F2 = -(_Tp(3) * __n2 - (__a + __b + _Tp(6)) * __n 34663d1a8abSmrg + _Tp(2) - __a * __b) * __npam1 * __npbm1 34763d1a8abSmrg / (_Tp(4) * __tnm1 * __tnm3 * __npcm2 * __npcm1); 34863d1a8abSmrg const _Tp __F3 = (__npam2 * __npam1 * __npbm2 * __npbm1 34963d1a8abSmrg * (_Tp(__n - 2) - __a) * (_Tp(__n - 2) - __b)) 35063d1a8abSmrg / (_Tp(8) * __tnm3 * __tnm3 * __tnm5 35163d1a8abSmrg * (_Tp(__n - 3) + __c) * __npcm2 * __npcm1); 35263d1a8abSmrg const _Tp __E = -__npam1 * __npbm1 * (_Tp(__n - 1) - __c) 35363d1a8abSmrg / (_Tp(2) * __tnm3 * __npcm2 * __npcm1); 35463d1a8abSmrg 35563d1a8abSmrg _Tp __An = (_Tp(1) + __F1 * __x) * __Anm1 35663d1a8abSmrg + (__E + __F2 * __x) * __x * __Anm2 + __F3 * __x3 * __Anm3; 35763d1a8abSmrg _Tp __Bn = (_Tp(1) + __F1 * __x) * __Bnm1 35863d1a8abSmrg + (__E + __F2 * __x) * __x * __Bnm2 + __F3 * __x3 * __Bnm3; 35963d1a8abSmrg const _Tp __r = __An / __Bn; 36063d1a8abSmrg 36163d1a8abSmrg const _Tp __prec = std::abs((__F - __r) / __F); 36263d1a8abSmrg __F = __r; 36363d1a8abSmrg 36463d1a8abSmrg if (__prec < __eps || __n > __nmax) 36563d1a8abSmrg break; 36663d1a8abSmrg 36763d1a8abSmrg if (std::abs(__An) > __big || std::abs(__Bn) > __big) 36863d1a8abSmrg { 36963d1a8abSmrg __An /= __big; 37063d1a8abSmrg __Bn /= __big; 37163d1a8abSmrg __Anm1 /= __big; 37263d1a8abSmrg __Bnm1 /= __big; 37363d1a8abSmrg __Anm2 /= __big; 37463d1a8abSmrg __Bnm2 /= __big; 37563d1a8abSmrg __Anm3 /= __big; 37663d1a8abSmrg __Bnm3 /= __big; 37763d1a8abSmrg } 37863d1a8abSmrg else if (std::abs(__An) < _Tp(1) / __big 37963d1a8abSmrg || std::abs(__Bn) < _Tp(1) / __big) 38063d1a8abSmrg { 38163d1a8abSmrg __An *= __big; 38263d1a8abSmrg __Bn *= __big; 38363d1a8abSmrg __Anm1 *= __big; 38463d1a8abSmrg __Bnm1 *= __big; 38563d1a8abSmrg __Anm2 *= __big; 38663d1a8abSmrg __Bnm2 *= __big; 38763d1a8abSmrg __Anm3 *= __big; 38863d1a8abSmrg __Bnm3 *= __big; 38963d1a8abSmrg } 39063d1a8abSmrg 39163d1a8abSmrg ++__n; 39263d1a8abSmrg __Bnm3 = __Bnm2; 39363d1a8abSmrg __Bnm2 = __Bnm1; 39463d1a8abSmrg __Bnm1 = __Bn; 39563d1a8abSmrg __Anm3 = __Anm2; 39663d1a8abSmrg __Anm2 = __Anm1; 39763d1a8abSmrg __Anm1 = __An; 39863d1a8abSmrg } 39963d1a8abSmrg 40063d1a8abSmrg if (__n >= __nmax) 40163d1a8abSmrg std::__throw_runtime_error(__N("Iteration failed to converge " 40263d1a8abSmrg "in __hyperg_luke.")); 40363d1a8abSmrg 40463d1a8abSmrg return __F; 40563d1a8abSmrg } 40663d1a8abSmrg 40763d1a8abSmrg 40863d1a8abSmrg /** 40963d1a8abSmrg * @brief Return the hypogeometric function @f$ _2F_1(a,b;c;x) @f$ 41063d1a8abSmrg * by the reflection formulae in Abramowitz & Stegun formula 41163d1a8abSmrg * 15.3.6 for d = c - a - b not integral and formula 15.3.11 for 41263d1a8abSmrg * d = c - a - b integral. This assumes a, b, c != negative 41363d1a8abSmrg * integer. 41463d1a8abSmrg * 41563d1a8abSmrg * The hypogeometric function is defined by 41663d1a8abSmrg * @f[ 41763d1a8abSmrg * _2F_1(a,b;c;x) = \frac{\Gamma(c)}{\Gamma(a)\Gamma(b)} 41863d1a8abSmrg * \sum_{n=0}^{\infty} 41963d1a8abSmrg * \frac{\Gamma(a+n)\Gamma(b+n)}{\Gamma(c+n)} 42063d1a8abSmrg * \frac{x^n}{n!} 42163d1a8abSmrg * @f] 42263d1a8abSmrg * 42363d1a8abSmrg * The reflection formula for nonintegral @f$ d = c - a - b @f$ is: 42463d1a8abSmrg * @f[ 42563d1a8abSmrg * _2F_1(a,b;c;x) = \frac{\Gamma(c)\Gamma(d)}{\Gamma(c-a)\Gamma(c-b)} 42663d1a8abSmrg * _2F_1(a,b;1-d;1-x) 42763d1a8abSmrg * + \frac{\Gamma(c)\Gamma(-d)}{\Gamma(a)\Gamma(b)} 42863d1a8abSmrg * _2F_1(c-a,c-b;1+d;1-x) 42963d1a8abSmrg * @f] 43063d1a8abSmrg * 43163d1a8abSmrg * The reflection formula for integral @f$ m = c - a - b @f$ is: 43263d1a8abSmrg * @f[ 43363d1a8abSmrg * _2F_1(a,b;a+b+m;x) = \frac{\Gamma(m)\Gamma(a+b+m)}{\Gamma(a+m)\Gamma(b+m)} 43463d1a8abSmrg * \sum_{k=0}^{m-1} \frac{(m+a)_k(m+b)_k}{k!(1-m)_k} 43563d1a8abSmrg * - 43663d1a8abSmrg * @f] 43763d1a8abSmrg */ 43863d1a8abSmrg template<typename _Tp> 43963d1a8abSmrg _Tp __hyperg_reflect(_Tp __a,_Tp __b,_Tp __c,_Tp __x)44063d1a8abSmrg __hyperg_reflect(_Tp __a, _Tp __b, _Tp __c, _Tp __x) 44163d1a8abSmrg { 44263d1a8abSmrg const _Tp __d = __c - __a - __b; 44363d1a8abSmrg const int __intd = std::floor(__d + _Tp(0.5L)); 44463d1a8abSmrg const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); 44563d1a8abSmrg const _Tp __toler = _Tp(1000) * __eps; 44663d1a8abSmrg const _Tp __log_max = std::log(std::numeric_limits<_Tp>::max()); 44763d1a8abSmrg const bool __d_integer = (std::abs(__d - __intd) < __toler); 44863d1a8abSmrg 44963d1a8abSmrg if (__d_integer) 45063d1a8abSmrg { 45163d1a8abSmrg const _Tp __ln_omx = std::log(_Tp(1) - __x); 45263d1a8abSmrg const _Tp __ad = std::abs(__d); 45363d1a8abSmrg _Tp __F1, __F2; 45463d1a8abSmrg 45563d1a8abSmrg _Tp __d1, __d2; 45663d1a8abSmrg if (__d >= _Tp(0)) 45763d1a8abSmrg { 45863d1a8abSmrg __d1 = __d; 45963d1a8abSmrg __d2 = _Tp(0); 46063d1a8abSmrg } 46163d1a8abSmrg else 46263d1a8abSmrg { 46363d1a8abSmrg __d1 = _Tp(0); 46463d1a8abSmrg __d2 = __d; 46563d1a8abSmrg } 46663d1a8abSmrg 46763d1a8abSmrg const _Tp __lng_c = __log_gamma(__c); 46863d1a8abSmrg 46963d1a8abSmrg // Evaluate F1. 47063d1a8abSmrg if (__ad < __eps) 47163d1a8abSmrg { 47263d1a8abSmrg // d = c - a - b = 0. 47363d1a8abSmrg __F1 = _Tp(0); 47463d1a8abSmrg } 47563d1a8abSmrg else 47663d1a8abSmrg { 47763d1a8abSmrg 47863d1a8abSmrg bool __ok_d1 = true; 47963d1a8abSmrg _Tp __lng_ad, __lng_ad1, __lng_bd1; 48063d1a8abSmrg __try 48163d1a8abSmrg { 48263d1a8abSmrg __lng_ad = __log_gamma(__ad); 48363d1a8abSmrg __lng_ad1 = __log_gamma(__a + __d1); 48463d1a8abSmrg __lng_bd1 = __log_gamma(__b + __d1); 48563d1a8abSmrg } 48663d1a8abSmrg __catch(...) 48763d1a8abSmrg { 48863d1a8abSmrg __ok_d1 = false; 48963d1a8abSmrg } 49063d1a8abSmrg 49163d1a8abSmrg if (__ok_d1) 49263d1a8abSmrg { 49363d1a8abSmrg /* Gamma functions in the denominator are ok. 49463d1a8abSmrg * Proceed with evaluation. 49563d1a8abSmrg */ 49663d1a8abSmrg _Tp __sum1 = _Tp(1); 49763d1a8abSmrg _Tp __term = _Tp(1); 49863d1a8abSmrg _Tp __ln_pre1 = __lng_ad + __lng_c + __d2 * __ln_omx 49963d1a8abSmrg - __lng_ad1 - __lng_bd1; 50063d1a8abSmrg 50163d1a8abSmrg /* Do F1 sum. 50263d1a8abSmrg */ 50363d1a8abSmrg for (int __i = 1; __i < __ad; ++__i) 50463d1a8abSmrg { 50563d1a8abSmrg const int __j = __i - 1; 50663d1a8abSmrg __term *= (__a + __d2 + __j) * (__b + __d2 + __j) 50763d1a8abSmrg / (_Tp(1) + __d2 + __j) / __i * (_Tp(1) - __x); 50863d1a8abSmrg __sum1 += __term; 50963d1a8abSmrg } 51063d1a8abSmrg 51163d1a8abSmrg if (__ln_pre1 > __log_max) 51263d1a8abSmrg std::__throw_runtime_error(__N("Overflow of gamma functions" 51363d1a8abSmrg " in __hyperg_luke.")); 51463d1a8abSmrg else 51563d1a8abSmrg __F1 = std::exp(__ln_pre1) * __sum1; 51663d1a8abSmrg } 51763d1a8abSmrg else 51863d1a8abSmrg { 51963d1a8abSmrg // Gamma functions in the denominator were not ok. 52063d1a8abSmrg // So the F1 term is zero. 52163d1a8abSmrg __F1 = _Tp(0); 52263d1a8abSmrg } 52363d1a8abSmrg } // end F1 evaluation 52463d1a8abSmrg 52563d1a8abSmrg // Evaluate F2. 52663d1a8abSmrg bool __ok_d2 = true; 52763d1a8abSmrg _Tp __lng_ad2, __lng_bd2; 52863d1a8abSmrg __try 52963d1a8abSmrg { 53063d1a8abSmrg __lng_ad2 = __log_gamma(__a + __d2); 53163d1a8abSmrg __lng_bd2 = __log_gamma(__b + __d2); 53263d1a8abSmrg } 53363d1a8abSmrg __catch(...) 53463d1a8abSmrg { 53563d1a8abSmrg __ok_d2 = false; 53663d1a8abSmrg } 53763d1a8abSmrg 53863d1a8abSmrg if (__ok_d2) 53963d1a8abSmrg { 54063d1a8abSmrg // Gamma functions in the denominator are ok. 54163d1a8abSmrg // Proceed with evaluation. 54263d1a8abSmrg const int __maxiter = 2000; 54363d1a8abSmrg const _Tp __psi_1 = -__numeric_constants<_Tp>::__gamma_e(); 54463d1a8abSmrg const _Tp __psi_1pd = __psi(_Tp(1) + __ad); 54563d1a8abSmrg const _Tp __psi_apd1 = __psi(__a + __d1); 54663d1a8abSmrg const _Tp __psi_bpd1 = __psi(__b + __d1); 54763d1a8abSmrg 54863d1a8abSmrg _Tp __psi_term = __psi_1 + __psi_1pd - __psi_apd1 54963d1a8abSmrg - __psi_bpd1 - __ln_omx; 55063d1a8abSmrg _Tp __fact = _Tp(1); 55163d1a8abSmrg _Tp __sum2 = __psi_term; 55263d1a8abSmrg _Tp __ln_pre2 = __lng_c + __d1 * __ln_omx 55363d1a8abSmrg - __lng_ad2 - __lng_bd2; 55463d1a8abSmrg 55563d1a8abSmrg // Do F2 sum. 55663d1a8abSmrg int __j; 55763d1a8abSmrg for (__j = 1; __j < __maxiter; ++__j) 55863d1a8abSmrg { 55963d1a8abSmrg // Values for psi functions use recurrence; 56063d1a8abSmrg // Abramowitz & Stegun 6.3.5 56163d1a8abSmrg const _Tp __term1 = _Tp(1) / _Tp(__j) 56263d1a8abSmrg + _Tp(1) / (__ad + __j); 56363d1a8abSmrg const _Tp __term2 = _Tp(1) / (__a + __d1 + _Tp(__j - 1)) 56463d1a8abSmrg + _Tp(1) / (__b + __d1 + _Tp(__j - 1)); 56563d1a8abSmrg __psi_term += __term1 - __term2; 56663d1a8abSmrg __fact *= (__a + __d1 + _Tp(__j - 1)) 56763d1a8abSmrg * (__b + __d1 + _Tp(__j - 1)) 56863d1a8abSmrg / ((__ad + __j) * __j) * (_Tp(1) - __x); 56963d1a8abSmrg const _Tp __delta = __fact * __psi_term; 57063d1a8abSmrg __sum2 += __delta; 57163d1a8abSmrg if (std::abs(__delta) < __eps * std::abs(__sum2)) 57263d1a8abSmrg break; 57363d1a8abSmrg } 57463d1a8abSmrg if (__j == __maxiter) 57563d1a8abSmrg std::__throw_runtime_error(__N("Sum F2 failed to converge " 57663d1a8abSmrg "in __hyperg_reflect")); 57763d1a8abSmrg 57863d1a8abSmrg if (__sum2 == _Tp(0)) 57963d1a8abSmrg __F2 = _Tp(0); 58063d1a8abSmrg else 58163d1a8abSmrg __F2 = std::exp(__ln_pre2) * __sum2; 58263d1a8abSmrg } 58363d1a8abSmrg else 58463d1a8abSmrg { 58563d1a8abSmrg // Gamma functions in the denominator not ok. 58663d1a8abSmrg // So the F2 term is zero. 58763d1a8abSmrg __F2 = _Tp(0); 58863d1a8abSmrg } // end F2 evaluation 58963d1a8abSmrg 59063d1a8abSmrg const _Tp __sgn_2 = (__intd % 2 == 1 ? -_Tp(1) : _Tp(1)); 59163d1a8abSmrg const _Tp __F = __F1 + __sgn_2 * __F2; 59263d1a8abSmrg 59363d1a8abSmrg return __F; 59463d1a8abSmrg } 59563d1a8abSmrg else 59663d1a8abSmrg { 59763d1a8abSmrg // d = c - a - b not an integer. 59863d1a8abSmrg 59963d1a8abSmrg // These gamma functions appear in the denominator, so we 60063d1a8abSmrg // catch their harmless domain errors and set the terms to zero. 60163d1a8abSmrg bool __ok1 = true; 60263d1a8abSmrg _Tp __sgn_g1ca = _Tp(0), __ln_g1ca = _Tp(0); 60363d1a8abSmrg _Tp __sgn_g1cb = _Tp(0), __ln_g1cb = _Tp(0); 60463d1a8abSmrg __try 60563d1a8abSmrg { 60663d1a8abSmrg __sgn_g1ca = __log_gamma_sign(__c - __a); 60763d1a8abSmrg __ln_g1ca = __log_gamma(__c - __a); 60863d1a8abSmrg __sgn_g1cb = __log_gamma_sign(__c - __b); 60963d1a8abSmrg __ln_g1cb = __log_gamma(__c - __b); 61063d1a8abSmrg } 61163d1a8abSmrg __catch(...) 61263d1a8abSmrg { 61363d1a8abSmrg __ok1 = false; 61463d1a8abSmrg } 61563d1a8abSmrg 61663d1a8abSmrg bool __ok2 = true; 61763d1a8abSmrg _Tp __sgn_g2a = _Tp(0), __ln_g2a = _Tp(0); 61863d1a8abSmrg _Tp __sgn_g2b = _Tp(0), __ln_g2b = _Tp(0); 61963d1a8abSmrg __try 62063d1a8abSmrg { 62163d1a8abSmrg __sgn_g2a = __log_gamma_sign(__a); 62263d1a8abSmrg __ln_g2a = __log_gamma(__a); 62363d1a8abSmrg __sgn_g2b = __log_gamma_sign(__b); 62463d1a8abSmrg __ln_g2b = __log_gamma(__b); 62563d1a8abSmrg } 62663d1a8abSmrg __catch(...) 62763d1a8abSmrg { 62863d1a8abSmrg __ok2 = false; 62963d1a8abSmrg } 63063d1a8abSmrg 63163d1a8abSmrg const _Tp __sgn_gc = __log_gamma_sign(__c); 63263d1a8abSmrg const _Tp __ln_gc = __log_gamma(__c); 63363d1a8abSmrg const _Tp __sgn_gd = __log_gamma_sign(__d); 63463d1a8abSmrg const _Tp __ln_gd = __log_gamma(__d); 63563d1a8abSmrg const _Tp __sgn_gmd = __log_gamma_sign(-__d); 63663d1a8abSmrg const _Tp __ln_gmd = __log_gamma(-__d); 63763d1a8abSmrg 63863d1a8abSmrg const _Tp __sgn1 = __sgn_gc * __sgn_gd * __sgn_g1ca * __sgn_g1cb; 63963d1a8abSmrg const _Tp __sgn2 = __sgn_gc * __sgn_gmd * __sgn_g2a * __sgn_g2b; 64063d1a8abSmrg 64163d1a8abSmrg _Tp __pre1, __pre2; 64263d1a8abSmrg if (__ok1 && __ok2) 64363d1a8abSmrg { 64463d1a8abSmrg _Tp __ln_pre1 = __ln_gc + __ln_gd - __ln_g1ca - __ln_g1cb; 64563d1a8abSmrg _Tp __ln_pre2 = __ln_gc + __ln_gmd - __ln_g2a - __ln_g2b 64663d1a8abSmrg + __d * std::log(_Tp(1) - __x); 64763d1a8abSmrg if (__ln_pre1 < __log_max && __ln_pre2 < __log_max) 64863d1a8abSmrg { 64963d1a8abSmrg __pre1 = std::exp(__ln_pre1); 65063d1a8abSmrg __pre2 = std::exp(__ln_pre2); 65163d1a8abSmrg __pre1 *= __sgn1; 65263d1a8abSmrg __pre2 *= __sgn2; 65363d1a8abSmrg } 65463d1a8abSmrg else 65563d1a8abSmrg { 65663d1a8abSmrg std::__throw_runtime_error(__N("Overflow of gamma functions " 65763d1a8abSmrg "in __hyperg_reflect")); 65863d1a8abSmrg } 65963d1a8abSmrg } 66063d1a8abSmrg else if (__ok1 && !__ok2) 66163d1a8abSmrg { 66263d1a8abSmrg _Tp __ln_pre1 = __ln_gc + __ln_gd - __ln_g1ca - __ln_g1cb; 66363d1a8abSmrg if (__ln_pre1 < __log_max) 66463d1a8abSmrg { 66563d1a8abSmrg __pre1 = std::exp(__ln_pre1); 66663d1a8abSmrg __pre1 *= __sgn1; 66763d1a8abSmrg __pre2 = _Tp(0); 66863d1a8abSmrg } 66963d1a8abSmrg else 67063d1a8abSmrg { 67163d1a8abSmrg std::__throw_runtime_error(__N("Overflow of gamma functions " 67263d1a8abSmrg "in __hyperg_reflect")); 67363d1a8abSmrg } 67463d1a8abSmrg } 67563d1a8abSmrg else if (!__ok1 && __ok2) 67663d1a8abSmrg { 67763d1a8abSmrg _Tp __ln_pre2 = __ln_gc + __ln_gmd - __ln_g2a - __ln_g2b 67863d1a8abSmrg + __d * std::log(_Tp(1) - __x); 67963d1a8abSmrg if (__ln_pre2 < __log_max) 68063d1a8abSmrg { 68163d1a8abSmrg __pre1 = _Tp(0); 68263d1a8abSmrg __pre2 = std::exp(__ln_pre2); 68363d1a8abSmrg __pre2 *= __sgn2; 68463d1a8abSmrg } 68563d1a8abSmrg else 68663d1a8abSmrg { 68763d1a8abSmrg std::__throw_runtime_error(__N("Overflow of gamma functions " 68863d1a8abSmrg "in __hyperg_reflect")); 68963d1a8abSmrg } 69063d1a8abSmrg } 69163d1a8abSmrg else 69263d1a8abSmrg { 69363d1a8abSmrg __pre1 = _Tp(0); 69463d1a8abSmrg __pre2 = _Tp(0); 69563d1a8abSmrg std::__throw_runtime_error(__N("Underflow of gamma functions " 69663d1a8abSmrg "in __hyperg_reflect")); 69763d1a8abSmrg } 69863d1a8abSmrg 69963d1a8abSmrg const _Tp __F1 = __hyperg_series(__a, __b, _Tp(1) - __d, 70063d1a8abSmrg _Tp(1) - __x); 70163d1a8abSmrg const _Tp __F2 = __hyperg_series(__c - __a, __c - __b, _Tp(1) + __d, 70263d1a8abSmrg _Tp(1) - __x); 70363d1a8abSmrg 70463d1a8abSmrg const _Tp __F = __pre1 * __F1 + __pre2 * __F2; 70563d1a8abSmrg 70663d1a8abSmrg return __F; 70763d1a8abSmrg } 70863d1a8abSmrg } 70963d1a8abSmrg 71063d1a8abSmrg 71163d1a8abSmrg /** 71263d1a8abSmrg * @brief Return the hypogeometric function @f$ _2F_1(a,b;c;x) @f$. 71363d1a8abSmrg * 71463d1a8abSmrg * The hypogeometric function is defined by 71563d1a8abSmrg * @f[ 71663d1a8abSmrg * _2F_1(a,b;c;x) = \frac{\Gamma(c)}{\Gamma(a)\Gamma(b)} 71763d1a8abSmrg * \sum_{n=0}^{\infty} 71863d1a8abSmrg * \frac{\Gamma(a+n)\Gamma(b+n)}{\Gamma(c+n)} 71963d1a8abSmrg * \frac{x^n}{n!} 72063d1a8abSmrg * @f] 72163d1a8abSmrg * 72263d1a8abSmrg * @param __a The first @a numerator parameter. 72363d1a8abSmrg * @param __a The second @a numerator parameter. 72463d1a8abSmrg * @param __c The @a denominator parameter. 72563d1a8abSmrg * @param __x The argument of the confluent hypergeometric function. 72663d1a8abSmrg * @return The confluent hypergeometric function. 72763d1a8abSmrg */ 72863d1a8abSmrg template<typename _Tp> 72963d1a8abSmrg _Tp __hyperg(_Tp __a,_Tp __b,_Tp __c,_Tp __x)73063d1a8abSmrg __hyperg(_Tp __a, _Tp __b, _Tp __c, _Tp __x) 73163d1a8abSmrg { 73263d1a8abSmrg #if _GLIBCXX_USE_C99_MATH_TR1 73363d1a8abSmrg const _Tp __a_nint = _GLIBCXX_MATH_NS::nearbyint(__a); 73463d1a8abSmrg const _Tp __b_nint = _GLIBCXX_MATH_NS::nearbyint(__b); 73563d1a8abSmrg const _Tp __c_nint = _GLIBCXX_MATH_NS::nearbyint(__c); 73663d1a8abSmrg #else 73763d1a8abSmrg const _Tp __a_nint = static_cast<int>(__a + _Tp(0.5L)); 73863d1a8abSmrg const _Tp __b_nint = static_cast<int>(__b + _Tp(0.5L)); 73963d1a8abSmrg const _Tp __c_nint = static_cast<int>(__c + _Tp(0.5L)); 74063d1a8abSmrg #endif 74163d1a8abSmrg const _Tp __toler = _Tp(1000) * std::numeric_limits<_Tp>::epsilon(); 74263d1a8abSmrg if (std::abs(__x) >= _Tp(1)) 74363d1a8abSmrg std::__throw_domain_error(__N("Argument outside unit circle " 74463d1a8abSmrg "in __hyperg.")); 74563d1a8abSmrg else if (__isnan(__a) || __isnan(__b) 74663d1a8abSmrg || __isnan(__c) || __isnan(__x)) 74763d1a8abSmrg return std::numeric_limits<_Tp>::quiet_NaN(); 74863d1a8abSmrg else if (__c_nint == __c && __c_nint <= _Tp(0)) 74963d1a8abSmrg return std::numeric_limits<_Tp>::infinity(); 75063d1a8abSmrg else if (std::abs(__c - __b) < __toler || std::abs(__c - __a) < __toler) 75163d1a8abSmrg return std::pow(_Tp(1) - __x, __c - __a - __b); 75263d1a8abSmrg else if (__a >= _Tp(0) && __b >= _Tp(0) && __c >= _Tp(0) 75363d1a8abSmrg && __x >= _Tp(0) && __x < _Tp(0.995L)) 75463d1a8abSmrg return __hyperg_series(__a, __b, __c, __x); 75563d1a8abSmrg else if (std::abs(__a) < _Tp(10) && std::abs(__b) < _Tp(10)) 75663d1a8abSmrg { 75763d1a8abSmrg // For integer a and b the hypergeometric function is a 75863d1a8abSmrg // finite polynomial. 75963d1a8abSmrg if (__a < _Tp(0) && std::abs(__a - __a_nint) < __toler) 76063d1a8abSmrg return __hyperg_series(__a_nint, __b, __c, __x); 76163d1a8abSmrg else if (__b < _Tp(0) && std::abs(__b - __b_nint) < __toler) 76263d1a8abSmrg return __hyperg_series(__a, __b_nint, __c, __x); 76363d1a8abSmrg else if (__x < -_Tp(0.25L)) 76463d1a8abSmrg return __hyperg_luke(__a, __b, __c, __x); 76563d1a8abSmrg else if (__x < _Tp(0.5L)) 76663d1a8abSmrg return __hyperg_series(__a, __b, __c, __x); 76763d1a8abSmrg else 76863d1a8abSmrg if (std::abs(__c) > _Tp(10)) 76963d1a8abSmrg return __hyperg_series(__a, __b, __c, __x); 77063d1a8abSmrg else 77163d1a8abSmrg return __hyperg_reflect(__a, __b, __c, __x); 77263d1a8abSmrg } 77363d1a8abSmrg else 77463d1a8abSmrg return __hyperg_luke(__a, __b, __c, __x); 77563d1a8abSmrg } 77663d1a8abSmrg } // namespace __detail 77763d1a8abSmrg #undef _GLIBCXX_MATH_NS 77863d1a8abSmrg #if ! _GLIBCXX_USE_STD_SPEC_FUNCS && defined(_GLIBCXX_TR1_CMATH) 77963d1a8abSmrg } // namespace tr1 78063d1a8abSmrg #endif 781c7a68eb7Smrg 782c7a68eb7Smrg _GLIBCXX_END_NAMESPACE_VERSION 78363d1a8abSmrg } 78463d1a8abSmrg 78563d1a8abSmrg #endif // _GLIBCXX_TR1_HYPERGEOMETRIC_TCC 786