// incomplete-gamma.h
//
// Kevin Karplus
//	12 Sept 2003

#ifndef INCOMPLETE_GAMMA_H
#define INCOMPLETE_GAMMA_H

// Computation of the incomplete gamma function P(a,x)
//
// The algorithm is based on the formulas and code as denoted in
// Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).


// The incomplete gamma function is useful, because it is
// the cumulative distribution function for a gamma distribution.
// In particular, for the gamma pdf with shape=a and scale=sigma
// the cdf is GammaP(a, x/sigma).
// Note: mean = a*sigma, var= a*sigma^2
double GammaP(double a,double x);
double LogGammaP(double a,double x);

 
// Because we often want to compute P-values for large x,
// it is computationally more sensible to use the GammaQ=1-GammaP function,
// when GammaQ(a,x) will be close to one.
// We can avoid some folating-point underflow problems by using LogGammaQ,
// which is log(GammaQ).
double GammaQ(double a, double x);
double LogGammaQ(double a, double x);



// Note:  chi^2 distribution with degrees of freedom df is
//	gamma distribution with a=df/2, scale=2 
// 	Mean = df, Var = 2*df

inline double chi2_Q(double df, double x)
{    return GammaQ(0.5*df, 0.5*x);
}

// 14 April 2004 Kevin Karplus
//	Added LogGammaP and LogGammaQ procedures.
#endif