/********************************************************************
   Returns the value of 
   Gamma(z) = int_0^infinity t^{z-1} e^{-t} dt      for z > 0
   C.A. Bertulani        May/15/2000
********************************************************************/

typedef double Number;

#include <math.h>
#include<iostream.h>

int main(){
	Number gamma_ln(Number);
	Number fact(int n);
	Number factln(int n);
	Number x, res;
	int j;
	j=1;
	cout << "\n Enter z. To stop, enter z=0.\n";
	cin >>x;
    while(x!=0){
		if(j>10){ cout << "\n My patience is over. Stop, please!\n";
			break;
		}
		if(j!=1){
			cout << "\n Enter z. To stop, enter z=0.\n";
			cin >>x;
		}
		cout << "\n OK. To check things, I will compare with the factorial \n" <<
                 "of the nearest integer of z-1. For z = integer this should be exact.\n";
		res = gamma_ln(x);
		cout << "\n Here is Gamma(z): " << exp(res)
			<< "\n The factorial of the (nearest integer - 1) is:  " << fact(int(x-1)) << endl;
		j=j+1;
	}
	return 0;
}
/********************************************************************
   Returns the value of ln[Gamma(xx)] for xx > 0
********************************************************************/

Number gamma_ln(Number xx)
{
	Number x,y,tmp,ser;
	static Number cof[6]={76.18009172947146,-86.50532032941677,
		24.01409824083091,-1.231739572450155,
		0.1208650973866179e-2,-0.5395239384953e-5};
	int j;

	y=x=xx;
	tmp=x+5.5;
	tmp -= (x+0.5)*log(tmp);
	ser=1.000000000190015;
	for (j=0;j<=5;j++) ser += cof[j]/++y;
	return -tmp+log(2.5066282746310005*ser/x);
}
/********************************************************************
   Returns the factorial of an integer
********************************************************************/
Number fact(int n)
{
	Number gamma_ln(Number xx);
	static int ntop=4;
	static Number a[33]={1.0,1.0,2.0,6.0,24.0}; /* Fill in table only as required */
	int j;

	if (n < 0) cerr << "Negative factorial in routine factrl\n";
	if (n > 32) return exp(gamma_ln(n+1.0));
	while (ntop<n) {		/* Fill in table uo t desired value */
		j=ntop++;
		a[ntop]=a[j]*ntop;
	}
	return a[n];
}
/********************************************************************
   Returns ln(n!)
********************************************************************/
Number factln(int n)
{
	Number gamma_ln(Number xx);
	static Number a[101];  /* a static array is automatically initialized to zero */

	if (n < 0) cerr << "Negative factorial in routine factl\n";
	if (n <= 1) return 0.0;
	if (n <= 100) return a[n] ? a[n] : (a[n]=gamma_ln(n+1.0));   /* In range of table */
	else return gamma_ln(n+1.0);  /* Out of range of table */
}
/********************************************************************/