www.pudn.com > ghs1.2.rar > gammaln.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;

double 
Gammaln(double x)
{

  double d1 = -0.5772156649015328605195174;

  double p1[] = {4.945235359296727046734888, 201.8112620856775083915565, 
				2290.838373831346393026739, 11319.67205903380828685045, 
				28557.24635671635335736389, 38484.96228443793359990269, 
				26377.48787624195437963534, 7225.813979700288197698961};

  double q1[] = {67.48212550303777196073036, 1113.332393857199323513008, 
				7738.757056935398733233834, 27639.87074403340708898585, 
				54993.10206226157329794414, 61611.22180066002127833352, 
				36351.27591501940507276287, 8785.536302431013170870835};

  double d2 = 0.4227843350984671393993777;

  double p2[] = {4.974607845568932035012064, 542.4138599891070494101986, 
           15506.93864978364947665077, 184793.2904445632425417223, 
           1088204.769468828767498470, 3338152.967987029735917223, 
           5106661.678927352456275255, 3074109.054850539556250927};

  double q2[] = {183.0328399370592604055942, 7765.049321445005871323047, 
				133190.3827966074194402448, 1136705.821321969608938755, 
				5267964.117437946917577538, 13467014.54311101692290052, 
				17827365.30353274213975932, 9533095.591844353613395747};

  double d4 = 1.791759469228055000094023;

  double p4[] = {14745.02166059939948905062, 2426813.369486704502836312, 
				121475557.4045093227939592, 2663432449.630976949898078, 
				29403789566.34553899906876, 170266573776.5398868392998, 
				492612579337.7430887588120, 560625185622.3951465078242};

  double q4[] = {2690.530175870899333379843, 639388.5654300092398984238, 
				41355999.30241388052042842, 1120872109.616147941376570, 
				14886137286.78813811542398, 101680358627.2438228077304, 
				341747634550.7377132798597, 446315818741.9713286462081};

  double c[] = {-0.001910444077728, 0.00084171387781295, 
				-0.0005952379913043012, 0.000793650793500350248, 
				-0.002777777777777681622553, 0.008333333333333333331554247, 
				0.0057083835261};
	
  double res;
  double eps = 0.0000001, y, xden, xnum, xm1, xm2, xm4, r;
  double ysq, corr, spi;
  int i = 0;

  if (x < 0) 
  {
    cout << "Input arguments must be real and nonnegative " << endl;
    exit(0);
  }
	
  res = x;

  if (x < eps) 
  {
    res = -1 * log(x);
  }

  /*
   * eps < x <= 0.5
   */

  if (x > eps && x <= 0.5)
  {

    y = x;
    xden = 1;
    xnum = 0;
		
    for (i = 0; i <= 7; i++) 
    {
      xnum = xnum * y + p1[i];
      xden = xden * y + q1[i];
    }
      
    //printf(" y= %f, d1=%f, xnum=%f, xden=%f \n", y, d1, xnum, xden);
    res = -1.0 * log(y) + (y * (d1 + y * (xnum / xden)));

    //printf("the log of gamma is %f \n", res);
  }

  /*
   *0.5 < x <= 0.6796875
   */
 
  if (x > 0.5 && x <= 0.6796875)
  {
    xm1 = (x - 0.5) - 0.5;
    xden = 1;
    xnum = 0;

    for (i = 0; i <= 7;i++)
    {
      xnum = xnum * xm1 + p2[i];
      xden = xden * xm1 + q2[i];
    }

    res = -log(x) + xm1 * (d2 + xm1 * (xnum / xden));

    //cout << "the log of gamma is " << res << endl;
  }

  /*
   * 0.6796875 < x <= 1.5
   */

  if (x > 0.6796875 && x <= 1.5)
  {
    xm1 = (x - 0.5) - 0.5;
    xden = 1;
    xnum = 0;

    for (i = 0; i <= 7; i++)
    {
      xnum = xnum * xm1 + p1[i];
      xden = xden * xm1 + q1[i];
    }
		
    res = xm1 * (d1 + xm1 * (xnum / xden));

    //cout << "the log of gamma is " << res << endl;
  }

  /*
   * 1.5 < x <= 4
   */

  if (x > 1.5 && x <= 4)
  {
    xm2 = x - 2;
    xden = 1;
    xnum = 0;

    for (i = 0; i <= 7; i++)
    {
      xnum = xnum * xm2 + p2[i];
      xden = xden * xm2 + q2[i];
    }

    res = xm2 * (d2 + xm2 * (xnum / xden));

    //cout << "the log of gamma is " << res << endl;
  }

  /*
   * 4 < x <= 12
   */

  if (x > 4 && x <= 12)
  {

    xm4 = x - 4;
    xden = -1;
    xnum = 0;

    for (i = 0; i <= 7; i++)
    {
      xnum = xnum * xm4 + p4[i];
      xden = xden * xm4 + q4[i];
    }

    res = d4 + xm4 * (xnum / xden);

    //cout << "the log of gamma is " << res << endl;
  }

  /*
   * x > 12
   */

  if (x > 12)
  {
    y = x;
    r = c[6] * 1;
    ysq = y * y;
      
    for (i = 0; i <= 5; i++)
    {
      r = r / ysq + c[i];
    }

    r = r / y;
    corr = log(y);
    spi = 0.9189385332046727417803297;
    res = r + spi - 0.5 * corr + y * (corr - 1);

    //cout << "the log of gamma is " << res << endl;
  }

  return(res);

}