|
Packit |
67cb25 |
/* specfunc/gamma_inc.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2007 Brian Gough
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Author: G. Jungman */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_erf.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_exp.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_log.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_gamma.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_expint.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "error.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* The dominant part,
|
|
Packit |
67cb25 |
* D(a,x) := x^a e^(-x) / Gamma(a+1)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gamma_inc_D(const double a, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(a < 10.0) {
|
|
Packit |
67cb25 |
double lnr;
|
|
Packit |
67cb25 |
gsl_sf_result lg;
|
|
Packit |
67cb25 |
gsl_sf_lngamma_e(a+1.0, &lg);
|
|
Packit |
67cb25 |
lnr = a * log(x) - x - lg.val;
|
|
Packit |
67cb25 |
result->val = exp(lnr);
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * (fabs(lnr) + 1.0) * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
gsl_sf_result gstar;
|
|
Packit |
67cb25 |
gsl_sf_result ln_term;
|
|
Packit |
67cb25 |
double term1;
|
|
Packit |
67cb25 |
if (x < 0.5*a) {
|
|
Packit |
67cb25 |
double u = x/a;
|
|
Packit |
67cb25 |
double ln_u = log(u);
|
|
Packit |
67cb25 |
ln_term.val = ln_u - u + 1.0;
|
|
Packit |
67cb25 |
ln_term.err = (fabs(ln_u) + fabs(u) + 1.0) * GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
} else {
|
|
Packit |
67cb25 |
double mu = (x-a)/a;
|
|
Packit |
67cb25 |
gsl_sf_log_1plusx_mx_e(mu, &ln_term); /* log(1+mu) - mu */
|
|
Packit |
67cb25 |
/* Propagate cancellation error from x-a, since the absolute
|
|
Packit |
67cb25 |
error of mu=x-a is DBL_EPSILON */
|
|
Packit |
67cb25 |
ln_term.err += GSL_DBL_EPSILON * fabs(mu);
|
|
Packit |
67cb25 |
};
|
|
Packit |
67cb25 |
gsl_sf_gammastar_e(a, &gstar);
|
|
Packit |
67cb25 |
term1 = exp(a*ln_term.val)/sqrt(2.0*M_PI*a);
|
|
Packit |
67cb25 |
result->val = term1/gstar.val;
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * (fabs(a*ln_term.val) + 1.0) * fabs(result->val);
|
|
Packit |
67cb25 |
/* Include propagated error from log term */
|
|
Packit |
67cb25 |
result->err += fabs(a) * ln_term.err * fabs(result->val);
|
|
Packit |
67cb25 |
result->err += gstar.err/fabs(gstar.val) * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* P series representation.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gamma_inc_P_series(const double a, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const int nmax = 10000;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_sf_result D;
|
|
Packit |
67cb25 |
int stat_D = gamma_inc_D(a, x, &D);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Approximating the terms of the series using Stirling's
|
|
Packit |
67cb25 |
approximation gives t_n = (x/a)^n * exp(-n(n+1)/(2a)), so the
|
|
Packit |
67cb25 |
convergence condition is n^2 / (2a) + (1-(x/a) + (1/2a)) n >>
|
|
Packit |
67cb25 |
-log(GSL_DBL_EPS) if we want t_n < O(1e-16) t_0. The condition
|
|
Packit |
67cb25 |
below detects cases where the minimum value of n is > 5000 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (x > 0.995 * a && a > 1e5) { /* Difficult case: try continued fraction */
|
|
Packit |
67cb25 |
gsl_sf_result cf_res;
|
|
Packit |
67cb25 |
int status = gsl_sf_exprel_n_CF_e(a, x, &cf_res);
|
|
Packit |
67cb25 |
result->val = D.val * cf_res.val;
|
|
Packit |
67cb25 |
result->err = fabs(D.val * cf_res.err) + fabs(D.err * cf_res.val);
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Series would require excessive number of terms */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (x > (a + nmax)) {
|
|
Packit |
67cb25 |
GSL_ERROR ("gamma_inc_P_series x>>a exceeds range", GSL_EMAXITER);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Normal case: sum the series */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double sum = 1.0;
|
|
Packit |
67cb25 |
double term = 1.0;
|
|
Packit |
67cb25 |
double remainder;
|
|
Packit |
67cb25 |
int n;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Handle lower part of the series where t_n is increasing, |x| > a+n */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int nlow = (x > a) ? (x - a): 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(n=1; n < nlow; n++) {
|
|
Packit |
67cb25 |
term *= x/(a+n);
|
|
Packit |
67cb25 |
sum += term;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Handle upper part of the series where t_n is decreasing, |x| < a+n */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (/* n = previous n */ ; n
|
|
Packit |
67cb25 |
term *= x/(a+n);
|
|
Packit |
67cb25 |
sum += term;
|
|
Packit |
67cb25 |
if(fabs(term/sum) < GSL_DBL_EPSILON) break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Estimate remainder of series ~ t_(n+1)/(1-x/(a+n+1)) */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double tnp1 = (x/(a+n)) * term;
|
|
Packit |
67cb25 |
remainder = tnp1 / (1.0 - x/(a + n + 1.0));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = D.val * sum;
|
|
Packit |
67cb25 |
result->err = D.err * fabs(sum) + fabs(D.val * remainder);
|
|
Packit |
67cb25 |
result->err += (1.0 + n) * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(n == nmax && fabs(remainder/sum) > GSL_SQRT_DBL_EPSILON)
|
|
Packit |
67cb25 |
GSL_ERROR ("gamma_inc_P_series failed to converge", GSL_EMAXITER);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return stat_D;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Q large x asymptotic
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gamma_inc_Q_large_x(const double a, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const int nmax = 5000;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_sf_result D;
|
|
Packit |
67cb25 |
const int stat_D = gamma_inc_D(a, x, &D);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double sum = 1.0;
|
|
Packit |
67cb25 |
double term = 1.0;
|
|
Packit |
67cb25 |
double last = 1.0;
|
|
Packit |
67cb25 |
int n;
|
|
Packit |
67cb25 |
for(n=1; n
|
|
Packit |
67cb25 |
term *= (a-n)/x;
|
|
Packit |
67cb25 |
if(fabs(term/last) > 1.0) break;
|
|
Packit |
67cb25 |
if(fabs(term/sum) < GSL_DBL_EPSILON) break;
|
|
Packit |
67cb25 |
sum += term;
|
|
Packit |
67cb25 |
last = term;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = D.val * (a/x) * sum;
|
|
Packit |
67cb25 |
result->err = D.err * fabs((a/x) * sum);
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(n == nmax)
|
|
Packit |
67cb25 |
GSL_ERROR ("error in large x asymptotic", GSL_EMAXITER);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return stat_D;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Uniform asymptotic for x near a, a and x large.
|
|
Packit |
67cb25 |
* See [Temme, p. 285]
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gamma_inc_Q_asymp_unif(const double a, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double rta = sqrt(a);
|
|
Packit |
67cb25 |
const double eps = (x-a)/a;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_sf_result ln_term;
|
|
Packit |
67cb25 |
const int stat_ln = gsl_sf_log_1plusx_mx_e(eps, &ln_term); /* log(1+eps) - eps */
|
|
Packit |
67cb25 |
const double eta = GSL_SIGN(eps) * sqrt(-2.0*ln_term.val);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_sf_result erfc;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double R;
|
|
Packit |
67cb25 |
double c0, c1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* This used to say erfc(eta*M_SQRT2*rta), which is wrong.
|
|
Packit |
67cb25 |
* The sqrt(2) is in the denominator. Oops.
|
|
Packit |
67cb25 |
* Fixed: [GJ] Mon Nov 15 13:25:32 MST 2004
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
gsl_sf_erfc_e(eta*rta/M_SQRT2, &erfc);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(fabs(eps) < GSL_ROOT5_DBL_EPSILON) {
|
|
Packit |
67cb25 |
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)));
|
|
Packit |
67cb25 |
c1 = -1.0/540.0 - eps/288.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
const double rt_term = sqrt(-2.0 * ln_term.val/(eps*eps));
|
|
Packit |
67cb25 |
const double lam = x/a;
|
|
Packit |
67cb25 |
c0 = (1.0 - 1.0/rt_term)/eps;
|
|
Packit |
67cb25 |
c1 = -(eta*eta*eta * (lam*lam + 10.0*lam + 1.0) - 12.0 * eps*eps*eps) / (12.0 * eta*eta*eta*eps*eps*eps);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
R = exp(-0.5*a*eta*eta)/(M_SQRT2*M_SQRTPI*rta) * (c0 + c1/a);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = 0.5 * erfc.val + R;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * fabs(R * 0.5 * a*eta*eta) + 0.5 * erfc.err;
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return stat_ln;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Continued fraction which occurs in evaluation
|
|
Packit |
67cb25 |
* of Q(a,x) or Gamma(a,x).
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* 1 (1-a)/x 1/x (2-a)/x 2/x (3-a)/x
|
|
Packit |
67cb25 |
* F(a,x) = ---- ------- ----- -------- ----- -------- ...
|
|
Packit |
67cb25 |
* 1 + 1 + 1 + 1 + 1 + 1 +
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Hans E. Plesser, 2002-01-22 (hans dot plesser at itf dot nlh dot no).
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Split out from gamma_inc_Q_CF() by GJ [Tue Apr 1 13:16:41 MST 2003].
|
|
Packit |
67cb25 |
* See gamma_inc_Q_CF() below.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
gamma_inc_F_CF(const double a, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const int nmax = 5000;
|
|
Packit |
67cb25 |
const double small = gsl_pow_3 (GSL_DBL_EPSILON);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double hn = 1.0; /* convergent */
|
|
Packit |
67cb25 |
double Cn = 1.0 / small;
|
|
Packit |
67cb25 |
double Dn = 1.0;
|
|
Packit |
67cb25 |
int n;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* n == 1 has a_1, b_1, b_0 independent of a,x,
|
|
Packit |
67cb25 |
so that has been done by hand */
|
|
Packit |
67cb25 |
for ( n = 2 ; n < nmax ; n++ )
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double an;
|
|
Packit |
67cb25 |
double delta;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(GSL_IS_ODD(n))
|
|
Packit |
67cb25 |
an = 0.5*(n-1)/x;
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
an = (0.5*n-a)/x;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Dn = 1.0 + an * Dn;
|
|
Packit |
67cb25 |
if ( fabs(Dn) < small )
|
|
Packit |
67cb25 |
Dn = small;
|
|
Packit |
67cb25 |
Cn = 1.0 + an/Cn;
|
|
Packit |
67cb25 |
if ( fabs(Cn) < small )
|
|
Packit |
67cb25 |
Cn = small;
|
|
Packit |
67cb25 |
Dn = 1.0 / Dn;
|
|
Packit |
67cb25 |
delta = Cn * Dn;
|
|
Packit |
67cb25 |
hn *= delta;
|
|
Packit |
67cb25 |
if(fabs(delta-1.0) < GSL_DBL_EPSILON) break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = hn;
|
|
Packit |
67cb25 |
result->err = 2.0*GSL_DBL_EPSILON * fabs(hn);
|
|
Packit |
67cb25 |
result->err += GSL_DBL_EPSILON * (2.0 + 0.5*n) * fabs(result->val);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(n == nmax)
|
|
Packit |
67cb25 |
GSL_ERROR ("error in CF for F(a,x)", GSL_EMAXITER);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Continued fraction for Q.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Q(a,x) = D(a,x) a/x F(a,x)
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Hans E. Plesser, 2002-01-22 (hans dot plesser at itf dot nlh dot no):
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Since the Gautschi equivalent series method for CF evaluation may lead
|
|
Packit |
67cb25 |
* to singularities, I have replaced it with the modified Lentz algorithm
|
|
Packit |
67cb25 |
* given in
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* I J Thompson and A R Barnett
|
|
Packit |
67cb25 |
* Coulomb and Bessel Functions of Complex Arguments and Order
|
|
Packit |
67cb25 |
* J Computational Physics 64:490-509 (1986)
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* In consequence, gamma_inc_Q_CF_protected() is now obsolete and has been
|
|
Packit |
67cb25 |
* removed.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Identification of terms between the above equation for F(a, x) and
|
|
Packit |
67cb25 |
* the first equation in the appendix of Thompson&Barnett is as follows:
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* b_0 = 0, b_n = 1 for all n > 0
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* a_1 = 1
|
|
Packit |
67cb25 |
* a_n = (n/2-a)/x for n even
|
|
Packit |
67cb25 |
* a_n = (n-1)/(2x) for n odd
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gamma_inc_Q_CF(const double a, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_sf_result D;
|
|
Packit |
67cb25 |
gsl_sf_result F;
|
|
Packit |
67cb25 |
const int stat_D = gamma_inc_D(a, x, &D);
|
|
Packit |
67cb25 |
const int stat_F = gamma_inc_F_CF(a, x, &F);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = D.val * (a/x) * F.val;
|
|
Packit |
67cb25 |
result->err = D.err * fabs((a/x) * F.val) + fabs(D.val * a/x * F.err);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_F, stat_D);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Useful for small a and x. Handles the subtraction analytically.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gamma_inc_Q_series(const double a, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double term1; /* 1 - x^a/Gamma(a+1) */
|
|
Packit |
67cb25 |
double sum; /* 1 + (a+1)/(a+2)(-x)/2! + (a+1)/(a+3)(-x)^2/3! + ... */
|
|
Packit |
67cb25 |
int stat_sum;
|
|
Packit |
67cb25 |
double term2; /* a temporary variable used at the end */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Evaluate series for 1 - x^a/Gamma(a+1), small a
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
const double pg21 = -2.404113806319188570799476; /* PolyGamma[2,1] */
|
|
Packit |
67cb25 |
const double lnx = log(x);
|
|
Packit |
67cb25 |
const double el = M_EULER+lnx;
|
|
Packit |
67cb25 |
const double c1 = -el;
|
|
Packit |
67cb25 |
const double c2 = M_PI*M_PI/12.0 - 0.5*el*el;
|
|
Packit |
67cb25 |
const double c3 = el*(M_PI*M_PI/12.0 - el*el/6.0) + pg21/6.0;
|
|
Packit |
67cb25 |
const double c4 = -0.04166666666666666667
|
|
Packit |
67cb25 |
* (-1.758243446661483480 + lnx)
|
|
Packit |
67cb25 |
* (-0.764428657272716373 + lnx)
|
|
Packit |
67cb25 |
* ( 0.723980571623507657 + lnx)
|
|
Packit |
67cb25 |
* ( 4.107554191916823640 + lnx);
|
|
Packit |
67cb25 |
const double c5 = -0.0083333333333333333
|
|
Packit |
67cb25 |
* (-2.06563396085715900 + lnx)
|
|
Packit |
67cb25 |
* (-1.28459889470864700 + lnx)
|
|
Packit |
67cb25 |
* (-0.27583535756454143 + lnx)
|
|
Packit |
67cb25 |
* ( 1.33677371336239618 + lnx)
|
|
Packit |
67cb25 |
* ( 5.17537282427561550 + lnx);
|
|
Packit |
67cb25 |
const double c6 = -0.0013888888888888889
|
|
Packit |
67cb25 |
* (-2.30814336454783200 + lnx)
|
|
Packit |
67cb25 |
* (-1.65846557706987300 + lnx)
|
|
Packit |
67cb25 |
* (-0.88768082560020400 + lnx)
|
|
Packit |
67cb25 |
* ( 0.17043847751371778 + lnx)
|
|
Packit |
67cb25 |
* ( 1.92135970115863890 + lnx)
|
|
Packit |
67cb25 |
* ( 6.22578557795474900 + lnx);
|
|
Packit |
67cb25 |
const double c7 = -0.00019841269841269841
|
|
Packit |
67cb25 |
* (-2.5078657901291800 + lnx)
|
|
Packit |
67cb25 |
* (-1.9478900888958200 + lnx)
|
|
Packit |
67cb25 |
* (-1.3194837322612730 + lnx)
|
|
Packit |
67cb25 |
* (-0.5281322700249279 + lnx)
|
|
Packit |
67cb25 |
* ( 0.5913834939078759 + lnx)
|
|
Packit |
67cb25 |
* ( 2.4876819633378140 + lnx)
|
|
Packit |
67cb25 |
* ( 7.2648160783762400 + lnx);
|
|
Packit |
67cb25 |
const double c8 = -0.00002480158730158730
|
|
Packit |
67cb25 |
* (-2.677341544966400 + lnx)
|
|
Packit |
67cb25 |
* (-2.182810448271700 + lnx)
|
|
Packit |
67cb25 |
* (-1.649350342277400 + lnx)
|
|
Packit |
67cb25 |
* (-1.014099048290790 + lnx)
|
|
Packit |
67cb25 |
* (-0.191366955370652 + lnx)
|
|
Packit |
67cb25 |
* ( 0.995403817918724 + lnx)
|
|
Packit |
67cb25 |
* ( 3.041323283529310 + lnx)
|
|
Packit |
67cb25 |
* ( 8.295966556941250 + lnx);
|
|
Packit |
67cb25 |
const double c9 = -2.75573192239859e-6
|
|
Packit |
67cb25 |
* (-2.8243487670469080 + lnx)
|
|
Packit |
67cb25 |
* (-2.3798494322701120 + lnx)
|
|
Packit |
67cb25 |
* (-1.9143674728689960 + lnx)
|
|
Packit |
67cb25 |
* (-1.3814529102920370 + lnx)
|
|
Packit |
67cb25 |
* (-0.7294312810261694 + lnx)
|
|
Packit |
67cb25 |
* ( 0.1299079285269565 + lnx)
|
|
Packit |
67cb25 |
* ( 1.3873333251885240 + lnx)
|
|
Packit |
67cb25 |
* ( 3.5857258865210760 + lnx)
|
|
Packit |
67cb25 |
* ( 9.3214237073814600 + lnx);
|
|
Packit |
67cb25 |
const double c10 = -2.75573192239859e-7
|
|
Packit |
67cb25 |
* (-2.9540329644556910 + lnx)
|
|
Packit |
67cb25 |
* (-2.5491366926991850 + lnx)
|
|
Packit |
67cb25 |
* (-2.1348279229279880 + lnx)
|
|
Packit |
67cb25 |
* (-1.6741881076349450 + lnx)
|
|
Packit |
67cb25 |
* (-1.1325949616098420 + lnx)
|
|
Packit |
67cb25 |
* (-0.4590034650618494 + lnx)
|
|
Packit |
67cb25 |
* ( 0.4399352987435699 + lnx)
|
|
Packit |
67cb25 |
* ( 1.7702236517651670 + lnx)
|
|
Packit |
67cb25 |
* ( 4.1231539047474080 + lnx)
|
|
Packit |
67cb25 |
* ( 10.342627908148680 + lnx);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
term1 = a*(c1+a*(c2+a*(c3+a*(c4+a*(c5+a*(c6+a*(c7+a*(c8+a*(c9+a*c10)))))))));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Evaluate the sum.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
const int nmax = 5000;
|
|
Packit |
67cb25 |
double t = 1.0;
|
|
Packit |
67cb25 |
int n;
|
|
Packit |
67cb25 |
sum = 1.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(n=1; n
|
|
Packit |
67cb25 |
t *= -x/(n+1.0);
|
|
Packit |
67cb25 |
sum += (a+1.0)/(a+n+1.0)*t;
|
|
Packit |
67cb25 |
if(fabs(t/sum) < GSL_DBL_EPSILON) break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(n == nmax)
|
|
Packit |
67cb25 |
stat_sum = GSL_EMAXITER;
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
stat_sum = GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
term2 = (1.0 - term1) * a/(a+1.0) * x * sum;
|
|
Packit |
67cb25 |
result->val = term1 + term2;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * (fabs(term1) + 2.0*fabs(term2));
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return stat_sum;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* series for small a and x, but not defined for a == 0 */
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
gamma_inc_series(double a, double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_sf_result Q;
|
|
Packit |
67cb25 |
gsl_sf_result G;
|
|
Packit |
67cb25 |
const int stat_Q = gamma_inc_Q_series(a, x, &Q);
|
|
Packit |
67cb25 |
const int stat_G = gsl_sf_gamma_e(a, &G);
|
|
Packit |
67cb25 |
result->val = Q.val * G.val;
|
|
Packit |
67cb25 |
result->err = fabs(Q.val * G.err) + fabs(Q.err * G.val);
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_Q, stat_G);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
gamma_inc_a_gt_0(double a, double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* x > 0 and a > 0; use result for Q */
|
|
Packit |
67cb25 |
gsl_sf_result Q;
|
|
Packit |
67cb25 |
gsl_sf_result G;
|
|
Packit |
67cb25 |
const int stat_Q = gsl_sf_gamma_inc_Q_e(a, x, &Q);
|
|
Packit |
67cb25 |
const int stat_G = gsl_sf_gamma_e(a, &G);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = G.val * Q.val;
|
|
Packit |
67cb25 |
result->err = fabs(G.val * Q.err) + fabs(G.err * Q.val);
|
|
Packit |
67cb25 |
result->err += 2.0*GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_G, stat_Q);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static int
|
|
Packit |
67cb25 |
gamma_inc_CF(double a, double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_sf_result F;
|
|
Packit |
67cb25 |
gsl_sf_result pre;
|
|
Packit |
67cb25 |
const double am1lgx = (a-1.0)*log(x);
|
|
Packit |
67cb25 |
const int stat_F = gamma_inc_F_CF(a, x, &F);
|
|
Packit |
67cb25 |
const int stat_E = gsl_sf_exp_err_e(am1lgx - x, GSL_DBL_EPSILON*fabs(am1lgx), &pre);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = F.val * pre.val;
|
|
Packit |
67cb25 |
result->err = fabs(F.err * pre.val) + fabs(F.val * pre.err);
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_F, stat_E);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* evaluate Gamma(0,x), x > 0 */
|
|
Packit |
67cb25 |
#define GAMMA_INC_A_0(x, result) gsl_sf_expint_E1_e(x, result)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_gamma_inc_Q_e(const double a, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(a < 0.0 || x < 0.0) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x == 0.0) {
|
|
Packit |
67cb25 |
result->val = 1.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a == 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
result->val = 0.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x <= 0.5*a) {
|
|
Packit |
67cb25 |
/* If the series is quick, do that. It is
|
|
Packit |
67cb25 |
* robust and simple.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
gsl_sf_result P;
|
|
Packit |
67cb25 |
int stat_P = gamma_inc_P_series(a, x, &P);
|
|
Packit |
67cb25 |
result->val = 1.0 - P.val;
|
|
Packit |
67cb25 |
result->err = P.err;
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return stat_P;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a >= 1.0e+06 && (x-a)*(x-a) < a) {
|
|
Packit |
67cb25 |
/* Then try the difficult asymptotic regime.
|
|
Packit |
67cb25 |
* This is the only way to do this region.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
return gamma_inc_Q_asymp_unif(a, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a < 0.2 && x < 5.0) {
|
|
Packit |
67cb25 |
/* Cancellations at small a must be handled
|
|
Packit |
67cb25 |
* analytically; x should not be too big
|
|
Packit |
67cb25 |
* either since the series terms grow
|
|
Packit |
67cb25 |
* with x and log(x).
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
return gamma_inc_Q_series(a, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a <= x) {
|
|
Packit |
67cb25 |
if(x <= 1.0e+06) {
|
|
Packit |
67cb25 |
/* Continued fraction is excellent for x >~ a.
|
|
Packit |
67cb25 |
* We do not let x be too large when x > a since
|
|
Packit |
67cb25 |
* it is somewhat pointless to try this there;
|
|
Packit |
67cb25 |
* the function is rapidly decreasing for
|
|
Packit |
67cb25 |
* x large and x > a, and it will just
|
|
Packit |
67cb25 |
* underflow in that region anyway. We
|
|
Packit |
67cb25 |
* catch that case in the standard
|
|
Packit |
67cb25 |
* large-x method.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
return gamma_inc_Q_CF(a, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
return gamma_inc_Q_large_x(a, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
if(x > a - sqrt(a)) {
|
|
Packit |
67cb25 |
/* Continued fraction again. The convergence
|
|
Packit |
67cb25 |
* is a little slower here, but that is fine.
|
|
Packit |
67cb25 |
* We have to trade that off against the slow
|
|
Packit |
67cb25 |
* convergence of the series, which is the
|
|
Packit |
67cb25 |
* only other option.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
return gamma_inc_Q_CF(a, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
gsl_sf_result P;
|
|
Packit |
67cb25 |
int stat_P = gamma_inc_P_series(a, x, &P);
|
|
Packit |
67cb25 |
result->val = 1.0 - P.val;
|
|
Packit |
67cb25 |
result->err = P.err;
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return stat_P;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_gamma_inc_P_e(const double a, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(a <= 0.0 || x < 0.0) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x == 0.0) {
|
|
Packit |
67cb25 |
result->val = 0.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < 20.0 || x < 0.5*a) {
|
|
Packit |
67cb25 |
/* Do the easy series cases. Robust and quick.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
return gamma_inc_P_series(a, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a > 1.0e+06 && (x-a)*(x-a) < a) {
|
|
Packit |
67cb25 |
/* Crossover region. Note that Q and P are
|
|
Packit |
67cb25 |
* roughly the same order of magnitude here,
|
|
Packit |
67cb25 |
* so the subtraction is stable.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
gsl_sf_result Q;
|
|
Packit |
67cb25 |
int stat_Q = gamma_inc_Q_asymp_unif(a, x, &Q);
|
|
Packit |
67cb25 |
result->val = 1.0 - Q.val;
|
|
Packit |
67cb25 |
result->err = Q.err;
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return stat_Q;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a <= x) {
|
|
Packit |
67cb25 |
/* Q <~ P in this area, so the
|
|
Packit |
67cb25 |
* subtractions are stable.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
gsl_sf_result Q;
|
|
Packit |
67cb25 |
int stat_Q;
|
|
Packit |
67cb25 |
if(a > 0.2*x) {
|
|
Packit |
67cb25 |
stat_Q = gamma_inc_Q_CF(a, x, &Q);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
stat_Q = gamma_inc_Q_large_x(a, x, &Q);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
result->val = 1.0 - Q.val;
|
|
Packit |
67cb25 |
result->err = Q.err;
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return stat_Q;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
if((x-a)*(x-a) < a) {
|
|
Packit |
67cb25 |
/* This condition is meant to insure
|
|
Packit |
67cb25 |
* that Q is not very close to 1,
|
|
Packit |
67cb25 |
* so the subtraction is stable.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
gsl_sf_result Q;
|
|
Packit |
67cb25 |
int stat_Q = gamma_inc_Q_CF(a, x, &Q);
|
|
Packit |
67cb25 |
result->val = 1.0 - Q.val;
|
|
Packit |
67cb25 |
result->err = Q.err;
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return stat_Q;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
return gamma_inc_P_series(a, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_gamma_inc_e(const double a, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(x < 0.0) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x == 0.0) {
|
|
Packit |
67cb25 |
return gsl_sf_gamma_e(a, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a == 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return GAMMA_INC_A_0(x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a > 0.0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return gamma_inc_a_gt_0(a, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x > 0.25)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* continued fraction seems to fail for x too small; otherwise
|
|
Packit |
67cb25 |
it is ok, independent of the value of |x/a|, because of the
|
|
Packit |
67cb25 |
non-oscillation in the expansion, i.e. the CF is
|
|
Packit |
67cb25 |
un-conditionally convergent for a < 0 and x > 0
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
return gamma_inc_CF(a, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(fabs(a) < 0.5)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return gamma_inc_series(a, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* a = fa + da; da >= 0 */
|
|
Packit |
67cb25 |
const double fa = floor(a);
|
|
Packit |
67cb25 |
const double da = a - fa;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_sf_result g_da;
|
|
Packit |
67cb25 |
const int stat_g_da = ( da > 0.0 ? gamma_inc_a_gt_0(da, x, &g_da)
|
|
Packit |
67cb25 |
: GAMMA_INC_A_0(x, &g_da));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double alpha = da;
|
|
Packit |
67cb25 |
double gax = g_da.val;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Gamma(alpha-1,x) = 1/(alpha-1) (Gamma(a,x) - x^(alpha-1) e^-x) */
|
|
Packit |
67cb25 |
do
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double shift = exp(-x + (alpha-1.0)*log(x));
|
|
Packit |
67cb25 |
gax = (gax - shift) / (alpha - 1.0);
|
|
Packit |
67cb25 |
alpha -= 1.0;
|
|
Packit |
67cb25 |
} while(alpha > a);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = gax;
|
|
Packit |
67cb25 |
result->err = 2.0*(1.0 + fabs(a))*GSL_DBL_EPSILON*fabs(gax);
|
|
Packit |
67cb25 |
return stat_g_da;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "eval.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_gamma_inc_P(const double a, const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_gamma_inc_P_e(a, x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_gamma_inc_Q(const double a, const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_gamma_inc_Q_e(a, x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_gamma_inc(const double a, const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_gamma_inc_e(a, x, &result));
|
|
Packit |
67cb25 |
}
|