[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] GSL : lower incomplete gamma function for negative x, not real
From: |
Matthias Schabel |
Subject: |
[Bug-gsl] GSL : lower incomplete gamma function for negative x, not really a bug... |
Date: |
Mon, 8 Aug 2005 23:04:17 -0600 |
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
- [Bug-gsl] GSL : lower incomplete gamma function for negative x, not really a bug...,
Matthias Schabel <=