26 |
#include <gsl/gsl_sf_exp.h> |
#include <gsl/gsl_sf_exp.h> |
27 |
#include <gsl/gsl_sf_log.h> |
#include <gsl/gsl_sf_log.h> |
28 |
#include <gsl/gsl_sf_gamma.h> |
#include <gsl/gsl_sf_gamma.h> |
29 |
|
#include <gsl/gsl_sf_expint.h> |
30 |
|
|
31 |
#include "error.h" |
#include "error.h" |
32 |
|
|
47 |
return GSL_SUCCESS; |
return GSL_SUCCESS; |
48 |
} |
} |
49 |
else { |
else { |
|
double mu = (x-a)/a; |
|
|
double term1; |
|
50 |
gsl_sf_result gstar; |
gsl_sf_result gstar; |
51 |
gsl_sf_result ln_term; |
gsl_sf_result ln_term; |
52 |
gsl_sf_log_1plusx_mx_e(mu, &ln_term); /* log(1+mu) - mu */ |
double term1; |
53 |
|
if (x < a) { |
54 |
|
double u = x/a; |
55 |
|
ln_term.val = log(u) - u + 1.0; |
56 |
|
ln_term.err = ln_term.val * GSL_DBL_EPSILON; |
57 |
|
} else { |
58 |
|
double mu = (x-a)/a; |
59 |
|
gsl_sf_log_1plusx_mx_e(mu, &ln_term); /* log(1+mu) - mu */ |
60 |
|
}; |
61 |
gsl_sf_gammastar_e(a, &gstar); |
gsl_sf_gammastar_e(a, &gstar); |
62 |
term1 = exp(a*ln_term.val)/sqrt(2.0*M_PI*a); |
term1 = exp(a*ln_term.val)/sqrt(2.0*M_PI*a); |
63 |
result->val = term1/gstar.val; |
result->val = term1/gstar.val; |
65 |
result->err += gstar.err/fabs(gstar.val) * fabs(result->val); |
result->err += gstar.err/fabs(gstar.val) * fabs(result->val); |
66 |
return GSL_SUCCESS; |
return GSL_SUCCESS; |
67 |
} |
} |
68 |
|
|
69 |
} |
} |
70 |
|
|
71 |
|
|
128 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
129 |
|
|
130 |
if(n == nmax) |
if(n == nmax) |
131 |
GSL_ERROR ("error", GSL_EMAXITER); |
GSL_ERROR ("error in large x asymptotic", GSL_EMAXITER); |
132 |
else |
else |
133 |
return stat_D; |
return stat_D; |
134 |
} |
} |
177 |
} |
} |
178 |
|
|
179 |
|
|
180 |
/* Continued fraction for Q. |
/* Continued fraction which occurs in evaluation |
181 |
* |
* of Q(a,x) or Gamma(a,x). |
|
* Q(a,x) = D(a,x) a/x F(a,x) |
|
|
* 1 (1-a)/x 1/x (2-a)/x 2/x (3-a)/x |
|
|
* F(a,x) = ---- ------- ----- -------- ----- -------- ... |
|
|
* 1 + 1 + 1 + 1 + 1 + 1 + |
|
182 |
* |
* |
183 |
* Hans E. Plesser, 2002-01-22 (hans dot plesser at itf dot nlh dot no): |
* 1 (1-a)/x 1/x (2-a)/x 2/x (3-a)/x |
184 |
|
* F(a,x) = ---- ------- ----- -------- ----- -------- ... |
185 |
|
* 1 + 1 + 1 + 1 + 1 + 1 + |
186 |
* |
* |
187 |
* Since the Gautschi equivalent series method for CF evaluation may lead |
* Hans E. Plesser, 2002-01-22 (hans dot plesser at itf dot nlh dot no). |
|
* to singularities, I have replaced it with the modified Lentz algorithm |
|
|
* given in |
|
188 |
* |
* |
189 |
* I J Thompson and A R Barnett |
* Split out from gamma_inc_Q_CF() by GJ [Tue Apr 1 13:16:41 MST 2003]. |
190 |
* Coulomb and Bessel Functions of Complex Arguments and Order |
* See gamma_inc_Q_CF() below. |
|
* J Computational Physics 64:490-509 (1986) |
|
|
* |
|
|
* In consequence, gamma_inc_Q_CF_protected() is now obsolete and has been |
|
|
* removed. |
|
|
* |
|
|
* Identification of terms between the above equation for F(a, x) and |
|
|
* the first equation in the appendix of Thompson&Barnett is as follows: |
|
|
* |
|
|
* b_0 = 0, b_n = 1 for all n > 0 |
|
|
* |
|
|
* a_1 = 1 |
|
|
* a_n = (n/2-a)/x for n even |
|
|
* a_n = (n-1)/(2x) for n odd |
|
191 |
* |
* |
192 |
*/ |
*/ |
193 |
|
static int |
194 |
static |
gamma_inc_F_CF(const double a, const double x, gsl_sf_result * result) |
|
int |
|
|
gamma_inc_Q_CF(const double a, const double x, gsl_sf_result * result) |
|
195 |
{ |
{ |
196 |
const int nmax = 5000; |
const int nmax = 5000; |
197 |
const double small = gsl_pow_3 (GSL_DBL_EPSILON); |
const double small = gsl_pow_3 (GSL_DBL_EPSILON); |
198 |
|
|
|
gsl_sf_result D; |
|
|
const int stat_D = gamma_inc_D(a, x, &D); |
|
|
|
|
199 |
double hn = 1.0; /* convergent */ |
double hn = 1.0; /* convergent */ |
200 |
double Cn = 1.0 / small; |
double Cn = 1.0 / small; |
201 |
double Dn = 1.0; |
double Dn = 1.0; |
203 |
|
|
204 |
/* n == 1 has a_1, b_1, b_0 independent of a,x, |
/* n == 1 has a_1, b_1, b_0 independent of a,x, |
205 |
so that has been done by hand */ |
so that has been done by hand */ |
206 |
for ( n = 2 ; n < nmax ; n++ ) { |
for ( n = 2 ; n < nmax ; n++ ) |
207 |
|
{ |
208 |
double an; |
double an; |
209 |
double delta; |
double delta; |
210 |
|
|
222 |
Dn = 1.0 / Dn; |
Dn = 1.0 / Dn; |
223 |
delta = Cn * Dn; |
delta = Cn * Dn; |
224 |
hn *= delta; |
hn *= delta; |
225 |
if(fabs(delta-1) < GSL_DBL_EPSILON) break; |
if(fabs(delta-1.0) < GSL_DBL_EPSILON) break; |
226 |
} |
} |
227 |
|
|
228 |
result->val = D.val * (a/x) * hn; |
result->val = hn; |
229 |
result->err = D.err * fabs((a/x) * hn); |
result->err = 2.0*GSL_DBL_EPSILON * fabs(hn); |
230 |
result->err += GSL_DBL_EPSILON * (2.0 + 0.5*n) * fabs(result->val); |
result->err += GSL_DBL_EPSILON * (2.0 + 0.5*n) * fabs(result->val); |
231 |
|
|
232 |
if(n == nmax) |
if(n == nmax) |
233 |
GSL_ERROR ("error", GSL_EMAXITER); |
GSL_ERROR ("error in CF for F(a,x)", GSL_EMAXITER); |
234 |
else |
else |
235 |
return stat_D; |
return GSL_SUCCESS; |
236 |
|
} |
237 |
|
|
238 |
|
|
239 |
|
/* Continued fraction for Q. |
240 |
|
* |
241 |
|
* Q(a,x) = D(a,x) a/x F(a,x) |
242 |
|
* |
243 |
|
* Hans E. Plesser, 2002-01-22 (hans dot plesser at itf dot nlh dot no): |
244 |
|
* |
245 |
|
* Since the Gautschi equivalent series method for CF evaluation may lead |
246 |
|
* to singularities, I have replaced it with the modified Lentz algorithm |
247 |
|
* given in |
248 |
|
* |
249 |
|
* I J Thompson and A R Barnett |
250 |
|
* Coulomb and Bessel Functions of Complex Arguments and Order |
251 |
|
* J Computational Physics 64:490-509 (1986) |
252 |
|
* |
253 |
|
* In consequence, gamma_inc_Q_CF_protected() is now obsolete and has been |
254 |
|
* removed. |
255 |
|
* |
256 |
|
* Identification of terms between the above equation for F(a, x) and |
257 |
|
* the first equation in the appendix of Thompson&Barnett is as follows: |
258 |
|
* |
259 |
|
* b_0 = 0, b_n = 1 for all n > 0 |
260 |
|
* |
261 |
|
* a_1 = 1 |
262 |
|
* a_n = (n/2-a)/x for n even |
263 |
|
* a_n = (n-1)/(2x) for n odd |
264 |
|
* |
265 |
|
*/ |
266 |
|
static |
267 |
|
int |
268 |
|
gamma_inc_Q_CF(const double a, const double x, gsl_sf_result * result) |
269 |
|
{ |
270 |
|
gsl_sf_result D; |
271 |
|
gsl_sf_result F; |
272 |
|
const int stat_D = gamma_inc_D(a, x, &D); |
273 |
|
const int stat_F = gamma_inc_F_CF(a, x, &F); |
274 |
|
|
275 |
|
result->val = D.val * (a/x) * F.val; |
276 |
|
result->err = D.err * fabs((a/x) * F.val) + fabs(D.val * a/x * F.err); |
277 |
|
|
278 |
|
return GSL_ERROR_SELECT_2(stat_F, stat_D); |
279 |
} |
} |
280 |
|
|
281 |
|
|
282 |
/* Useful for small a and x. Handles the subtraction analytically. |
/* Useful for small a and x. Handles the subtraction analytically. |
283 |
*/ |
*/ |
284 |
static |
static |
302 |
const double c4 = -0.04166666666666666667 |
const double c4 = -0.04166666666666666667 |
303 |
* (-1.758243446661483480 + lnx) |
* (-1.758243446661483480 + lnx) |
304 |
* (-0.764428657272716373 + lnx) |
* (-0.764428657272716373 + lnx) |
305 |
* ( 0.723980571623507657 + lnx) |
* ( 0.723980571623507657 + lnx) |
306 |
* ( 4.107554191916823640 + lnx); |
* ( 4.107554191916823640 + lnx); |
307 |
const double c5 = -0.0083333333333333333 |
const double c5 = -0.0083333333333333333 |
308 |
* (-2.06563396085715900 + lnx) |
* (-2.06563396085715900 + lnx) |
309 |
* (-1.28459889470864700 + lnx) |
* (-1.28459889470864700 + lnx) |
310 |
* (-0.27583535756454143 + lnx) |
* (-0.27583535756454143 + lnx) |
311 |
* ( 1.33677371336239618 + lnx) |
* ( 1.33677371336239618 + lnx) |
312 |
* ( 5.17537282427561550 + lnx); |
* ( 5.17537282427561550 + lnx); |
313 |
const double c6 = -0.0013888888888888889 |
const double c6 = -0.0013888888888888889 |
314 |
* (-2.30814336454783200 + lnx) |
* (-2.30814336454783200 + lnx) |
315 |
* (-1.65846557706987300 + lnx) |
* (-1.65846557706987300 + lnx) |
387 |
} |
} |
388 |
|
|
389 |
|
|
390 |
|
/* series for small a and x, but not defined for a == 0 */ |
391 |
|
static int |
392 |
|
gamma_inc_series(double a, double x, gsl_sf_result * result) |
393 |
|
{ |
394 |
|
gsl_sf_result Q; |
395 |
|
gsl_sf_result G; |
396 |
|
const int stat_Q = gamma_inc_Q_series(a, x, &Q); |
397 |
|
const int stat_G = gsl_sf_gamma_e(a, &G); |
398 |
|
result->val = Q.val * G.val; |
399 |
|
result->err = fabs(Q.val * G.err) + fabs(Q.err * G.val); |
400 |
|
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
401 |
|
|
402 |
|
return GSL_ERROR_SELECT_2(stat_Q, stat_G); |
403 |
|
} |
404 |
|
|
405 |
|
|
406 |
|
static int |
407 |
|
gamma_inc_a_gt_0(double a, double x, gsl_sf_result * result) |
408 |
|
{ |
409 |
|
/* x > 0 and a > 0; use result for Q */ |
410 |
|
gsl_sf_result Q; |
411 |
|
gsl_sf_result G; |
412 |
|
const int stat_Q = gsl_sf_gamma_inc_Q_e(a, x, &Q); |
413 |
|
const int stat_G = gsl_sf_gamma_e(a, &G); |
414 |
|
|
415 |
|
result->val = G.val * Q.val; |
416 |
|
result->err = fabs(G.val * Q.err) + fabs(G.err * Q.val); |
417 |
|
result->err += 2.0*GSL_DBL_EPSILON * fabs(result->val); |
418 |
|
|
419 |
|
return GSL_ERROR_SELECT_2(stat_G, stat_Q); |
420 |
|
} |
421 |
|
|
422 |
|
|
423 |
|
static int |
424 |
|
gamma_inc_CF(double a, double x, gsl_sf_result * result) |
425 |
|
{ |
426 |
|
gsl_sf_result F; |
427 |
|
gsl_sf_result pre; |
428 |
|
const int stat_F = gamma_inc_F_CF(a, x, &F); |
429 |
|
const int stat_E = gsl_sf_exp_e((a-1.0)*log(x) - x, &pre); |
430 |
|
|
431 |
|
result->val = F.val * pre.val; |
432 |
|
result->err = fabs(F.err * pre.val) + fabs(F.val * pre.err); |
433 |
|
result->err += (2.0 + fabs(a)) * GSL_DBL_EPSILON * fabs(result->val); |
434 |
|
|
435 |
|
return GSL_ERROR_SELECT_2(stat_F, stat_E); |
436 |
|
} |
437 |
|
|
438 |
|
|
439 |
|
/* evaluate Gamma(0,x), x > 0 */ |
440 |
|
#define GAMMA_INC_A_0(x, result) gsl_sf_expint_E1_e(x, result) |
441 |
|
|
442 |
|
|
443 |
/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ |
/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ |
444 |
|
|
445 |
int |
int |
446 |
gsl_sf_gamma_inc_Q_e(const double a, const double x, gsl_sf_result * result) |
gsl_sf_gamma_inc_Q_e(const double a, const double x, gsl_sf_result * result) |
447 |
{ |
{ |
448 |
if(a <= 0.0 || x < 0.0) { |
if(a < 0.0 || x < 0.0) { |
449 |
DOMAIN_ERROR(result); |
DOMAIN_ERROR(result); |
450 |
} |
} |
451 |
else if(x == 0.0) { |
else if(x == 0.0) { |
453 |
result->err = 0.0; |
result->err = 0.0; |
454 |
return GSL_SUCCESS; |
return GSL_SUCCESS; |
455 |
} |
} |
456 |
|
else if(a == 0.0) |
457 |
|
{ |
458 |
|
result->val = 0.0; |
459 |
|
result->err = 0.0; |
460 |
|
return GSL_SUCCESS; |
461 |
|
} |
462 |
else if(x <= 0.5*a) { |
else if(x <= 0.5*a) { |
463 |
/* If the series is quick, do that. It is |
/* If the series is quick, do that. It is |
464 |
* robust and simple. |
* robust and simple. |
588 |
} |
} |
589 |
|
|
590 |
|
|
591 |
|
int |
592 |
|
gsl_sf_gamma_inc_e(const double a, const double x, gsl_sf_result * result) |
593 |
|
{ |
594 |
|
if(x < 0.0) { |
595 |
|
DOMAIN_ERROR(result); |
596 |
|
} |
597 |
|
else if(x == 0.0) { |
598 |
|
return gsl_sf_gamma_e(a, result); |
599 |
|
} |
600 |
|
else if(a == 0.0) |
601 |
|
{ |
602 |
|
return GAMMA_INC_A_0(x, result); |
603 |
|
} |
604 |
|
else if(a > 0.0) |
605 |
|
{ |
606 |
|
return gamma_inc_a_gt_0(a, x, result); |
607 |
|
} |
608 |
|
else if(x > 0.25) |
609 |
|
{ |
610 |
|
/* continued fraction seems to fail for x too small; otherwise |
611 |
|
it is ok, independent of the value of |x/a|, because of the |
612 |
|
non-oscillation in the expansion, i.e. the CF is |
613 |
|
un-conditionally convergent for a < 0 and x > 0 |
614 |
|
*/ |
615 |
|
return gamma_inc_CF(a, x, result); |
616 |
|
} |
617 |
|
else if(fabs(a) < 0.5) |
618 |
|
{ |
619 |
|
return gamma_inc_series(a, x, result); |
620 |
|
} |
621 |
|
else |
622 |
|
{ |
623 |
|
/* a = fa + da; da >= 0 */ |
624 |
|
const double fa = floor(a); |
625 |
|
const double da = a - fa; |
626 |
|
|
627 |
|
gsl_sf_result g_da; |
628 |
|
const int stat_g_da = ( da > 0.0 ? gamma_inc_a_gt_0(da, x, &g_da) |
629 |
|
: GAMMA_INC_A_0(x, &g_da)); |
630 |
|
|
631 |
|
double alpha = da; |
632 |
|
double gax = g_da.val; |
633 |
|
|
634 |
|
/* Gamma(alpha-1,x) = 1/(alpha-1) (Gamma(a,x) - x^(alpha-1) e^-x) */ |
635 |
|
do |
636 |
|
{ |
637 |
|
const double shift = exp(-x + (alpha-1.0)*log(x)); |
638 |
|
gax = (gax - shift) / (alpha - 1.0); |
639 |
|
alpha -= 1.0; |
640 |
|
} while(alpha > a); |
641 |
|
|
642 |
|
result->val = gax; |
643 |
|
result->err = 2.0*(1.0 + fabs(a))*GSL_DBL_EPSILON*fabs(gax); |
644 |
|
return stat_g_da; |
645 |
|
} |
646 |
|
|
647 |
|
} |
648 |
|
|
649 |
|
|
650 |
/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/ |
/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/ |
651 |
|
|
652 |
#include "eval.h" |
#include "eval.h" |
660 |
{ |
{ |
661 |
EVAL_RESULT(gsl_sf_gamma_inc_Q_e(a, x, &result)); |
EVAL_RESULT(gsl_sf_gamma_inc_Q_e(a, x, &result)); |
662 |
} |
} |
663 |
|
|
664 |
|
double gsl_sf_gamma_inc(const double a, const double x) |
665 |
|
{ |
666 |
|
EVAL_RESULT(gsl_sf_gamma_inc_e(a, x, &result)); |
667 |
|
} |