]>
Commit | Line | Data |
---|---|---|
7c62b943 BK |
1 | // Special functions -*- C++ -*- |
2 | ||
bc2631e0 | 3 | // Copyright (C) 2006, 2007, 2008, 2009 |
7c62b943 BK |
4 | // Free Software Foundation, Inc. |
5 | // | |
6 | // This file is part of the GNU ISO C++ Library. This library is free | |
7 | // software; you can redistribute it and/or modify it under the | |
8 | // terms of the GNU General Public License as published by the | |
748086b7 | 9 | // Free Software Foundation; either version 3, or (at your option) |
7c62b943 BK |
10 | // any later version. |
11 | // | |
12 | // This library is distributed in the hope that it will be useful, | |
13 | // but WITHOUT ANY WARRANTY; without even the implied warranty of | |
14 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
15 | // GNU General Public License for more details. | |
16 | // | |
748086b7 JJ |
17 | // Under Section 7 of GPL version 3, you are granted additional |
18 | // permissions described in the GCC Runtime Library Exception, version | |
19 | // 3.1, as published by the Free Software Foundation. | |
20 | ||
21 | // You should have received a copy of the GNU General Public License and | |
22 | // a copy of the GCC Runtime Library Exception along with this program; | |
23 | // see the files COPYING3 and COPYING.RUNTIME respectively. If not, see | |
24 | // <http://www.gnu.org/licenses/>. | |
7c62b943 BK |
25 | |
26 | /** @file tr1/hypergeometric.tcc | |
27 | * This is an internal header file, included by other library headers. | |
28 | * You should not attempt to use it directly. | |
29 | */ | |
30 | ||
31 | // | |
32 | // ISO C++ 14882 TR1: 5.2 Special functions | |
33 | // | |
34 | ||
35 | // Written by Edward Smith-Rowland based: | |
36 | // (1) Handbook of Mathematical Functions, | |
37 | // ed. Milton Abramowitz and Irene A. Stegun, | |
38 | // Dover Publications, | |
39 | // Section 6, pp. 555-566 | |
40 | // (2) The Gnu Scientific Library, http://www.gnu.org/software/gsl | |
41 | ||
e133ace8 PC |
42 | #ifndef _GLIBCXX_TR1_HYPERGEOMETRIC_TCC |
43 | #define _GLIBCXX_TR1_HYPERGEOMETRIC_TCC 1 | |
7c62b943 BK |
44 | |
45 | namespace std | |
46 | { | |
e133ace8 PC |
47 | namespace tr1 |
48 | { | |
7c62b943 BK |
49 | |
50 | // [5.2] Special functions | |
51 | ||
7c62b943 | 52 | // Implementation-space details. |
7c62b943 BK |
53 | namespace __detail |
54 | { | |
55 | ||
56 | /** | |
28dac70a | 57 | * @brief This routine returns the confluent hypergeometric function |
7c62b943 BK |
58 | * by series expansion. |
59 | * | |
60 | * @f[ | |
61 | * _1F_1(a;c;x) = \frac{\Gamma(c)}{\Gamma(a)} | |
62 | * \sum_{n=0}^{\infty} | |
63 | * \frac{\Gamma(a+n)}{\Gamma(c+n)} | |
64 | * \frac{x^n}{n!} | |
65 | * @f] | |
66 | * | |
67 | * If a and b are integers and a < 0 and either b > 0 or b < a then the | |
68 | * series is a polynomial with a finite number of terms. If b is an integer | |
28dac70a | 69 | * and b <= 0 the confluent hypergeometric function is undefined. |
7c62b943 BK |
70 | * |
71 | * @param __a The "numerator" parameter. | |
72 | * @param __c The "denominator" parameter. | |
73 | * @param __x The argument of the confluent hypergeometric function. | |
74 | * @return The confluent hypergeometric function. | |
75 | */ | |
76 | template<typename _Tp> | |
77 | _Tp | |
78 | __conf_hyperg_series(const _Tp __a, const _Tp __c, const _Tp __x) | |
79 | { | |
80 | const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); | |
81 | ||
82 | _Tp __term = _Tp(1); | |
83 | _Tp __Fac = _Tp(1); | |
84 | const unsigned int __max_iter = 100000; | |
85 | unsigned int __i; | |
86 | for (__i = 0; __i < __max_iter; ++__i) | |
87 | { | |
88 | __term *= (__a + _Tp(__i)) * __x | |
89 | / ((__c + _Tp(__i)) * _Tp(1 + __i)); | |
90 | if (std::abs(__term) < __eps) | |
91 | { | |
92 | break; | |
93 | } | |
94 | __Fac += __term; | |
95 | } | |
96 | if (__i == __max_iter) | |
97 | std::__throw_runtime_error(__N("Series failed to converge " | |
98 | "in __conf_hyperg_series.")); | |
99 | ||
100 | return __Fac; | |
101 | } | |
102 | ||
103 | ||
104 | /** | |
105 | * @brief Return the hypogeometric function @f$ _2F_1(a,b;c;x) @f$ | |
106 | * by an iterative procedure described in | |
107 | * Luke, Algorithms for the Computation of Mathematical Functions. | |
108 | * | |
109 | * Like the case of the 2F1 rational approximations, these are | |
110 | * probably guaranteed to converge for x < 0, barring gross | |
111 | * numerical instability in the pre-asymptotic regime. | |
112 | */ | |
113 | template<typename _Tp> | |
114 | _Tp | |
115 | __conf_hyperg_luke(const _Tp __a, const _Tp __c, const _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 "numerator" parameter. | |
216 | * @param __c The "denominator" parameter. | |
217 | * @param __x The argument of the confluent hypergeometric function. | |
218 | * @return The confluent hypergeometric function. | |
219 | */ | |
220 | template<typename _Tp> | |
221 | inline _Tp | |
222 | __conf_hyperg(const _Tp __a, const _Tp __c, const _Tp __x) | |
223 | { | |
224 | #if _GLIBCXX_USE_C99_MATH_TR1 | |
e133ace8 | 225 | const _Tp __c_nint = std::tr1::nearbyint(__c); |
7c62b943 BK |
226 | #else |
227 | const _Tp __c_nint = static_cast<int>(__c + _Tp(0.5L)); | |
228 | #endif | |
229 | if (__isnan(__a) || __isnan(__c) || __isnan(__x)) | |
230 | return std::numeric_limits<_Tp>::quiet_NaN(); | |
231 | else if (__c_nint == __c && __c_nint <= 0) | |
232 | return std::numeric_limits<_Tp>::infinity(); | |
233 | else if (__a == _Tp(0)) | |
234 | return _Tp(1); | |
235 | else if (__c == __a) | |
236 | return std::exp(__x); | |
237 | else if (__x < _Tp(0)) | |
238 | return __conf_hyperg_luke(__a, __c, __x); | |
239 | else | |
240 | return __conf_hyperg_series(__a, __c, __x); | |
241 | } | |
242 | ||
243 | ||
244 | /** | |
245 | * @brief Return the hypogeometric function @f$ _2F_1(a,b;c;x) @f$ | |
246 | * by series expansion. | |
247 | * | |
248 | * The hypogeometric function is defined by | |
249 | * @f[ | |
250 | * _2F_1(a,b;c;x) = \frac{\Gamma(c)}{\Gamma(a)\Gamma(b)} | |
251 | * \sum_{n=0}^{\infty} | |
252 | * \frac{\Gamma(a+n)\Gamma(b+n)}{\Gamma(c+n)} | |
253 | * \frac{x^n}{n!} | |
254 | * @f] | |
255 | * | |
256 | * This works and it's pretty fast. | |
257 | * | |
258 | * @param __a The first "numerator" parameter. | |
259 | * @param __a The second "numerator" parameter. | |
260 | * @param __c The "denominator" parameter. | |
261 | * @param __x The argument of the confluent hypergeometric function. | |
262 | * @return The confluent hypergeometric function. | |
263 | */ | |
264 | template<typename _Tp> | |
265 | _Tp | |
266 | __hyperg_series(const _Tp __a, const _Tp __b, | |
267 | const _Tp __c, const _Tp __x) | |
268 | { | |
269 | const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); | |
270 | ||
271 | _Tp __term = _Tp(1); | |
272 | _Tp __Fabc = _Tp(1); | |
273 | const unsigned int __max_iter = 100000; | |
274 | unsigned int __i; | |
275 | for (__i = 0; __i < __max_iter; ++__i) | |
276 | { | |
277 | __term *= (__a + _Tp(__i)) * (__b + _Tp(__i)) * __x | |
278 | / ((__c + _Tp(__i)) * _Tp(1 + __i)); | |
279 | if (std::abs(__term) < __eps) | |
280 | { | |
281 | break; | |
282 | } | |
283 | __Fabc += __term; | |
284 | } | |
285 | if (__i == __max_iter) | |
286 | std::__throw_runtime_error(__N("Series failed to converge " | |
287 | "in __hyperg_series.")); | |
288 | ||
289 | return __Fabc; | |
290 | } | |
291 | ||
292 | ||
293 | /** | |
294 | * @brief Return the hypogeometric function @f$ _2F_1(a,b;c;x) @f$ | |
295 | * by an iterative procedure described in | |
296 | * Luke, Algorithms for the Computation of Mathematical Functions. | |
297 | */ | |
298 | template<typename _Tp> | |
299 | _Tp | |
300 | __hyperg_luke(const _Tp __a, const _Tp __b, const _Tp __c, | |
301 | const _Tp __xin) | |
302 | { | |
303 | const _Tp __big = std::pow(std::numeric_limits<_Tp>::max(), _Tp(0.16L)); | |
304 | const int __nmax = 20000; | |
305 | const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); | |
306 | const _Tp __x = -__xin; | |
307 | const _Tp __x3 = __x * __x * __x; | |
308 | const _Tp __t0 = __a * __b / __c; | |
309 | const _Tp __t1 = (__a + _Tp(1)) * (__b + _Tp(1)) / (_Tp(2) * __c); | |
310 | const _Tp __t2 = (__a + _Tp(2)) * (__b + _Tp(2)) | |
311 | / (_Tp(2) * (__c + _Tp(1))); | |
312 | ||
313 | _Tp __F = _Tp(1); | |
314 | ||
315 | _Tp __Bnm3 = _Tp(1); | |
316 | _Tp __Bnm2 = _Tp(1) + __t1 * __x; | |
317 | _Tp __Bnm1 = _Tp(1) + __t2 * __x * (_Tp(1) + __t1 / _Tp(3) * __x); | |
318 | ||
319 | _Tp __Anm3 = _Tp(1); | |
320 | _Tp __Anm2 = __Bnm2 - __t0 * __x; | |
321 | _Tp __Anm1 = __Bnm1 - __t0 * (_Tp(1) + __t2 * __x) * __x | |
322 | + __t0 * __t1 * (__c / (__c + _Tp(1))) * __x * __x; | |
323 | ||
324 | int __n = 3; | |
325 | while (1) | |
326 | { | |
327 | const _Tp __npam1 = _Tp(__n - 1) + __a; | |
328 | const _Tp __npbm1 = _Tp(__n - 1) + __b; | |
329 | const _Tp __npcm1 = _Tp(__n - 1) + __c; | |
330 | const _Tp __npam2 = _Tp(__n - 2) + __a; | |
331 | const _Tp __npbm2 = _Tp(__n - 2) + __b; | |
332 | const _Tp __npcm2 = _Tp(__n - 2) + __c; | |
333 | const _Tp __tnm1 = _Tp(2 * __n - 1); | |
334 | const _Tp __tnm3 = _Tp(2 * __n - 3); | |
335 | const _Tp __tnm5 = _Tp(2 * __n - 5); | |
336 | const _Tp __n2 = __n * __n; | |
337 | const _Tp __F1 = (_Tp(3) * __n2 + (__a + __b - _Tp(6)) * __n | |
338 | + _Tp(2) - __a * __b - _Tp(2) * (__a + __b)) | |
339 | / (_Tp(2) * __tnm3 * __npcm1); | |
340 | const _Tp __F2 = -(_Tp(3) * __n2 - (__a + __b + _Tp(6)) * __n | |
341 | + _Tp(2) - __a * __b) * __npam1 * __npbm1 | |
342 | / (_Tp(4) * __tnm1 * __tnm3 * __npcm2 * __npcm1); | |
343 | const _Tp __F3 = (__npam2 * __npam1 * __npbm2 * __npbm1 | |
344 | * (_Tp(__n - 2) - __a) * (_Tp(__n - 2) - __b)) | |
345 | / (_Tp(8) * __tnm3 * __tnm3 * __tnm5 | |
346 | * (_Tp(__n - 3) + __c) * __npcm2 * __npcm1); | |
347 | const _Tp __E = -__npam1 * __npbm1 * (_Tp(__n - 1) - __c) | |
348 | / (_Tp(2) * __tnm3 * __npcm2 * __npcm1); | |
349 | ||
350 | _Tp __An = (_Tp(1) + __F1 * __x) * __Anm1 | |
351 | + (__E + __F2 * __x) * __x * __Anm2 + __F3 * __x3 * __Anm3; | |
352 | _Tp __Bn = (_Tp(1) + __F1 * __x) * __Bnm1 | |
353 | + (__E + __F2 * __x) * __x * __Bnm2 + __F3 * __x3 * __Bnm3; | |
354 | const _Tp __r = __An / __Bn; | |
355 | ||
356 | const _Tp __prec = std::abs((__F - __r) / __F); | |
357 | __F = __r; | |
358 | ||
359 | if (__prec < __eps || __n > __nmax) | |
360 | break; | |
361 | ||
362 | if (std::abs(__An) > __big || std::abs(__Bn) > __big) | |
363 | { | |
364 | __An /= __big; | |
365 | __Bn /= __big; | |
366 | __Anm1 /= __big; | |
367 | __Bnm1 /= __big; | |
368 | __Anm2 /= __big; | |
369 | __Bnm2 /= __big; | |
370 | __Anm3 /= __big; | |
371 | __Bnm3 /= __big; | |
372 | } | |
373 | else if (std::abs(__An) < _Tp(1) / __big | |
374 | || std::abs(__Bn) < _Tp(1) / __big) | |
375 | { | |
376 | __An *= __big; | |
377 | __Bn *= __big; | |
378 | __Anm1 *= __big; | |
379 | __Bnm1 *= __big; | |
380 | __Anm2 *= __big; | |
381 | __Bnm2 *= __big; | |
382 | __Anm3 *= __big; | |
383 | __Bnm3 *= __big; | |
384 | } | |
385 | ||
386 | ++__n; | |
387 | __Bnm3 = __Bnm2; | |
388 | __Bnm2 = __Bnm1; | |
389 | __Bnm1 = __Bn; | |
390 | __Anm3 = __Anm2; | |
391 | __Anm2 = __Anm1; | |
392 | __Anm1 = __An; | |
393 | } | |
394 | ||
395 | if (__n >= __nmax) | |
396 | std::__throw_runtime_error(__N("Iteration failed to converge " | |
397 | "in __hyperg_luke.")); | |
398 | ||
399 | return __F; | |
400 | } | |
401 | ||
402 | ||
403 | /** | |
404 | * @brief Return the hypogeometric function @f$ _2F_1(a,b;c;x) @f$ by the reflection | |
405 | * formulae in Abramowitz & Stegun formula 15.3.6 for d = c - a - b not integral | |
406 | * and formula 15.3.11 for d = c - a - b integral. | |
407 | * This assumes a, b, c != negative integer. | |
408 | * | |
409 | * The hypogeometric function is defined by | |
410 | * @f[ | |
411 | * _2F_1(a,b;c;x) = \frac{\Gamma(c)}{\Gamma(a)\Gamma(b)} | |
412 | * \sum_{n=0}^{\infty} | |
413 | * \frac{\Gamma(a+n)\Gamma(b+n)}{\Gamma(c+n)} | |
414 | * \frac{x^n}{n!} | |
415 | * @f] | |
416 | * | |
417 | * The reflection formula for nonintegral @f$ d = c - a - b @f$ is: | |
418 | * @f[ | |
419 | * _2F_1(a,b;c;x) = \frac{\Gamma(c)\Gamma(d)}{\Gamma(c-a)\Gamma(c-b)} | |
420 | * _2F_1(a,b;1-d;1-x) | |
421 | * + \frac{\Gamma(c)\Gamma(-d)}{\Gamma(a)\Gamma(b)} | |
422 | * _2F_1(c-a,c-b;1+d;1-x) | |
423 | * @f] | |
424 | * | |
425 | * The reflection formula for integral @f$ m = c - a - b @f$ is: | |
426 | * @f[ | |
427 | * _2F_1(a,b;a+b+m;x) = \frac{\Gamma(m)\Gamma(a+b+m)}{\Gamma(a+m)\Gamma(b+m)} | |
428 | * \sum_{k=0}^{m-1} \frac{(m+a)_k(m+b)_k}{k!(1-m)_k} | |
429 | * - | |
430 | * @f] | |
431 | */ | |
432 | template<typename _Tp> | |
433 | _Tp | |
434 | __hyperg_reflect(const _Tp __a, const _Tp __b, const _Tp __c, | |
435 | const _Tp __x) | |
436 | { | |
437 | const _Tp __d = __c - __a - __b; | |
438 | const int __intd = std::floor(__d + _Tp(0.5L)); | |
439 | const _Tp __eps = std::numeric_limits<_Tp>::epsilon(); | |
440 | const _Tp __toler = _Tp(1000) * __eps; | |
441 | const _Tp __log_max = std::log(std::numeric_limits<_Tp>::max()); | |
442 | const bool __d_integer = (std::abs(__d - __intd) < __toler); | |
443 | ||
444 | if (__d_integer) | |
445 | { | |
446 | const _Tp __ln_omx = std::log(_Tp(1) - __x); | |
447 | const _Tp __ad = std::abs(__d); | |
448 | _Tp __F1, __F2; | |
449 | ||
450 | _Tp __d1, __d2; | |
451 | if (__d >= _Tp(0)) | |
452 | { | |
453 | __d1 = __d; | |
454 | __d2 = _Tp(0); | |
455 | } | |
456 | else | |
457 | { | |
458 | __d1 = _Tp(0); | |
459 | __d2 = __d; | |
460 | } | |
461 | ||
462 | const _Tp __lng_c = __log_gamma(__c); | |
463 | ||
464 | // Evaluate F1. | |
465 | if (__ad < __eps) | |
466 | { | |
467 | // d = c - a - b = 0. | |
468 | __F1 = _Tp(0); | |
469 | } | |
470 | else | |
471 | { | |
472 | ||
473 | bool __ok_d1 = true; | |
474 | _Tp __lng_ad, __lng_ad1, __lng_bd1; | |
bc2631e0 | 475 | __try |
7c62b943 BK |
476 | { |
477 | __lng_ad = __log_gamma(__ad); | |
478 | __lng_ad1 = __log_gamma(__a + __d1); | |
479 | __lng_bd1 = __log_gamma(__b + __d1); | |
480 | } | |
bc2631e0 | 481 | __catch(...) |
7c62b943 BK |
482 | { |
483 | __ok_d1 = false; | |
484 | } | |
485 | ||
486 | if (__ok_d1) | |
487 | { | |
488 | /* Gamma functions in the denominator are ok. | |
489 | * Proceed with evaluation. | |
490 | */ | |
491 | _Tp __sum1 = _Tp(1); | |
492 | _Tp __term = _Tp(1); | |
493 | _Tp __ln_pre1 = __lng_ad + __lng_c + __d2 * __ln_omx | |
494 | - __lng_ad1 - __lng_bd1; | |
495 | ||
496 | /* Do F1 sum. | |
497 | */ | |
498 | for (int __i = 1; __i < __ad; ++__i) | |
499 | { | |
500 | const int __j = __i - 1; | |
501 | __term *= (__a + __d2 + __j) * (__b + __d2 + __j) | |
502 | / (_Tp(1) + __d2 + __j) / __i * (_Tp(1) - __x); | |
503 | __sum1 += __term; | |
504 | } | |
505 | ||
506 | if (__ln_pre1 > __log_max) | |
507 | std::__throw_runtime_error(__N("Overflow of gamma functions " | |
508 | "in __hyperg_luke.")); | |
509 | else | |
510 | __F1 = std::exp(__ln_pre1) * __sum1; | |
511 | } | |
512 | else | |
513 | { | |
514 | // Gamma functions in the denominator were not ok. | |
515 | // So the F1 term is zero. | |
516 | __F1 = _Tp(0); | |
517 | } | |
518 | } // end F1 evaluation | |
519 | ||
520 | // Evaluate F2. | |
521 | bool __ok_d2 = true; | |
522 | _Tp __lng_ad2, __lng_bd2; | |
bc2631e0 | 523 | __try |
7c62b943 BK |
524 | { |
525 | __lng_ad2 = __log_gamma(__a + __d2); | |
526 | __lng_bd2 = __log_gamma(__b + __d2); | |
527 | } | |
bc2631e0 | 528 | __catch(...) |
7c62b943 BK |
529 | { |
530 | __ok_d2 = false; | |
531 | } | |
532 | ||
533 | if (__ok_d2) | |
534 | { | |
535 | // Gamma functions in the denominator are ok. | |
536 | // Proceed with evaluation. | |
537 | const int __maxiter = 2000; | |
538 | const _Tp __psi_1 = -__numeric_constants<_Tp>::__gamma_e(); | |
539 | const _Tp __psi_1pd = __psi(_Tp(1) + __ad); | |
540 | const _Tp __psi_apd1 = __psi(__a + __d1); | |
541 | const _Tp __psi_bpd1 = __psi(__b + __d1); | |
542 | ||
543 | _Tp __psi_term = __psi_1 + __psi_1pd - __psi_apd1 | |
544 | - __psi_bpd1 - __ln_omx; | |
545 | _Tp __fact = _Tp(1); | |
546 | _Tp __sum2 = __psi_term; | |
547 | _Tp __ln_pre2 = __lng_c + __d1 * __ln_omx | |
548 | - __lng_ad2 - __lng_bd2; | |
549 | ||
550 | // Do F2 sum. | |
551 | int __j; | |
552 | for (__j = 1; __j < __maxiter; ++__j) | |
553 | { | |
554 | // Values for psi functions use recurrence; Abramowitz & Stegun 6.3.5 | |
555 | const _Tp __term1 = _Tp(1) / _Tp(__j) | |
556 | + _Tp(1) / (__ad + __j); | |
557 | const _Tp __term2 = _Tp(1) / (__a + __d1 + _Tp(__j - 1)) | |
558 | + _Tp(1) / (__b + __d1 + _Tp(__j - 1)); | |
559 | __psi_term += __term1 - __term2; | |
560 | __fact *= (__a + __d1 + _Tp(__j - 1)) | |
561 | * (__b + __d1 + _Tp(__j - 1)) | |
562 | / ((__ad + __j) * __j) * (_Tp(1) - __x); | |
563 | const _Tp __delta = __fact * __psi_term; | |
564 | __sum2 += __delta; | |
565 | if (std::abs(__delta) < __eps * std::abs(__sum2)) | |
566 | break; | |
567 | } | |
568 | if (__j == __maxiter) | |
569 | std::__throw_runtime_error(__N("Sum F2 failed to converge " | |
570 | "in __hyperg_reflect")); | |
571 | ||
572 | if (__sum2 == _Tp(0)) | |
573 | __F2 = _Tp(0); | |
574 | else | |
575 | __F2 = std::exp(__ln_pre2) * __sum2; | |
576 | } | |
577 | else | |
578 | { | |
579 | // Gamma functions in the denominator not ok. | |
580 | // So the F2 term is zero. | |
581 | __F2 = _Tp(0); | |
582 | } // end F2 evaluation | |
583 | ||
584 | const _Tp __sgn_2 = (__intd % 2 == 1 ? -_Tp(1) : _Tp(1)); | |
585 | const _Tp __F = __F1 + __sgn_2 * __F2; | |
586 | ||
587 | return __F; | |
588 | } | |
589 | else | |
590 | { | |
591 | // d = c - a - b not an integer. | |
592 | ||
593 | // These gamma functions appear in the denominator, so we | |
594 | // catch their harmless domain errors and set the terms to zero. | |
595 | bool __ok1 = true; | |
596 | _Tp __sgn_g1ca = _Tp(0), __ln_g1ca = _Tp(0); | |
597 | _Tp __sgn_g1cb = _Tp(0), __ln_g1cb = _Tp(0); | |
bc2631e0 | 598 | __try |
7c62b943 BK |
599 | { |
600 | __sgn_g1ca = __log_gamma_sign(__c - __a); | |
601 | __ln_g1ca = __log_gamma(__c - __a); | |
602 | __sgn_g1cb = __log_gamma_sign(__c - __b); | |
603 | __ln_g1cb = __log_gamma(__c - __b); | |
604 | } | |
bc2631e0 | 605 | __catch(...) |
7c62b943 BK |
606 | { |
607 | __ok1 = false; | |
608 | } | |
609 | ||
610 | bool __ok2 = true; | |
611 | _Tp __sgn_g2a = _Tp(0), __ln_g2a = _Tp(0); | |
612 | _Tp __sgn_g2b = _Tp(0), __ln_g2b = _Tp(0); | |
bc2631e0 | 613 | __try |
7c62b943 BK |
614 | { |
615 | __sgn_g2a = __log_gamma_sign(__a); | |
616 | __ln_g2a = __log_gamma(__a); | |
617 | __sgn_g2b = __log_gamma_sign(__b); | |
618 | __ln_g2b = __log_gamma(__b); | |
619 | } | |
bc2631e0 | 620 | __catch(...) |
7c62b943 BK |
621 | { |
622 | __ok2 = false; | |
623 | } | |
624 | ||
625 | const _Tp __sgn_gc = __log_gamma_sign(__c); | |
626 | const _Tp __ln_gc = __log_gamma(__c); | |
627 | const _Tp __sgn_gd = __log_gamma_sign(__d); | |
628 | const _Tp __ln_gd = __log_gamma(__d); | |
629 | const _Tp __sgn_gmd = __log_gamma_sign(-__d); | |
630 | const _Tp __ln_gmd = __log_gamma(-__d); | |
631 | ||
632 | const _Tp __sgn1 = __sgn_gc * __sgn_gd * __sgn_g1ca * __sgn_g1cb; | |
633 | const _Tp __sgn2 = __sgn_gc * __sgn_gmd * __sgn_g2a * __sgn_g2b; | |
634 | ||
635 | _Tp __pre1, __pre2; | |
636 | if (__ok1 && __ok2) | |
637 | { | |
638 | _Tp __ln_pre1 = __ln_gc + __ln_gd - __ln_g1ca - __ln_g1cb; | |
639 | _Tp __ln_pre2 = __ln_gc + __ln_gmd - __ln_g2a - __ln_g2b | |
640 | + __d * std::log(_Tp(1) - __x); | |
641 | if (__ln_pre1 < __log_max && __ln_pre2 < __log_max) | |
642 | { | |
643 | __pre1 = std::exp(__ln_pre1); | |
644 | __pre2 = std::exp(__ln_pre2); | |
645 | __pre1 *= __sgn1; | |
646 | __pre2 *= __sgn2; | |
647 | } | |
648 | else | |
649 | { | |
650 | std::__throw_runtime_error(__N("Overflow of gamma functions " | |
651 | "in __hyperg_reflect")); | |
652 | } | |
653 | } | |
654 | else if (__ok1 && !__ok2) | |
655 | { | |
656 | _Tp __ln_pre1 = __ln_gc + __ln_gd - __ln_g1ca - __ln_g1cb; | |
657 | if (__ln_pre1 < __log_max) | |
658 | { | |
659 | __pre1 = std::exp(__ln_pre1); | |
660 | __pre1 *= __sgn1; | |
661 | __pre2 = _Tp(0); | |
662 | } | |
663 | else | |
664 | { | |
665 | std::__throw_runtime_error(__N("Overflow of gamma functions " | |
666 | "in __hyperg_reflect")); | |
667 | } | |
668 | } | |
669 | else if (!__ok1 && __ok2) | |
670 | { | |
671 | _Tp __ln_pre2 = __ln_gc + __ln_gmd - __ln_g2a - __ln_g2b | |
672 | + __d * std::log(_Tp(1) - __x); | |
673 | if (__ln_pre2 < __log_max) | |
674 | { | |
675 | __pre1 = _Tp(0); | |
676 | __pre2 = std::exp(__ln_pre2); | |
677 | __pre2 *= __sgn2; | |
678 | } | |
679 | else | |
680 | { | |
681 | std::__throw_runtime_error(__N("Overflow of gamma functions " | |
682 | "in __hyperg_reflect")); | |
683 | } | |
684 | } | |
685 | else | |
686 | { | |
687 | __pre1 = _Tp(0); | |
688 | __pre2 = _Tp(0); | |
689 | std::__throw_runtime_error(__N("Underflow of gamma functions " | |
690 | "in __hyperg_reflect")); | |
691 | } | |
692 | ||
693 | const _Tp __F1 = __hyperg_series(__a, __b, _Tp(1) - __d, | |
694 | _Tp(1) - __x); | |
695 | const _Tp __F2 = __hyperg_series(__c - __a, __c - __b, _Tp(1) + __d, | |
696 | _Tp(1) - __x); | |
697 | ||
698 | const _Tp __F = __pre1 * __F1 + __pre2 * __F2; | |
699 | ||
700 | return __F; | |
701 | } | |
702 | } | |
703 | ||
704 | ||
6165bbdd | 705 | /** |
7c62b943 BK |
706 | * @brief Return the hypogeometric function @f$ _2F_1(a,b;c;x) @f$. |
707 | * | |
708 | * The hypogeometric function is defined by | |
709 | * @f[ | |
710 | * _2F_1(a,b;c;x) = \frac{\Gamma(c)}{\Gamma(a)\Gamma(b)} | |
711 | * \sum_{n=0}^{\infty} | |
712 | * \frac{\Gamma(a+n)\Gamma(b+n)}{\Gamma(c+n)} | |
713 | * \frac{x^n}{n!} | |
714 | * @f] | |
715 | * | |
716 | * @param __a The first "numerator" parameter. | |
717 | * @param __a The second "numerator" parameter. | |
718 | * @param __c The "denominator" parameter. | |
719 | * @param __x The argument of the confluent hypergeometric function. | |
720 | * @return The confluent hypergeometric function. | |
721 | */ | |
722 | template<typename _Tp> | |
723 | inline _Tp | |
724 | __hyperg(const _Tp __a, const _Tp __b, const _Tp __c, const _Tp __x) | |
725 | { | |
726 | #if _GLIBCXX_USE_C99_MATH_TR1 | |
e133ace8 PC |
727 | const _Tp __a_nint = std::tr1::nearbyint(__a); |
728 | const _Tp __b_nint = std::tr1::nearbyint(__b); | |
729 | const _Tp __c_nint = std::tr1::nearbyint(__c); | |
7c62b943 BK |
730 | #else |
731 | const _Tp __a_nint = static_cast<int>(__a + _Tp(0.5L)); | |
732 | const _Tp __b_nint = static_cast<int>(__b + _Tp(0.5L)); | |
733 | const _Tp __c_nint = static_cast<int>(__c + _Tp(0.5L)); | |
734 | #endif | |
735 | const _Tp __toler = _Tp(1000) * std::numeric_limits<_Tp>::epsilon(); | |
736 | if (std::abs(__x) >= _Tp(1)) | |
737 | std::__throw_domain_error(__N("Argument outside unit circle " | |
738 | "in __hyperg.")); | |
739 | else if (__isnan(__a) || __isnan(__b) | |
740 | || __isnan(__c) || __isnan(__x)) | |
741 | return std::numeric_limits<_Tp>::quiet_NaN(); | |
742 | else if (__c_nint == __c && __c_nint <= _Tp(0)) | |
743 | return std::numeric_limits<_Tp>::infinity(); | |
744 | else if (std::abs(__c - __b) < __toler || std::abs(__c - __a) < __toler) | |
745 | return std::pow(_Tp(1) - __x, __c - __a - __b); | |
746 | else if (__a >= _Tp(0) && __b >= _Tp(0) && __c >= _Tp(0) | |
747 | && __x >= _Tp(0) && __x < _Tp(0.995L)) | |
748 | return __hyperg_series(__a, __b, __c, __x); | |
749 | else if (std::abs(__a) < _Tp(10) && std::abs(__b) < _Tp(10)) | |
750 | { | |
751 | // For integer a and b the hypergeometric function is a 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 | } // namespace std::tr1::__detail | |
e133ace8 | 771 | } |
7c62b943 BK |
772 | } |
773 | ||
e133ace8 | 774 | #endif // _GLIBCXX_TR1_HYPERGEOMETRIC_TCC |