| 5134 |
serge |
1 |
// Special functions -*- C++ -*- |
|
|
2 |
|
|
|
3 |
// Copyright (C) 2006-2013 Free Software Foundation, Inc. |
|
|
4 |
// |
|
|
5 |
// This file is part of the GNU ISO C++ Library. This library is free |
|
|
6 |
// software; you can redistribute it and/or modify it under the |
|
|
7 |
// terms of the GNU General Public License as published by the |
|
|
8 |
// Free Software Foundation; either version 3, or (at your option) |
|
|
9 |
// any later version. |
|
|
10 |
// |
|
|
11 |
// This library is distributed in the hope that it will be useful, |
|
|
12 |
// but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
13 |
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
|
14 |
// GNU General Public License for more details. |
|
|
15 |
// |
|
|
16 |
// Under Section 7 of GPL version 3, you are granted additional |
|
|
17 |
// permissions described in the GCC Runtime Library Exception, version |
|
|
18 |
// 3.1, as published by the Free Software Foundation. |
|
|
19 |
|
|
|
20 |
// You should have received a copy of the GNU General Public License and |
|
|
21 |
// a copy of the GCC Runtime Library Exception along with this program; |
|
|
22 |
// see the files COPYING3 and COPYING.RUNTIME respectively. If not, see |
|
|
23 |
// . |
|
|
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 |
|
|
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 |
|
|
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 |
|
|
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(__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 |
|
|
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 |
|
|
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 |
|
|
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 |
|
|
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(__a + _Tp(0.5L)); |
|
|
731 |
const _Tp __b_nint = static_cast(__b + _Tp(0.5L)); |
|
|
732 |
const _Tp __c_nint = static_cast(__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 |