1 /* specfunc/gamma_inc.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18 */
19
20 /* Author: G. Jungman */
21
22 #include <config.h>
23 #include <gsl/gsl_math.h>
24 #include <gsl/gsl_errno.h>
25 #include <gsl/gsl_sf_erf.h>
26 #include <gsl/gsl_sf_exp.h>
27 #include <gsl/gsl_sf_log.h>
28 #include <gsl/gsl_sf_gamma.h>
29 #include <gsl/gsl_sf_expint.h>
30
31 #include "error.h"
32
33 /* The dominant part,
34 * D(a,x) := x^a e^(-x) / Gamma(a+1)
35 */
36 static
37 int
gamma_inc_D(const double a,const double x,gsl_sf_result * result)38 gamma_inc_D(const double a, const double x, gsl_sf_result * result)
39 {
40 if(a < 10.0) {
41 double lnr;
42 gsl_sf_result lg;
43 gsl_sf_lngamma_e(a+1.0, &lg);
44 lnr = a * log(x) - x - lg.val;
45 result->val = exp(lnr);
46 result->err = 2.0 * GSL_DBL_EPSILON * (fabs(lnr) + 1.0) * fabs(result->val);
47 return GSL_SUCCESS;
48 }
49 else {
50 gsl_sf_result gstar;
51 gsl_sf_result ln_term;
52 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);
62 term1 = exp(a*ln_term.val)/sqrt(2.0*M_PI*a);
63 result->val = term1/gstar.val;
64 result->err = 2.0 * GSL_DBL_EPSILON * (fabs(a*ln_term.val) + 1.0) * fabs(result->val);
65 result->err += gstar.err/fabs(gstar.val) * fabs(result->val);
66 return GSL_SUCCESS;
67 }
68
69 }
70
71
72 /* P series representation.
73 */
74 static
75 int
gamma_inc_P_series(const double a,const double x,gsl_sf_result * result)76 gamma_inc_P_series(const double a, const double x, gsl_sf_result * result)
77 {
78 const int nmax = 5000;
79
80 gsl_sf_result D;
81 int stat_D = gamma_inc_D(a, x, &D);
82
83 double sum = 1.0;
84 double term = 1.0;
85 int n;
86 for(n=1; n<nmax; n++) {
87 term *= x/(a+n);
88 sum += term;
89 if(fabs(term/sum) < GSL_DBL_EPSILON) break;
90 }
91
92 result->val = D.val * sum;
93 result->err = D.err * fabs(sum);
94 result->err += (1.0 + n) * GSL_DBL_EPSILON * fabs(result->val);
95
96 if(n == nmax)
97 GSL_ERROR ("error", GSL_EMAXITER);
98 else
99 return stat_D;
100 }
101
102
103 /* Q large x asymptotic
104 */
105 static
106 int
gamma_inc_Q_large_x(const double a,const double x,gsl_sf_result * result)107 gamma_inc_Q_large_x(const double a, const double x, gsl_sf_result * result)
108 {
109 const int nmax = 5000;
110
111 gsl_sf_result D;
112 const int stat_D = gamma_inc_D(a, x, &D);
113
114 double sum = 1.0;
115 double term = 1.0;
116 double last = 1.0;
117 int n;
118 for(n=1; n<nmax; n++) {
119 term *= (a-n)/x;
120 if(fabs(term/last) > 1.0) break;
121 if(fabs(term/sum) < GSL_DBL_EPSILON) break;
122 sum += term;
123 last = term;
124 }
125
126 result->val = D.val * (a/x) * sum;
127 result->err = D.err * fabs((a/x) * sum);
128 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
129
130 if(n == nmax)
131 GSL_ERROR ("error in large x asymptotic", GSL_EMAXITER);
132 else
133 return stat_D;
134 }
135
136
137 /* Uniform asymptotic for x near a, a and x large.
138 * See [Temme, p. 285]
139 * FIXME: need c1 coefficient
140 */
141 static
142 int
gamma_inc_Q_asymp_unif(const double a,const double x,gsl_sf_result * result)143 gamma_inc_Q_asymp_unif(const double a, const double x, gsl_sf_result * result)
144 {
145 const double rta = sqrt(a);
146 const double eps = (x-a)/a;
147
148 gsl_sf_result ln_term;
149 const int stat_ln = gsl_sf_log_1plusx_mx_e(eps, &ln_term); /* log(1+eps) - eps */
150 const double eta = eps * sqrt(-2.0*ln_term.val/(eps*eps));
151
152 gsl_sf_result erfc;
153
154 double R;
155 double c0, c1;
156
157 gsl_sf_erfc_e(eta*M_SQRT2*rta, &erfc);
158
159 if(fabs(eps) < GSL_ROOT5_DBL_EPSILON) {
160 c0 = -1.0/3.0 + eps*(1.0/12.0 - eps*(23.0/540.0 - eps*(353.0/12960.0 - eps*589.0/30240.0)));
161 c1 = 0.0;
162 }
163 else {
164 double rt_term;
165 rt_term = sqrt(-2.0 * ln_term.val/(eps*eps));
166 c0 = (1.0 - 1.0/rt_term)/eps;
167 c1 = 0.0;
168 }
169
170 R = exp(-0.5*a*eta*eta)/(M_SQRT2*M_SQRTPI*rta) * (c0 + c1/a);
171
172 result->val = 0.5 * erfc.val + R;
173 result->err = GSL_DBL_EPSILON * fabs(R * 0.5 * a*eta*eta) + 0.5 * erfc.err;
174 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
175
176 return stat_ln;
177 }
178
179
180 /* Continued fraction which occurs in evaluation
181 * of Q(a,x) or Gamma(a,x).
182 *
183 * 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 * Hans E. Plesser, 2002-01-22 (hans dot plesser at itf dot nlh dot no).
188 *
189 * Split out from gamma_inc_Q_CF() by GJ [Tue Apr 1 13:16:41 MST 2003].
190 * See gamma_inc_Q_CF() below.
191 *
192 */
193 static int
gamma_inc_F_CF(const double a,const double x,gsl_sf_result * result)194 gamma_inc_F_CF(const double a, const double x, gsl_sf_result * result)
195 {
196 const int nmax = 5000;
197 const double small = gsl_pow_3 (GSL_DBL_EPSILON);
198
199 double hn = 1.0; /* convergent */
200 double Cn = 1.0 / small;
201 double Dn = 1.0;
202 int n;
203
204 /* n == 1 has a_1, b_1, b_0 independent of a,x,
205 so that has been done by hand */
206 for ( n = 2 ; n < nmax ; n++ )
207 {
208 double an;
209 double delta;
210
211 if(GSL_IS_ODD(n))
212 an = 0.5*(n-1)/x;
213 else
214 an = (0.5*n-a)/x;
215
216 Dn = 1.0 + an * Dn;
217 if ( fabs(Dn) < small )
218 Dn = small;
219 Cn = 1.0 + an/Cn;
220 if ( fabs(Cn) < small )
221 Cn = small;
222 Dn = 1.0 / Dn;
223 delta = Cn * Dn;
224 hn *= delta;
225 if(fabs(delta-1.0) < GSL_DBL_EPSILON) break;
226 }
227
228 result->val = hn;
229 result->err = 2.0*GSL_DBL_EPSILON * fabs(hn);
230 result->err += GSL_DBL_EPSILON * (2.0 + 0.5*n) * fabs(result->val);
231
232 if(n == nmax)
233 GSL_ERROR ("error in CF for F(a,x)", GSL_EMAXITER);
234 else
235 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
gamma_inc_Q_CF(const double a,const double x,gsl_sf_result * result)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.
283 */
284 static
285 int
gamma_inc_Q_series(const double a,const double x,gsl_sf_result * result)286 gamma_inc_Q_series(const double a, const double x, gsl_sf_result * result)
287 {
288 double term1; /* 1 - x^a/Gamma(a+1) */
289 double sum; /* 1 + (a+1)/(a+2)(-x)/2! + (a+1)/(a+3)(-x)^2/3! + ... */
290 int stat_sum;
291 double term2; /* a temporary variable used at the end */
292
293 {
294 /* Evaluate series for 1 - x^a/Gamma(a+1), small a
295 */
296 const double pg21 = -2.404113806319188570799476; /* PolyGamma[2,1] */
297 const double lnx = log(x);
298 const double el = M_EULER+lnx;
299 const double c1 = -el;
300 const double c2 = M_PI*M_PI/12.0 - 0.5*el*el;
301 const double c3 = el*(M_PI*M_PI/12.0 - el*el/6.0) + pg21/6.0;
302 const double c4 = -0.04166666666666666667
303 * (-1.758243446661483480 + lnx)
304 * (-0.764428657272716373 + lnx)
305 * ( 0.723980571623507657 + lnx)
306 * ( 4.107554191916823640 + lnx);
307 const double c5 = -0.0083333333333333333
308 * (-2.06563396085715900 + lnx)
309 * (-1.28459889470864700 + lnx)
310 * (-0.27583535756454143 + lnx)
311 * ( 1.33677371336239618 + lnx)
312 * ( 5.17537282427561550 + lnx);
313 const double c6 = -0.0013888888888888889
314 * (-2.30814336454783200 + lnx)
315 * (-1.65846557706987300 + lnx)
316 * (-0.88768082560020400 + lnx)
317 * ( 0.17043847751371778 + lnx)
318 * ( 1.92135970115863890 + lnx)
319 * ( 6.22578557795474900 + lnx);
320 const double c7 = -0.00019841269841269841
321 * (-2.5078657901291800 + lnx)
322 * (-1.9478900888958200 + lnx)
323 * (-1.3194837322612730 + lnx)
324 * (-0.5281322700249279 + lnx)
325 * ( 0.5913834939078759 + lnx)
326 * ( 2.4876819633378140 + lnx)
327 * ( 7.2648160783762400 + lnx);
328 const double c8 = -0.00002480158730158730
329 * (-2.677341544966400 + lnx)
330 * (-2.182810448271700 + lnx)
331 * (-1.649350342277400 + lnx)
332 * (-1.014099048290790 + lnx)
333 * (-0.191366955370652 + lnx)
334 * ( 0.995403817918724 + lnx)
335 * ( 3.041323283529310 + lnx)
336 * ( 8.295966556941250 + lnx);
337 const double c9 = -2.75573192239859e-6
338 * (-2.8243487670469080 + lnx)
339 * (-2.3798494322701120 + lnx)
340 * (-1.9143674728689960 + lnx)
341 * (-1.3814529102920370 + lnx)
342 * (-0.7294312810261694 + lnx)
343 * ( 0.1299079285269565 + lnx)
344 * ( 1.3873333251885240 + lnx)
345 * ( 3.5857258865210760 + lnx)
346 * ( 9.3214237073814600 + lnx);
347 const double c10 = -2.75573192239859e-7
348 * (-2.9540329644556910 + lnx)
349 * (-2.5491366926991850 + lnx)
350 * (-2.1348279229279880 + lnx)
351 * (-1.6741881076349450 + lnx)
352 * (-1.1325949616098420 + lnx)
353 * (-0.4590034650618494 + lnx)
354 * ( 0.4399352987435699 + lnx)
355 * ( 1.7702236517651670 + lnx)
356 * ( 4.1231539047474080 + lnx)
357 * ( 10.342627908148680 + lnx);
358
359 term1 = a*(c1+a*(c2+a*(c3+a*(c4+a*(c5+a*(c6+a*(c7+a*(c8+a*(c9+a*c10)))))))));
360 }
361
362 {
363 /* Evaluate the sum.
364 */
365 const int nmax = 5000;
366 double t = 1.0;
367 int n;
368 sum = 1.0;
369
370 for(n=1; n<nmax; n++) {
371 t *= -x/(n+1.0);
372 sum += (a+1.0)/(a+n+1.0)*t;
373 if(fabs(t/sum) < GSL_DBL_EPSILON) break;
374 }
375
376 if(n == nmax)
377 stat_sum = GSL_EMAXITER;
378 else
379 stat_sum = GSL_SUCCESS;
380 }
381
382 term2 = (1.0 - term1) * a/(a+1.0) * x * sum;
383 result->val = term1 + term2;
384 result->err = GSL_DBL_EPSILON * (fabs(term1) + 2.0*fabs(term2));
385 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
386 return stat_sum;
387 }
388
389
390 /* series for small a and x, but not defined for a == 0 */
391 static int
gamma_inc_series(double a,double x,gsl_sf_result * result)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
gamma_inc_a_gt_0(double a,double x,gsl_sf_result * result)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
gamma_inc_CF(double a,double x,gsl_sf_result * result)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 *-*-*-*-*-*-*-*-*-*-*-*/
444
445 int
gsl_sf_gamma_inc_Q_e(const double a,const double x,gsl_sf_result * result)446 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) {
449 DOMAIN_ERROR(result);
450 }
451 else if(x == 0.0) {
452 result->val = 1.0;
453 result->err = 0.0;
454 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) {
463 /* If the series is quick, do that. It is
464 * robust and simple.
465 */
466 gsl_sf_result P;
467 int stat_P = gamma_inc_P_series(a, x, &P);
468 result->val = 1.0 - P.val;
469 result->err = P.err;
470 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
471 return stat_P;
472 }
473 else if(a >= 1.0e+06 && (x-a)*(x-a) < a) {
474 /* Then try the difficult asymptotic regime.
475 * This is the only way to do this region.
476 */
477 return gamma_inc_Q_asymp_unif(a, x, result);
478 }
479 else if(a < 0.2 && x < 5.0) {
480 /* Cancellations at small a must be handled
481 * analytically; x should not be too big
482 * either since the series terms grow
483 * with x and log(x).
484 */
485 return gamma_inc_Q_series(a, x, result);
486 }
487 else if(a <= x) {
488 if(x <= 1.0e+06) {
489 /* Continued fraction is excellent for x >~ a.
490 * We do not let x be too large when x > a since
491 * it is somewhat pointless to try this there;
492 * the function is rapidly decreasing for
493 * x large and x > a, and it will just
494 * underflow in that region anyway. We
495 * catch that case in the standard
496 * large-x method.
497 */
498 return gamma_inc_Q_CF(a, x, result);
499 }
500 else {
501 return gamma_inc_Q_large_x(a, x, result);
502 }
503 }
504 else {
505 if(a < 0.8*x) {
506 /* Continued fraction again. The convergence
507 * is a little slower here, but that is fine.
508 * We have to trade that off against the slow
509 * convergence of the series, which is the
510 * only other option.
511 */
512 return gamma_inc_Q_CF(a, x, result);
513 }
514 else {
515 gsl_sf_result P;
516 int stat_P = gamma_inc_P_series(a, x, &P);
517 result->val = 1.0 - P.val;
518 result->err = P.err;
519 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
520 return stat_P;
521 }
522 }
523 }
524
525
526 int
gsl_sf_gamma_inc_P_e(const double a,const double x,gsl_sf_result * result)527 gsl_sf_gamma_inc_P_e(const double a, const double x, gsl_sf_result * result)
528 {
529 if(a <= 0.0 || x < 0.0) {
530 DOMAIN_ERROR(result);
531 }
532 else if(x == 0.0) {
533 result->val = 0.0;
534 result->err = 0.0;
535 return GSL_SUCCESS;
536 }
537 else if(x < 20.0 || x < 0.5*a) {
538 /* Do the easy series cases. Robust and quick.
539 */
540 return gamma_inc_P_series(a, x, result);
541 }
542 else if(a > 1.0e+06 && (x-a)*(x-a) < a) {
543 /* Crossover region. Note that Q and P are
544 * roughly the same order of magnitude here,
545 * so the subtraction is stable.
546 */
547 gsl_sf_result Q;
548 int stat_Q = gamma_inc_Q_asymp_unif(a, x, &Q);
549 result->val = 1.0 - Q.val;
550 result->err = Q.err;
551 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
552 return stat_Q;
553 }
554 else if(a <= x) {
555 /* Q <~ P in this area, so the
556 * subtractions are stable.
557 */
558 gsl_sf_result Q;
559 int stat_Q;
560 if(a > 0.2*x) {
561 stat_Q = gamma_inc_Q_CF(a, x, &Q);
562 }
563 else {
564 stat_Q = gamma_inc_Q_large_x(a, x, &Q);
565 }
566 result->val = 1.0 - Q.val;
567 result->err = Q.err;
568 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
569 return stat_Q;
570 }
571 else {
572 if((x-a)*(x-a) < a) {
573 /* This condition is meant to insure
574 * that Q is not very close to 1,
575 * so the subtraction is stable.
576 */
577 gsl_sf_result Q;
578 int stat_Q = gamma_inc_Q_CF(a, x, &Q);
579 result->val = 1.0 - Q.val;
580 result->err = Q.err;
581 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
582 return stat_Q;
583 }
584 else {
585 return gamma_inc_P_series(a, x, result);
586 }
587 }
588 }
589
590
591 int
gsl_sf_gamma_inc_e(const double a,const double x,gsl_sf_result * result)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 *-*-*-*-*-*-*-*-*-*-*/
651
652 #include "eval.h"
653
gsl_sf_gamma_inc_P(const double a,const double x)654 double gsl_sf_gamma_inc_P(const double a, const double x)
655 {
656 EVAL_RESULT(gsl_sf_gamma_inc_P_e(a, x, &result));
657 }
658
gsl_sf_gamma_inc_Q(const double a,const double x)659 double gsl_sf_gamma_inc_Q(const double a, const double x)
660 {
661 EVAL_RESULT(gsl_sf_gamma_inc_Q_e(a, x, &result));
662 }
663
gsl_sf_gamma_inc(const double a,const double x)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 }
668