rgamma.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h> // DBL_EPSILON, DBL_MAX_10_EXP
#include <time.h>
#include <string.h>
/*****************************************************************************/
#include "tpcextensions.h"
/*****************************************************************************/
#include "tpcfunc.h"
/*****************************************************************************/
 
/*****************************************************************************/
double igam(
  double a,
  double x
) {
  /* Check parameters */
  if(x==0.0) return(0.0);
  if((x<0.0) || (a<=0.0)) return(nan(""));
 
  if((x>1.0) && (x>a)) return(1.0 - igamc(a, x));
 
  /* Left tail of incomplete Gamma function:
    x^a * e^-x * Sum(k=0,Inf)(x^k / Gamma(a+k+1)) */
 
  double ans, ax, c, r;
 
  /* Compute  x**a * exp(-x) / Gamma(a) */
  ax=a*log(x) - x - lgamma(a);
  if(ax<-DBL_MAX_10_EXP) return(0.0); // underflow
  ax=exp(ax);
 
  /* power series */
  r=a; c=1.0; ans=1.0;
  do {
    r+=1.0;
    c*=x/r;
    ans+=c;
  } while(c/ans > DBL_EPSILON);
  return(ans*ax/a);
}
/*****************************************************************************/
 
/*****************************************************************************/
double igamc(
  double a,
  double x
) {
  /* Check parameters */
  if((x<0.0) || (a<=0.0)) return(nan(""));
 
  if((x<1.0) || (x<a)) return(1.0-igam(a, x));
 
  double ans, ax, c, yc, r, t, y, z;
  double pk, pkm1, pkm2, qk, qkm1, qkm2;
  double big=4.503599627370496E+015;
  double biginv=2.22044604925031308085E-016;
 
  ax = a*log(x) - x - lgamma(a);
  if(ax < -DBL_MAX_10_EXP) return(0.0); // underflow
  ax=exp(ax);
 
  /* continued fraction */
  y=1.0-a; z=x+y+1.0; c=0.0;
  pkm2=1.0; qkm2=x; pkm1=x+1.0; qkm1=z*x;
  ans=pkm1/qkm1;
  do {
    c+=1.0; y+=1.0; z+=2.0;
    yc=y*c; pk=pkm1*z - pkm2*yc; qk=qkm1*z - qkm2*yc;
    if(qk!=0.0) {r=pk/qk; t=fabs((ans-r)/r); ans=r;}
    else t=1.0;
    pkm2=pkm1; pkm1=pk; qkm2=qkm1; qkm1=qk;
    if(fabs(pk)>big) {pkm2*=biginv; pkm1*=biginv; qkm2*=biginv; qkm1*=biginv;}
  } while(t>DBL_EPSILON);
  return(ans*ax);
}
/*****************************************************************************/
 
/*****************************************************************************/