rgamma.c File Reference

Regularized gamma function. More...

#include "tpcclibConfig.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <time.h>
#include <string.h>
#include "tpcextensions.h"
#include "tpcfunc.h"

Go to the source code of this file.

Functions
double	igam (double a, double x)
double	igamc (double a, double x)

Detailed Description

Regularized gamma function.

Note

C99 is required because of Gamma function, lgamma().

These functions could be used to fit the plasma parent fractions (Naganawa et al 2014).

Definition in file rgamma.c.

Function Documentation

igam()

double igam	(	double	a,
		double	x )

Cumulative gamma distribution, or Regularized gamma function, more specifically, lower incomplete gamma function divided by gamma function.

Standard gamma distribution is assumed (Beta=1). f(a,x) = (1/Gamma(a)) * Integral(0,x)(e^-t * t^(a-1))dt

Returns

Returns the value of regularized gamma function, or NaN in case of an error.

See also

igamc, inverfc, factorial

Parameters

a	Shape parameter alpha; must be > 0.
x	Integral from 0 to x; must be >= 0.

Definition at line 29 of file rgamma.c.

  {
  /* 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);
}

Referenced by igamc(), and mfEvalY().

igamc()

double igamc	(	double	a,
		double	x )

Regularized gamma function, more specifically, upper incomplete gamma function divided by gamma function.

f(a,x) = (1/Gamma(a)) * Integral(x,Inf)(e^-t * t^(a-1))dt Standard gamma distribution is assumed (Beta=1).

Returns

Returns the value of regularized gamma function, or NaN in case of an error.

See also

igam, inverfc, factorial

Parameters

a	Shape parameter alpha; must be > 0.
x	Integral from x to infinity; must be >= 0

Definition at line 70 of file rgamma.c.

  {
  /* 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);
}

Referenced by igam().