Re: gamma(a,x): forwarded message from Mark Galassi

[Gerard Jungman <jungman at lanl dot gov>]

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