gammai.c File Reference

Routines for computing the incomplete gamma functions. More...

#include "mdb.h"

Go to the source code of this file.

Macros
#define	GAMMAI_ACCURACY   1e-12
#define	MAX_SERIES   1000

Functions
double	gammaIncSeries (double a, double x)
double	gammaIncCFrac (double a, double x)
double	gammaP (double a, double x)
	Compute the regularized lower incomplete gamma function P(a,x).
double	gammaQ (double a, double x)
	Compute the regularized upper incomplete gamma function Q(a,x).

Detailed Description

Routines for computing the incomplete gamma functions.

Copyright

• (c) 2002 The University of Chicago, as Operator of Argonne National Laboratory.
• (c) 2002 The Regents of the University of California, as Operator of Los Alamos National Laboratory.

License

This file is distributed under the terms of the Software License Agreement found in the file LICENSE included with this distribution.

Author

M. Borland, C. Saunders, R. Soliday

Definition in file gammai.c.

Macro Definition Documentation

GAMMAI_ACCURACY

#define GAMMAI_ACCURACY   1e-12

Definition at line 27 of file gammai.c.

MAX_SERIES

#define MAX_SERIES   1000

Definition at line 78 of file gammai.c.

Function Documentation

gammaIncCFrac()

double gammaIncCFrac	(	double	a,
		double	x )

Definition at line 99 of file gammai.c.

                                         {
  double A0, B0, A1, B1, A_1, B_1;
  double an, bn, f1, f2, accuracy = 0, factor = 0;
  long m;
 
#if defined(vxWorks)
  fprintf(stderr, "lgamma function not implemented in vxWorks");
  exit(1);
#else
  accuracy = GAMMAI_ACCURACY / (factor = exp(-x - lgamma(a) + a * log(x)));
#endif
  A0 = 0;
  B0 = 1;
  A_1 = 1;
  B_1 = 0;
  an = 1;
  bn = x + 1 - a;
  A1 = bn * A0 + an * A_1;
  B1 = bn * B0 + an * B_1;
  f2 = A1 / B1;
  m = 1;
  do {
    A_1 = A0;
    B_1 = B0;
    A0 = A1;
    B0 = B1;
    f1 = f2;
    an = -m * (m - a);
    bn += 2;
    A1 = bn * A0 + an * A_1;
    B1 = bn * B0 + an * B_1;
    f2 = A1 / B1;
    if (B1) {
      /* rescale the recurrence to avoid over/underflow */
      A0 /= B1;
      B0 /= B1;
      A1 /= B1;
      B1 = 1;
    }
    m++;
  } while (m < MAX_SERIES && fabs(f1 - f2) > accuracy);
  return factor * f2;
}

gammaIncSeries()

double gammaIncSeries	(	double	a,
		double	x )

Definition at line 80 of file gammai.c.

                                          {
  double term = 0, sum;
  long n;
#if defined(vxWorks)
  fprintf(stderr, "lgamma function not implemented in vxWorks");
  exit(1);
#else
  term = exp(-x) * pow(x, a) / exp(lgamma(a + 1));
#endif
  sum = 0;
  n = 0;
  do {
    sum += term;
    n++;
    term *= x / (a + n);
  } while (term > GAMMAI_ACCURACY && n < MAX_SERIES);
  return (sum + term);
}

gammaP()

double gammaP	(	double	a,
		double	x )

Compute the regularized lower incomplete gamma function P(a,x).

For given parameters a > 0 and x ≥ 0, this function returns the value of P(a,x) = γ(a,x)/Γ(a), where γ(a,x) is the lower incomplete gamma function. Uses series expansion for x < a+1, and a continued fraction expansion otherwise.

Parameters

a	Shape parameter (a > 0).
x	The upper limit of the integral (x ≥ 0).

Returns

The value of P(a,x), or -1 if inputs are invalid.

Definition at line 40 of file gammai.c.

                                  {
  if (a <= 0 || x < 0)
    return -1;
  if (x == 0)
    return 0;
  if (x < a + 1) {
    /* use series */
    return gammaIncSeries(a, x);
  } else {
    /* use continued fraction */
    return 1 - gammaIncCFrac(a, x);
  }
}

gammaQ()

double gammaQ	(	double	a,
		double	x )

Compute the regularized upper incomplete gamma function Q(a,x).

For given parameters a > 0 and x ≥ 0, this function returns Q(a,x) = 1 - P(a,x). It uses a series expansion when x < a+1 and a continued fraction expansion otherwise.

Parameters

a	Shape parameter (a > 0).
x	The upper limit of the integral (x ≥ 0).

Returns

The value of Q(a,x), or -1 if inputs are invalid.

Definition at line 64 of file gammai.c.

                                  {
  if (a <= 0 || x < 0)
    return -1;
  if (x == 0)
    return 0;
  if (x < a + 1) {
    /* use series */
    return 1 - gammaIncSeries(a, x);
  } else {
    /* use continued fraction */
    return gammaIncCFrac(a, x);
  }
}