This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Re: gamma(a,x): forwarded message from Mark Galassi
- From: Gerard Jungman <jungman at lanl dot gov>
- To: gsl-discuss at sources dot redhat dot com
- Cc: Carlo Graziani <carlo at oddjob dot uchicago dot edu>
- Date: 01 Apr 2003 21:27:51 -0700
- Subject: Re: gamma(a,x): forwarded message from Mark Galassi
- Organization: Los Alamos National Laboratory
- References: <16009.50738.562839.444105@hppav.local>
Hi Carlo. Fresh in CVS: an implementation of the non-normalized
incomplete gamma function. Prototypes:
int gsl_sf_gamma_inc_e(double a, double x, gsl_sf_result & result);
double gsl_sf_gamma_inc(double a, double x);
It tests ok on non-crazy inputs, but let me know if you have
any problem with it.
Note to developers (or anybody who cares): that file, gamma_inc.c,
is becoming a mess...
--
Gerard Jungman <jungman at lanl dot gov>
Los Alamos National Laboratory
________________________________________________________________________
>
> From: Mark Galassi <rosalia at galassi dot org>
> To: Brian Gough <bjg at network-theory dot co dot uk>
> Subject: [Carlo Graziani <carlo at oddjob dot uchicago dot edu>] Incomplete Gamma function in GSL
> Date: 01 Apr 2003 09:10:03 -0700
>
> From: Carlo Graziani <carlo at oddjob dot uchicago dot edu>
> To: Mark Galassi <rosalia at lanl dot gov>
> Subject: Incomplete Gamma function in GSL
> Date: 31 Mar 2003 22:20:55 -0600
>
> Dude,
>
> Ran into a problem with gsl_sf_gamma_inc_Q, the incomplete gamma
> function. Actually, my problem is with missing functionality in the
> gamma function support.
>
> I'm doing integrals of Band spectra, the low-energy component of which,
> you may recall, has the form (E^alpha)*Exp(-E/E0). Typical alphas are
> around -1.0 - -1.5.
>
> If I had an *un-normalized* incomplete gamma function
>
> g(a,x) = \int_x^\infty dy y^(a-1) Exp(-y)
>
> then the integral over the low-energy part of the Band spectrum could be
> neatly expressed as the difference of two invocations of g. Note that
> this integral is defined for all real a, and for positive x.
>
> The gsl offers a normalized variant of this incomplete gamma function in
> gsl_sf_gamma_inc_Q:
>
> gsl_sf_gamma_inc_Q = (1/Gamma(a)) * \int_x^\infty dy y^(a-1) Exp(-y)
>
> The problem is, the normalization factor is defined only for a>0. In my
> case, that corresponds to alpha>-1. That's a disaster, because of the
> typical range of alpha.
>
> The thing is, that normalization factor is gratuitous, and removing it
> would increase the parameter range over which the function may be
> validly invoked. This is how the (fortran) SLATEC routine DGAMIC
> operates, using Gautschi's algorithm (see
> http://gams.nist.gov/serve.cgi/Module/SLATEC/DGAMIC/7066/).
>
> So, if you're taking feature requests for future releases, an
> unnormalized incomplete gamma function would be nice.
>
> Carlo