[Bug-gsl] GSL : lower incomplete gamma function for negative x, not really a bug...

[Matthias Schabel]

Dear Dr. Jungman,

I have recently been working on a problem for which evaluation of the  
lower incomplete gamma function, \gamma(\alpha,x), is needed for  
negative values of $x$.  As the current implementation of this  
function in GSL
is restricted to positive values, (with the proviso that I am most  
definitely not a numericist) I have made a quick implementation of  
the series expression from Abramowitz and Stegun (Eq. 6.5.4 and  
6.5.29) :

P(\alpha,x) = x^\alpha e^{-x} \sum_{n=0}^\infty \frac{x^n}{\Gamma(a+n 
+1)}

which appears to work well - I have compared the results with those  
given by Mathematica.  Perhaps it would be worthwhile to include this  
expression in the GSL implementation to extend its utility?

Best Regards,

Matthias

double gamma_star(double a,double x)
{
    // Use Abramowitz and Stegun equations 6.5.4 and 6.5.29
    const double    tol = 1e-10;
    const int            maxit = 100;

    // use Gamma(z+1) = z Gamma(z) to eliminate repeated calls to Gamma
    // n = 0 term
    int            n = 0;
    double    gstar = 1.0/Gamma(a+1);
    double    sum = gstar,
                    fac;

    do
    {
        fac = z/(a+n+1);
        gstar *= fac;
        sum += gstar;
        ++n;
    }
    while (fabs(gstar) >= tol & n <= maxit)

    return exp(-z)*sum;
}

double gamma(double a,double x)
{
    if (x >= 0)
    {
        return Gamma(a)*P(a,x);
    }
    else
    {
        return Gamma(a)*x^a*gamma_star(a,x);
    }

    // should never get here...
    return NaN;
}

------------------------------------------------------------------------ 
---------------------------
Matthias Schabel, Ph.D.
Assistant Professor, Department of Radiology
Utah Center for Advanced Imaging Research
729 Arapeen Drive
Salt Lake City, UT  84108
801-587-9413 (work)
801-585-3592 (fax)
801-706-5760 (cell)
801-484-0811 (home)
matthias dot schabel at hsc dot utah dot edu

Matthias Schabel.vcf
Description: Text Data