/* Logarithm of Gamma function, Gamma itself and Beta -- for the value 
   of argument(s) >0. Double precision.
   Copyright (c) Nikitin V.F. 2000

   Calculate logarithm of gamma-function for x>0
   double GammLn(double x);
       
   Calculate gamma-function for x>0
   double Gamma(double x);

   Calculate beta-function for x>0, y>0
   double Beta(double x,double y);

   Note: no domain check is implemented. 
*/

/* math.h included for log() and exp() */
#include <math.h>

/* the coefficients table */
#define CN 8
static double cof[CN]={
       2.5066282746310005,
       1.0000000000190015,
       76.18009172947146,
      -86.50532032941677,
       24.01409824083091,
      -1.231739572450155,
       0.1208650973866179e-2,
      -0.5395239384953e-5,
};

/* logarithm of gamma-function */
double GammLn(double x) {
  double y,ser,*co;
  int j;
  /* calculate the series */
  ser=cof[1]; y=x; co=cof+2;
  for(j=2;j<CN;j++) {
    y+=1.; ser+=(*co)/y; co++;
  }
  /* and the other parts of the function */
  y=x+5.5;
  y-=(x+0.5)*log(y);
  return(-y+log(cof[0]*ser/x));
}

/* the gamma-function itself */
double Gamma(double x) {
  return(exp(GammLn(x)));
}

/* the beta function */
double Beta(double x,double y) {
  return(exp(GammLn(x)+GammLn(y)-GammLn(x+y)));
}