Gamma function: implementation?

[mccaugh at s.cs.uiuc.edu]

 Brian:

I am posting this in case mailing fails (also in case anyone wants to correct):


#include <math.h>

/* Gamma  -  xx = argument,                                           */
/*           ier = return-code (0: no error; 1: negative; 2: overflow */

 double  Gamma (xx,ier)
         double xx;
            int *ier;
 {
  double  x,y, Gx,Gy, err;

#define  MAX_FLOAT  1.0E75

  if(xx > 57.0) { *ier = 2;  /* overflow */
		   return(MAX_FLOAT);
	        }
  /* else xx <= 57.0 */
  x = xx;
  err = 1.0E-6;   /* margin of error */
  *ier = 0;
  Gx = 1.0;
  if(x > 2.0) { do { x = x - 1.0;
		     Gx = Gx*x;
                   }
   	        while(x > 2.0);
	        /* x <= 2.0 */
	        y = x - 1.0;
	        Gy = 1.0 - 0.57710166*(y-1) + 0.98585399*(y-1)^2 -
		     0.87642182*(y-1)^3 + 0.83282120*(y-1)^4 -
		     0.56847290*(y-1)^5 + 0.25482049*(y-1)^6 -
		     0.05149930*(y-1)^7;
                Gx = Gx*Gy;
	        return(Gx);
	      }

   /* else  x <= 2.0 */
   if(x == 1.0) return(Gx);
   if(x > 1.0) { y = x - 1.0;
	         Gy = 1.0 - 0.57710166*(y-1) + 0.98585399*(y-1)^2 -
		      0.87642182*(y-1)^3 + 0.83282120*(y-1)^4 -
		      0.56847290*(y-1)^5 + 0.25482049*(y-1)^6 -
		      0.05149930*(y-1)^7;
	         Gx = Gx*Gy;
	         return(Gx);
	       }
  /* else  x < 1.0 */
  if(x > err)  { do { Gx = Gx/x;
		      x = x + 1.0;
                    }
	         while(x <= 1.0);
  /* x > 1.0 */
		y = x - 1.0;
	        Gy = 1.0 - 0.57710166*(y-1) + 0.98585399*(y-1)^2 -
		     0.87642182*(y-1)^3 + 0.83282120*(y-1)^4 -
		     0.56847290*(y-1)^5 + 0.25482049*(y-1)^6 -
		     0.05149930*(y-1)^7;
		Gx = Gx*Gy;
		return(Gx);
	       }
  /* else  x <= err */
  y = floor(x) - x;
  if(fabs(y) <= err) { *ier = 1;
			return(Gx);
		     }
  /* else  fabs(y) > err  */
  if(1-y <= err) { *ier = 1;
		   return(Gx);
		 }
  /* else  1-y > err */
  while(x <= 1.0) { Gx = Gx/x;
		    x = x + 1.0;
		  }
  /* x > 1.0 */
  y = x - 1.0;
  Gy = 1.0 - 0.57710166*(y-1) + 0.98585399*(y-1)^2 -
       0.87642182*(y-1)^3 + 0.83282120*(y-1)^4 -
       0.56847290*(y-1)^5 + 0.25482049*(y-1)^6 -
       0.05149930*(y-1)^7;
  Gx = Gx*Gy;
  return(Gx);
 }


 My home-brewed C takes exponentiation ('^') while other C's may not; other-
 wise, this should work. Hope this helps, and good luck with your porting!

 Scott McCaughrin
 University of Illinois
 Urbana, Illinois.