]>
Commit | Line | Data |
---|---|---|
1 | // Special functions -*- C++ -*- | |
2 | ||
3 | // Copyright (C) 2006-2024 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 |