]>
Commit | Line | Data |
---|---|---|
b34f60ac | 1 | // Special functions -*- C++ -*- |
2 | ||
f1717362 | 3 | // Copyright (C) 2006-2016 Free Software Foundation, Inc. |
b34f60ac | 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 | |
6bc9506f | 8 | // Free Software Foundation; either version 3, or (at your option) |
b34f60ac | 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 | // | |
6bc9506f | 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/>. | |
b34f60ac | 24 | |
25 | /** @file tr1/hypergeometric.tcc | |
26 | * This is an internal header file, included by other library headers. | |
5846aeac | 27 | * Do not attempt to use it directly. @headername{tr1/cmath} |
b34f60ac | 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 | ||
c17b0a1c | 41 | #ifndef _GLIBCXX_TR1_HYPERGEOMETRIC_TCC |
42 | #define _GLIBCXX_TR1_HYPERGEOMETRIC_TCC 1 | |
b34f60ac | 43 | |
2948dd21 | 44 | namespace std _GLIBCXX_VISIBILITY(default) |
b34f60ac | 45 | { |
c17b0a1c | 46 | namespace tr1 |
47 | { | |
b34f60ac | 48 | // [5.2] Special functions |
49 | ||
b34f60ac | 50 | // Implementation-space details. |
b34f60ac | 51 | namespace __detail |
52 | { | |
2948dd21 | 53 | _GLIBCXX_BEGIN_NAMESPACE_VERSION |
b34f60ac | 54 | |
55 | /** | |
9fc1117c | 56 | * @brief This routine returns the confluent hypergeometric function |
b34f60ac | 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 | * | |
72117d76 | 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. | |
b34f60ac | 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 | |
cd7f5f45 | 78 | __conf_hyperg_series(_Tp __a, _Tp __c, _Tp __x) |
b34f60ac | 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 | |
cd7f5f45 | 115 | __conf_hyperg_luke(_Tp __a, _Tp __c, _Tp __xin) |
b34f60ac | 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 | * | |
72117d76 | 215 | * @param __a The @a numerator parameter. |
216 | * @param __c The @a denominator parameter. | |
b34f60ac | 217 | * @param __x The argument of the confluent hypergeometric function. |
218 | * @return The confluent hypergeometric function. | |
219 | */ | |
220 | template<typename _Tp> | |
cd7f5f45 | 221 | _Tp |
222 | __conf_hyperg(_Tp __a, _Tp __c, _Tp __x) | |
b34f60ac | 223 | { |
224 | #if _GLIBCXX_USE_C99_MATH_TR1 | |
c17b0a1c | 225 | const _Tp __c_nint = std::tr1::nearbyint(__c); |
b34f60ac | 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 | * | |
72117d76 | 258 | * @param __a The first @a numerator parameter. |
259 | * @param __a The second @a numerator parameter. | |
260 | * @param __c The @a denominator parameter. | |
b34f60ac | 261 | * @param __x The argument of the confluent hypergeometric function. |
262 | * @return The confluent hypergeometric function. | |
263 | */ | |
264 | template<typename _Tp> | |
265 | _Tp | |
cd7f5f45 | 266 | __hyperg_series(_Tp __a, _Tp __b, _Tp __c, _Tp __x) |
b34f60ac | 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<typename _Tp> | |
298 | _Tp | |
cd7f5f45 | 299 | __hyperg_luke(_Tp __a, _Tp __b, _Tp __c, _Tp __xin) |
b34f60ac | 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 | /** | |
72117d76 | 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. | |
b34f60ac | 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<typename _Tp> | |
432 | _Tp | |
cd7f5f45 | 433 | __hyperg_reflect(_Tp __a, _Tp __b, _Tp __c, _Tp __x) |
b34f60ac | 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; | |
b1a09c64 | 473 | __try |
b34f60ac | 474 | { |
475 | __lng_ad = __log_gamma(__ad); | |
476 | __lng_ad1 = __log_gamma(__a + __d1); | |
477 | __lng_bd1 = __log_gamma(__b + __d1); | |
478 | } | |
b1a09c64 | 479 | __catch(...) |
b34f60ac | 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) | |
72117d76 | 505 | std::__throw_runtime_error(__N("Overflow of gamma functions" |
506 | " in __hyperg_luke.")); | |
b34f60ac | 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; | |
b1a09c64 | 521 | __try |
b34f60ac | 522 | { |
523 | __lng_ad2 = __log_gamma(__a + __d2); | |
524 | __lng_bd2 = __log_gamma(__b + __d2); | |
525 | } | |
b1a09c64 | 526 | __catch(...) |
b34f60ac | 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 | { | |
72117d76 | 552 | // Values for psi functions use recurrence; |
553 | // Abramowitz & Stegun 6.3.5 | |
b34f60ac | 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); | |
b1a09c64 | 597 | __try |
b34f60ac | 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 | } | |
b1a09c64 | 604 | __catch(...) |
b34f60ac | 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); | |
b1a09c64 | 612 | __try |
b34f60ac | 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 | } | |
b1a09c64 | 619 | __catch(...) |
b34f60ac | 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 | ||
048ff85f | 704 | /** |
b34f60ac | 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 | * | |
72117d76 | 715 | * @param __a The first @a numerator parameter. |
716 | * @param __a The second @a numerator parameter. | |
717 | * @param __c The @a denominator parameter. | |
b34f60ac | 718 | * @param __x The argument of the confluent hypergeometric function. |
719 | * @return The confluent hypergeometric function. | |
720 | */ | |
721 | template<typename _Tp> | |
cd7f5f45 | 722 | _Tp |
723 | __hyperg(_Tp __a, _Tp __b, _Tp __c, _Tp __x) | |
b34f60ac | 724 | { |
725 | #if _GLIBCXX_USE_C99_MATH_TR1 | |
c17b0a1c | 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); | |
b34f60ac | 729 | #else |
730 | const _Tp __a_nint = static_cast<int>(__a + _Tp(0.5L)); | |
731 | const _Tp __b_nint = static_cast<int>(__b + _Tp(0.5L)); | |
732 | const _Tp __c_nint = static_cast<int>(__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 | { | |
72117d76 | 750 | // For integer a and b the hypergeometric function is a |
751 | // finite polynomial. | |
b34f60ac | 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 | ||
2948dd21 | 770 | _GLIBCXX_END_NAMESPACE_VERSION |
b34f60ac | 771 | } // namespace std::tr1::__detail |
c17b0a1c | 772 | } |
b34f60ac | 773 | } |
774 | ||
c17b0a1c | 775 | #endif // _GLIBCXX_TR1_HYPERGEOMETRIC_TCC |