www.pudn.com > ghs1.2.rar > gammain.cpp

/*
 * Copyright (c) 2000-2003 Illinois Institute of Technology.
 *                         All rights reserved.
 *
 * This file is part of the GHS software package.  For license
 * information, see the LICENSE file in the top level directory of the
 * GHS source distribution.
 *
 * Released Date: August 18 2005
 * This is the GHS I release by the
 *   SCS Group,
 *   http://www.cs.iit.edu/~scs
 *   Department of Computer Science,
 *   Illinois Institute of Technology.
 *
 * This release has been tested under SUN OS 5.9.
 *
 * Supervisor
 * ----------
 * Illiniois Institute of Technology:
 * -Dr. Xian-He Sun
 *
 * Author(s)
 * -----------------
 * Illinois Institute of Technology:
 * -Ming Wu
 * -Xian-He Sun
 *
 */

#include
#include

using namespace std;

extern double Gammaln(double x);

double 
Gammainc(double x, double a)
{

  double gam, b;
  double ap, sum, del;
  double eps = 0.000001;
  double a0 = 0, a1 = 0, b0 = 0, b1 = 0, n = 0, ana = 0, g = 0, gold = 0, anf = 0, fac = 0;
  double realmin = 0.000000000001;

  gam = Gammaln(a + realmin);

  //gam = 12.8018;
	
  b = x;

  if (x == 0) b = 0;
  if (a == 0) b = 1;

  // Series expansion for x < a+1

  if ( (a != 0) && (x != 0) && (x < (a + 1)) )
  {
    ap = a;
    sum = 1 / ap;
    del = sum;
		
    while ( fabs(del) >= 10 * eps * sum )
    {
      ap = ap + 1;
      del = x * del / ap;
      sum = sum + del;
    }
		
    b = sum * exp(-1 * x + a * log(x) - gam);

  }


  //Continued fraction for x >= a+1
  if ( (a != 0) && (x != 0) && (x >= (a + 1)) )
  {
    a0 = 1;
    a1 = x;
    b0 = 0;
    b1 = a0;
    fac = 1;
    n = 1;
    g = b1;
    gold = b0;
   
    while ( fabs(g - gold) >= 10 * eps * g)
    {
      gold = g;
      ana = n - a;
      a0 = (a1 + a0 * ana) * fac;
      b0 = (b1 + b0 * ana) * fac;
      anf = n * fac;
      a1 = x * a0 + anf * a1;
      b1 = x * b0 + anf * b1;
      fac = 1 / a1;
      g = b1 * fac;
      n = n + 1;
    }
		
    b = 1 - exp(-1.0 * x + a * log(x) - gam) * g;

  }

  return(b);

}

double 
GammaCDF(double x, double a, double b)
{

  double res;

  if ( (x >= 0) && (a > 0) && (b > 0) )
  {
    res = Gammainc(x / b, a);
  } 
  else 
  {
    cout << "the Gammainc input parameters are out of range" << endl;
    cout << "x=" << x << " a=" << a << " b=" << b <